1105 lines
113 KiB
Plaintext
1105 lines
113 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
|
|
"grade": false,
|
|
"grade_id": "cell-98f724ece1aacb67",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"# CDS: Numerical Methods Assignments\n",
|
|
"\n",
|
|
"- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
|
|
"\n",
|
|
"- Solutions must be submitted via the Jupyter Hub.\n",
|
|
"\n",
|
|
"- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
|
|
"\n",
|
|
"## Submission\n",
|
|
"\n",
|
|
"1. Name all team members in the the cell below\n",
|
|
"2. make sure everything runs as expected\n",
|
|
"3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
|
|
"4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
|
|
"5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
|
|
"6. submit your solutions **in time (before the deadline)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "raw",
|
|
"metadata": {},
|
|
"source": [
|
|
"team_members = \"Koen Vendrig, Kees van Kempen\""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "2900c663e4345c7f0707be39cc0cc7f4",
|
|
"grade": false,
|
|
"grade_id": "cell-b6b5c93e38567117",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Linear Equation Systems\n",
|
|
"\n",
|
|
"In the following you will implement the Gauss-Seidel (GS), Steepest Descent (SD) and the Conjugate Gradient (CG) algorithms to solve linear equation systems of the form \n",
|
|
"\n",
|
|
"$$A \\mathbf{x} = \\mathbf{b},$$ \n",
|
|
"\n",
|
|
"with $A$ being an $n \\times n$ matrix."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "55d777356d474cb7327bb087dfe5e644",
|
|
"grade": true,
|
|
"grade_id": "cell-bcf2dd8f9a194be8",
|
|
"locked": false,
|
|
"points": 0,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import numpy.linalg as linalg\n",
|
|
"from matplotlib import pyplot as plt\n",
|
|
"import time"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "9717dea99ca91d38963fc25b2ef95e03",
|
|
"grade": false,
|
|
"grade_id": "cell-595bb99f65f79ca7",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 1\n",
|
|
"First, you need to implement a Python function $\\text{diff(a,b)}$, which returns the difference $\\text{d}$ between two $n$-dimensional vectors $\\text{a}$ and $\\text{b}$ according to \n",
|
|
"\n",
|
|
"$$ d = || \\mathbf{a} - \\mathbf{b}||_\\infty = \\underset{i=1,2,\\dots,n}{\\operatorname{max}} |a_i - b_i|. $$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "437477eff36ed8c083b4903a6921591c",
|
|
"grade": true,
|
|
"grade_id": "cell-bcd26d8a9f0e6447",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def diff(a, b):\n",
|
|
" return np.max(np.abs(a - b))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "a10d859fa1b26c11ae00c7b846fdd1e8",
|
|
"grade": false,
|
|
"grade_id": "cell-19ae2e6fe5c9264c",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 2 \n",
|
|
"\n",
|
|
"The Gauss-Seidel iteration scheme to solve the linear equation system \n",
|
|
"\n",
|
|
"$$A \\mathbf{x} = \\mathbf{b}$$\n",
|
|
"\n",
|
|
"is defined by \n",
|
|
"\n",
|
|
"$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right].$$\n",
|
|
"\n",
|
|
"Note especially the difference in the sums: the first one involves $x_j^{(k)}$ and the second one $x_j^{(k-1)}$.\n",
|
|
"\n",
|
|
"\n",
|
|
"Give the outline of the derivation in LaTeX math notation in the markdown cell below. (Double click on \"YOUR ANSWER HERE\" to open the cell, and ctrl+enter to compile.) \n",
|
|
"\n",
|
|
"Hint: Similar to the Jacobi scheme, start by seperating the matrix $A$ into its diagonal ($D$), lower triangular ($L$) and upper triangular ($U$) forms, such that $A = D - L - U$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "3c0c1bc83919a390739fa5802fc75743",
|
|
"grade": true,
|
|
"grade_id": "cell-3a8403fee85d9582",
|
|
"locked": false,
|
|
"points": 2,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"---\n",
|
|
"\n",
|
|
"We start from our linear equations:\n",
|
|
"\n",
|
|
"$$Ax = b$$\n",
|
|
"\n",
|
|
"We separate A into different components (diagonal, strictly lower triangular and strictly upper triangular):\n",
|
|
"\n",
|
|
"$$A = D - L - U$$\n",
|
|
"\n",
|
|
"\n",
|
|
"We write $D - L$ as $L'$ to get:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$(L' - U)x = b$$\n",
|
|
"\n",
|
|
"\n",
|
|
"We take the iterative process of the Gauss-Seidel method to write:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$\n",
|
|
"L'x^k = b + Ux^{k-1}\\\\\n",
|
|
"x^k = L'^{-1}(b + Ux^{k-1})\\\\\n",
|
|
"$$\n",
|
|
"If we write every component of the matrix $A$ as $a_{ij}$, we can use forward substitution to rewrite our previous equation to:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$x^k _i = \\frac{1}{a_{ii}}\\left[-\\sum_{j=0}^{i-1}a_{ij}x_{j}^{k} -\\sum_{j=i+1}^{n-1}a_{ij}x_{j}^{k-1} + b_i\\right].$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "dfbc2086b4294a7d942d35d9e1a22d5b",
|
|
"grade": false,
|
|
"grade_id": "cell-4ebd4c774509912d",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 3\n",
|
|
"\n",
|
|
"Implement the Gauss-Seidel iteration scheme derived above\n",
|
|
"\n",
|
|
"$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right],$$\n",
|
|
"\n",
|
|
"where $a_{ij}$ are the elements of the matrix $A$, and $x_i$ and $b_i$ the elements of vectors $\\mathbf{x}$ and $\\mathbf{b}$, respectively.\n",
|
|
"\n",
|
|
"Write a Python function $\\text{GS(A, b, eps)}$, where $\\text{A}$ represents the $n \\times n$ $A$ matrix, $\\text{b}$ represents the $n$-dimensional solution vector $\\mathbf{b}$, and $\\text{eps}$ is a scalar $\\varepsilon$ defining the accuracy up to which the iteration is performed. Your function should return both the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. \n",
|
|
"\n",
|
|
"Use an assertion to make sure the diagonal elements of $A$ are all non-zero. Initialize your iteration with $\\mathbf{x}^{(0)} = \\mathbf{0}$ (or with $\\mathbf{x}^{(1)} = D^{-1}\\mathbf{b}$, with $D$ the diagonal of $A$) and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "65b01c4e4c650dc6971b7d6dde62044b",
|
|
"grade": true,
|
|
"grade_id": "cell-ec2390c3cafa2f8e",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def GS(A, b, eps, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Gauss-Seidel algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" # First we decompose A = D - L - U.\n",
|
|
" D = np.diag(np.diag(A))\n",
|
|
" U = -np.triu(A) + D\n",
|
|
" L = -np.tril(A) + D\n",
|
|
" \n",
|
|
" # We need non-zero diagonals elements.\n",
|
|
" assert np.all(np.diag(D) != 0)\n",
|
|
" \n",
|
|
" x_prev = np.zeros(n)\n",
|
|
" x_cur = np.dot(linalg.inv(D), b)\n",
|
|
" \n",
|
|
" k = 1\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max:\n",
|
|
" k += 1\n",
|
|
" # We will have to copy, as the array elements will point to the same\n",
|
|
" # memory otherwise, and changes to one array will change the other aswell.\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" for i in range(n):\n",
|
|
" x_cur[i] = 1/A[i, i]*(-np.dot(A[i, :i], x_cur[:i]) - np.dot(A[i, i + 1:], x_prev[i + 1:]) + b[i])\n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "badd9e6110882657c17c42e3d9ad8370",
|
|
"grade": false,
|
|
"grade_id": "cell-78d8cad45f5350b9",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 4\n",
|
|
"\n",
|
|
"Verify your implementation by comparing your approximate result to an exact solution. Use $\\text{numpy.linalg.solve()}$ to obtain the exact solution of the system\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{align*}\n",
|
|
" \\begin{pmatrix}\n",
|
|
" 10 & -1 & 2 & 0 \\\\ \n",
|
|
" -1 & 11 &-1 & 3 \\\\\n",
|
|
" 2 & -1 & 10&-1 \\\\\n",
|
|
" 0 & 3 & -1& 8\n",
|
|
" \\end{pmatrix} \\mathbf{x}^*\n",
|
|
" =\n",
|
|
" \\begin{pmatrix}\n",
|
|
" 6 \\\\\n",
|
|
" 25 \\\\\n",
|
|
" -11\\\\\n",
|
|
" 15\n",
|
|
" \\end{pmatrix}\n",
|
|
"\\end{align*}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Then compare you approximate result $\\mathbf{\\tilde{x}}$ to the exact result $\\mathbf{x^*}$ by plotting $|| \\mathbf{x}^* - \\mathbf{\\tilde{x}}||_\\infty$ for different accuracies $\\varepsilon = 10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}$. \n",
|
|
"\n",
|
|
"Implement a unit test for your function using this system."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "b54e0939053f9afee0f0ba46e490ba17",
|
|
"grade": true,
|
|
"grade_id": "cell-0585da22f1d304ab",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQxUlEQVR4nO3df2id133H8c9nipOoGZVL445YiWcXBdGwEAzCgXgEN1mmFKLEhKxLukFLTEwysrF/RG3S9b/WYWJjLQlkLg1uyvLDpEbYiVcxSF2HEDYr0aiTGXUmECI5zMkWqU25UNv97g9dx7IiyffqPI+ee+95v0D4Pue599xvOFw+Oec897mOCAEAkOL3qi4AAND+CBMAQDLCBACQjDABACQjTAAAyQgTAECyy6ouoApXX311bNy4seoyAKCtvPHGGx9GxLrFzmUZJhs3btT4+HjVZQBAW7H97lLnslrmsj1ke+/s7GzVpQBAR8kqTCLiUETs7OnpqboUAOgoWYUJAKAchAkAIFlWG/C2hyQN9fX1VV0KAKyq0YlpjYxN6tRMTevXdmt4sF/bN/cW1n9WMxP2TADkaHRiWrsPHNf0TE0haXqmpt0Hjmt0Yrqw98gqTAAgRyNjk6qdOXdRW+3MOY2MTRb2HoQJAHS4UzO1ptpXgjABgA63fm13U+0rkVWY8KVFADkaHuxX95qui9q613RpeLC/sPfIKkzYgAeQo+2be7Xn3hvVu7ZbltS7tlt77r2x0Ku5sro0GABytX1zb6HhsVBWMxMAQDkIEwBAMsIEAJCMMAEAJMsqTLg0GADKkVWYcGkwAJQjqzABAJSDMAEAJCNMAADJCBMAQDLCBACQjDABACTLKkz4ngkAlCOrMOF7JgBQjqzCBABQDsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAECyrMKE26kAQDmyChNupwIA5cgqTAAA5SBMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJCBMAQLKOCBPbX7L9lO0XbT9SdT0AkJvKw8T207ZP235rQfudtidtn7S9a7k+IuJERDws6auStpZZLwDg0yoPE0n7JN05v8F2l6QnJX1F0g2SHrB9g+0bbb+04O8L9dfcLellSYdXt3wAwGVVFxARR21vXNC8RdLJiHhHkmw/L+meiNgj6a4l+jko6aDtlyU9u/C87Z2SdkrShg0bivsPAABUHyZL6JX03rzjKUk3L/Vk29sk3SvpCi0xM4mIvZL2StLAwEAUVCcAQK0bJk2JiCOSjlRcBgBkqxX2TBYzLem6ecfX1tuS2B6yvXd2dja1KwDAPK0aJsckXW97k+3LJd0v6WBqpxFxKCJ29vT0JBcIALig8jCx/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE21XWCQBYWuV7JhHxwBLth8VlvgDQFiqfmawm9kwAoBxZhQl7JgBQjqzCBABQjqzChGUuAChHVmHCMhdQvtGJaW19/BVt2vWytj7+ikYnkr8ihjZQ+dVcADrH6MS0dh84rtqZc5Kk6Zmadh84Lknavrm3ytJQsqxmJgDKNTI2+UmQnFc7c04jY5MVVYTVQpgAKMypmVpT7egcWYUJG/BAudav7W6qHZ0jqzBhAx4o1/Bgv7rXdF3U1r2mS8OD/RVVhNXCBjyAwpzfZB8Zm9SpmZrWr+3W8GA/m+8ZIEwAFGr75l7CI0NZLXMBAMqRVZiwAQ8A5cgqTNiAB4ByZBUmAIByECYAgGSECQAgGWECAEiWVZhwNRcAlCOrMOFqLgAoR6FhYvuzRfYHAGgPhd1OxfYjks7ZvjUi/rKofgEAra/Ie3Odqf97tsA+AQBtoMhlrvcl9Uo6XWCfAIA2UGSY3CzpNUmbCuwTANAGCguTiPi2pP+T9FBRfQIA2kOhV3NFxJsRMVNkn0XieyYAUI6mN+Btb2jwqTMR8atm+y9TRBySdGhgYIDZEwAUaCVXc/2ogeeEpH2SnllB/wCANtN0mETEl8soBADQvrJa5gIAlINlLgBAsqbCxHYXy1wAgIWWvTTYdq/t79r+TL3p720/UT/3vdKrAwC0hUt9z+RGSTdJGqgffyTp3frjX9s+ZPsqSbI9aPu1csoEALSyZZe5IuKntr8eEUfrTVsk/Vv93Ldsf03SEdu/lfSxpF2lVgsAaElN7ZlExN22PydJtm/X3K1TfiPpGkkPRsRk8SUCAFpdI7dTeXP+QUR8VH/4mKS/i4htku6T9ILt24otDwDQDi45M4mIkSXab5v3+Ljtr0j6iaRbiiuvWLaHJA319fVVXQoAdJQi7xr8vqTbi+qvDPwGPACUo+i7BteK7A8A0B6Sw6S+dAQAyFgRM5PvFNAHAKCNFREmLqAPAEAbKyJMooA+AABtrNANeABAnggTAECyIsLkfwroAwDQxpLDJCLuKKIQAED7YpkLAJCMMAEAJCNMAADJmgqTeb+q+PvllAMAaEfNzkw+Z/tRSX9cRjEAgPbUbJjcLukbkr5o+wvFl7Nytq+yPW77rqprAYDcNBsm/yHpQUnvRsTpIgqw/bTt07bfWtB+p+1J2ydtN/Lb8t+UtL+ImgAAzWn2N+BP1B/+osAa9kl6QtIz5xtsd0l6UtIdkqYkHbN9UFKXpD0LXv+gpJsk/ZekKwusCwDQoKbCpAwRcdT2xgXNWySdjIh3JMn285LuiYg9kj61jGV7m6SrJN0gqWb7cET8rsy6AQAXNLzMZbvX9ndtf6Z+/D3bZd1+vlfSe/OOp+pti4qIxyLibyU9K+kHiwWJ7Z31PZXxDz74oOh6ASBrzeyZ3Ki55aSB+vGvJR2cd7nwoO3XCq6vKRGxLyJeWuLc3ogYiIiBdevWrXZpANDRGl7mioif2v56RBytH3/L9tckHbH9W0kfS2pko7wR05Kum3d8bb0NANCCVvwNeNu3S3pI0m8kXS3pbyLi1YLqOibpetubbF8u6X5JB1M7tT1ke+/s7GxygQCAC5oNkzfnPX5M0rcjYpuk+yS9YPu2Zguw/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE2832vVBEHIqInT09PaldAQDmcUQxv7pr+xpJP4mIWwrpsEQDAwMxPj5edRkA0FZsvxERA4uda/rSYNsbljm9Y975mYj4VbP9l8n2kKShvr6+qksBgI7S9MzE9s+WOR2SXP93X0Q8s8xzK8PMBACaV+jMJCK+nF4SAKCTFL3MNV/LLXMBAMqxktup/KiB54Tm7rnVUstc7JkAQDkKu5qrnbBn0jlGJ6Y1MjapUzM1rV/breHBfm3fvOSddwAkKHTPBGgVoxPT2n3guGpnzkmSpmdq2n3guCQRKMAq4zfg0bZGxiY/CZLzamfOaWRssqKKgHxlFSbcTqWznJqpNdUOoDxZhQm3U+ks69d2N9UOoDxZhQk6y/Bgv7rXdF3U1r2mS8OD/RVVBOSLDXi0rfOb7FzNBVSPMEFb2765l/AAWkBWy1xswANAObIKEzbgAaAcWYUJAKAchAkAIBlhAgBIRpgAAJIRJgCAZFmFCZcGA0A5sgoTLg0GgHJkFSYAgHIQJgCAZIQJACAZYQIASEaYAACSESYAgGRZhQnfMwGAcmQVJnzPBADKkVWYAADKQZgAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEiWVZhwOxUAKEdWYcLtVACgHFmFCQCgHIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASNYRYWJ7m+1XbT9le1vV9QBAbioPE9tP2z5t+60F7XfanrR90vauS3QTkj6WdKWkqbJqBQAs7rKqC5C0T9ITkp4532C7S9KTku7QXDgcs31QUpekPQte/6CkVyPi57b/QNI/SvqLVagbAFBXeZhExFHbGxc0b5F0MiLekSTbz0u6JyL2SLprme4+knTFYids75S0U5I2bNiQWjYAYJ7Kl7mW0CvpvXnHU/W2Rdm+1/Y/S/qx5mY5nxIReyNiICIG1q1bV2ixAJC7ymcmRYiIA5IOVF0HAOSqVWcm05Kum3d8bb0tie0h23tnZ2dTuwIAzNOqYXJM0vW2N9m+XNL9kg6mdhoRhyJiZ09PT3KBAIALKg8T289Jel1Sv+0p2zsi4qykRyWNSTohaX9EvF1lnQCApVW+ZxIRDyzRfljS4SLfy/aQpKG+vr4iuwWA7FU+M1lNLHMBQDmyChMAQDkIEwBAsqzChEuDAaAcWYUJeyYAUI6swgQAUA7CBACQjDABACTLKkzYgAeAcmQVJqkb8KMT09r6+CvatOtlbX38FY1OJN97EgA6QuW3U2kXoxPT2n3guGpnzkmSpmdq2n3guCRp++Ylf2oFALKQ1cwkxcjY5CdBcl7tzDmNjE1WVBEAtA7CpEGnZmpNtQNATrIKk5QN+PVru5tqB4CcZBUmKRvww4P96l7TdVFb95ouDQ/2F1UeALQtNuAbdH6TfWRsUqdmalq/tlvDg/1svgOACJOmbN/cS3gAwCKyWuYCAJSDMAEAJMsqTLidCgCUI6sw4fdMAKAcWYUJAKAcjoiqa1h1tj+Q9O4ip3okLVwDW6ztakkfllDapSxWy2r108hrLvWc5c4vda7Vx0QqZlzKGpNGnlfWuLT7mKy0n07+rPxhRKxb9ExE8Ff/k7S3wbbxVqlvtfpp5DWXes5y55c61+pjUtS4lDUmVY5Lu49JmePSiZ8VlrkudqjBtqoUVctK+mnkNZd6znLnlzrX6mMiFVNPWWPSyPM6cVz4rDReSyGyXOZKZXs8IgaqrgMXMCathzFpTWWNCzOTldlbdQH4FMak9TAmramUcWFmAgBIxswEAJCMMAEAJCNMAADJCJOC2b7K9rjtu6quBXNsf8n2U7ZftP1I1fVAsr3d9g9sv2D7T6uuB3Nsf9H2D22/2OxrCZM620/bPm37rQXtd9qetH3S9q4GuvqmpP3lVJmfIsYlIk5ExMOSvippa5n15qCgMRmNiIckPSzpz8usNxcFjcs7EbFjRe/P1VxzbN8q6WNJz0TEH9XbuiT9UtIdkqYkHZP0gKQuSXsWdPGgpJskfV7SlZI+jIiXVqf6zlXEuETEadt3S3pE0o8j4tnVqr8TFTUm9df9g6R/iYg3V6n8jlXwuLwYEfc18/780mJdRBy1vXFB8xZJJyPiHUmy/bykeyJij6RPLWPZ3ibpKkk3SKrZPhwRvyuz7k5XxLjU+zko6aDtlyURJgkK+qxY0uOS/pUgKUZRn5WVIkyW1yvpvXnHU5JuXurJEfGYJNn+huZmJgRJOZoal3rI3yvpCkmHyywsY02NiaS/lvQnknps90XEU2UWl7FmPyufl/QdSZtt766HTkMIkxJExL6qa8AFEXFE0pGKy8A8EfF9Sd+vug5cLCL+V3P7WE1jA35505Kum3d8bb0N1WJcWg9j0ppWbVwIk+Udk3S97U22L5d0v6SDFdcExqUVMSatadXGhTCps/2cpNcl9duesr0jIs5KelTSmKQTkvZHxNtV1pkbxqX1MCatqepx4dJgAEAyZiYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJkCLqP+I0X/W//7dNp9PtA2+AQ+0CNv/LenWiHi/6lqAZvF/PkDrOCzpF7b/qepCgGbxeyZAC7B9iyRLuqZ+cz6grTAzAVrDn0n6ZUSc9ZzPVl0Q0Az2TIAWYHuLpB9KCkk1SX8VEW9UWxXQOMIEAJCMZS4AQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMn+H8is0a3Xr+ppAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"A = np.array([[ 10, - 1, 2, 0],\n",
|
|
" [- 1, 11, - 1, 3],\n",
|
|
" [ 2, - 1, 10, - 1],\n",
|
|
" [ 0, 3, - 1, 8]])\n",
|
|
"b = np.array( [ 6, 25, -11, 15] )\n",
|
|
"x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
"eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n",
|
|
"diff_list = []\n",
|
|
"for eps in eps_list:\n",
|
|
" x, k = GS(A, b, eps)\n",
|
|
" diff_list.append(diff(x_exact, x))\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"\n",
|
|
"ax.scatter(eps_list, diff_list)\n",
|
|
"ax.set_xscale(\"log\")\n",
|
|
"ax.set_yscale(\"log\")\n",
|
|
"ax.set_xlabel(\"$\\epsilon$\")\n",
|
|
"ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n",
|
|
"\n",
|
|
"fig.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# As the three algorithm functions will have the same signature,\n",
|
|
"# it makes sense to only write the test function once.\n",
|
|
"\n",
|
|
"def test_alg(alg, alg_name):\n",
|
|
" \"\"\"\n",
|
|
" Check that function alg returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" A = np.array([[ 10, - 1, 2, 0],\n",
|
|
" [- 1, 11, - 1, 3],\n",
|
|
" [ 2, - 1, 10, - 1],\n",
|
|
" [ 0, 3, - 1, 8]])\n",
|
|
" b = np.array( [ 6, 25, -11, 15] )\n",
|
|
" x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
" print(\"Starting with A =\")\n",
|
|
" print(A)\n",
|
|
" print(\"and b =\", b)\n",
|
|
" print(\"We apply the {} algorithm to solve Ax = b.\".format(alg_name))\n",
|
|
" print()\n",
|
|
"\n",
|
|
" eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations:\".format(eps, k))\n",
|
|
" print(\"x =\\t\\t\\t\", x)\n",
|
|
" print(\"Ax =\\t\\t\\t\", np.dot(A, x))\n",
|
|
" print(\"diff(Ax, b) =\\t\\t\", diff(A @ x, b))\n",
|
|
" print(\"diff(x, x_exact) =\\t\", diff(x, x_exact))\n",
|
|
" print()\n",
|
|
" \n",
|
|
" assert diff(x, x_exact) < eps\n",
|
|
" "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "d8e54a45261b013b9dd07a401d86401c",
|
|
"grade": true,
|
|
"grade_id": "cell-6524fb6d322a6ea0",
|
|
"locked": false,
|
|
"points": 0,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Gauss-Seidel algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 0.99463393 1.99776509 -0.99803257 1.00108402]\n",
|
|
"Ax =\t\t\t [ 5.95250909 24.98206671 -10.98990699 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.04749090616931895\n",
|
|
"diff(x, x_exact) =\t 0.005366066491359844\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 0.99938302 1.99982713 -0.99978549 1.00009164]\n",
|
|
"Ax =\t\t\t [ 5.99443213 24.99877578 -10.99900762 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.005567865937722516\n",
|
|
"diff(x, x_exact) =\t 0.000616975874427883\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 6\t iterations:\n",
|
|
"x =\t\t\t [ 0.99993981 1.99998904 -0.99997989 1.00000662]\n",
|
|
"Ax =\t\t\t [ 5.99944928 24.99993935 -10.99991498 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.000550717702960668\n",
|
|
"diff(x, x_exact) =\t 6.018928065554263e-05\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 7\t iterations:\n",
|
|
"x =\t\t\t [ 0.99999488 1.99999956 -0.99999836 1.00000037]\n",
|
|
"Ax =\t\t\t [ 5.99995255 24.99999971 -10.99999375 15. ]\n",
|
|
"diff(Ax, b) =\t\t 4.744782363452771e-05\n",
|
|
"diff(x, x_exact) =\t 5.11751035947583e-06\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_GS():\n",
|
|
" \"\"\"\n",
|
|
" Check that GS returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" return test_alg(GS, \"Gauss-Seidel\")\n",
|
|
" \n",
|
|
"test_GS()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "f88ad91685d2bd73471acc5b600b5177",
|
|
"grade": false,
|
|
"grade_id": "cell-59514f17e337138f",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 5\n",
|
|
"\n",
|
|
"Next, implement the Steepest Descent algorithm in a similar Python function $\\text{SD(A, b, eps)}$, which calculates\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\mathbf{v}^{(k)} &= \\mathbf{b} - A \\mathbf{x}^{(k-1)} \\\\\n",
|
|
" t_k &= \\frac{ \\langle \\mathbf{v}^{(k)}, \\mathbf{v}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)}\\rangle } \\\\\n",
|
|
" \\mathbf{x}^{(k)} &= \\mathbf{x}^{(k-1)} + t_k \\mathbf{v}^{(k)} .\n",
|
|
"\\end{align*}\n",
|
|
"\n",
|
|
"Initialize your iteration again with $\\mathbf{x}^{(0)} = \\mathbf{0}$ and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n",
|
|
"Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products. "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "215218647aa6a6a314c093e82a4a24fb",
|
|
"grade": true,
|
|
"grade_id": "cell-a105ce7e1a34ce3c",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def SD(A, b, eps, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Steepest Descent algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" \n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" x_cur = np.zeros(n)\n",
|
|
" x_prev = np.zeros(n)\n",
|
|
" \n",
|
|
" k = 0\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
|
|
" k += 1\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" \n",
|
|
" v = b - A @ x_prev\n",
|
|
" t = np.dot(v, v)/np.dot(v, A @ v)\n",
|
|
" x_cur = x_prev.copy() + t*v\n",
|
|
" \n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "f554d4f402c9e1a4cbff1f4e9018e9d0",
|
|
"grade": true,
|
|
"grade_id": "cell-62100bf55800145b",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Steepest Descent algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 0.99748613 1.98300329 -0.98904751 1.01283183]\n",
|
|
"Ax =\t\t\t [ 6.01376302 24.84309304 -10.89133793 15.04071202]\n",
|
|
"diff(Ax, b) =\t\t 0.15690696195356324\n",
|
|
"diff(x, x_exact) =\t 0.016996711694861055\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 6\t iterations:\n",
|
|
"x =\t\t\t [ 0.99983175 1.99716093 -0.99850509 1.00217552]\n",
|
|
"Ax =\t\t\t [ 6.00414638 24.97397012 -10.98472385 15.00739201]\n",
|
|
"diff(Ax, b) =\t\t 0.02602988052923294\n",
|
|
"diff(x, x_exact) =\t 0.002839069878101119\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 9\t iterations:\n",
|
|
"x =\t\t\t [ 0.99991029 1.99994247 -0.99999645 1.00023877]\n",
|
|
"Ax =\t\t\t [ 5.99916754 25.00016961 -11.00032515 15.00173404]\n",
|
|
"diff(Ax, b) =\t\t 0.0017340408511579142\n",
|
|
"diff(x, x_exact) =\t 0.00023877424067331177\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 11\t iterations:\n",
|
|
"x =\t\t\t [ 0.99998551 1.99999065 -0.99999949 1.00003874]\n",
|
|
"Ax =\t\t\t [ 5.9998655 25.00002733 -11.00005329 15.00028137]\n",
|
|
"diff(Ax, b) =\t\t 0.0002813662634846281\n",
|
|
"diff(x, x_exact) =\t 3.874139053650083e-05\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_SD():\n",
|
|
" \"\"\"\n",
|
|
" Check that SD returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" return test_alg(SD, \"Steepest Descent\")\n",
|
|
" \n",
|
|
"test_SD()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "695745368fdf84d75691832839fb2834",
|
|
"grade": false,
|
|
"grade_id": "cell-3277aa3696a5c87f",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 6\n",
|
|
"\n",
|
|
"Finally, based on your steepest decent implementation from task 5, implement the Conjugate Gradient algorithm in a Python function $\\text{CG(A, b, eps)}$ in the following way: \n",
|
|
"\n",
|
|
"Initialize your procedure with:\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" \\mathbf{x}^{(0)} &= \\mathbf{0} \\\\\n",
|
|
" \\mathbf{r}^{(0)} &= \\mathbf{b} - A \\mathbf{x}^{(0)} \\\\\n",
|
|
" \\mathbf{v}^{(0)} &= \\mathbf{r}^{(0)}\n",
|
|
"\\end{align*}\n",
|
|
"\n",
|
|
"Then increase $k$ and repeat the following until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$.\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
" t_k &= \\frac{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)} \\rangle } \\\\\n",
|
|
" \\mathbf{x}^{(k+1)} &= \\mathbf{x}^{(k)} + t_k \\mathbf{v}^{(k)} \\\\\n",
|
|
" \\mathbf{r}^{(k+1)} &= \\mathbf{r}^{(k)} - t_k A \\mathbf{v}^{(k)} \\\\\n",
|
|
" s_k &= \\frac{ \\langle \\mathbf{r}^{(k+1)}, \\mathbf{r}^{(k+1)} \\rangle }{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle } \\\\\n",
|
|
" \\mathbf{v}^{(k+1)} &= \\mathbf{r}^{(k+1)} + s_k \\mathbf{v}^{(k)}\n",
|
|
"\\end{align*}\n",
|
|
"\n",
|
|
"Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n",
|
|
"Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products.\n",
|
|
"\n",
|
|
"How do you expect the number of needed iteration steps to behave when changing the accuracy $\\epsilon$? What do you see?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "f709b7bfc8af2ba2d607f60ed03e1074",
|
|
"grade": true,
|
|
"grade_id": "cell-8ac3d07c9d237e47",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def CG(A, b, eps, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Conjugate Gradient algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" \n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" x_cur = np.zeros(n)\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" r_cur = b - A @ x_cur\n",
|
|
" v = r_cur\n",
|
|
" \n",
|
|
" k = 0\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
|
|
" k += 1\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" r_prev = r_cur\n",
|
|
" \n",
|
|
" t = np.dot(r_prev, r_prev)/np.dot(v, A @ v)\n",
|
|
" x_cur = x_prev + t*v\n",
|
|
" r_cur = r_prev - t*A @ v\n",
|
|
" s = np.dot(r_cur, r_cur)/np.dot(r_prev, r_prev)\n",
|
|
" v = r_cur + s*v\n",
|
|
" \n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "0d5b479bee6c0297974c003b4ebd99f3",
|
|
"grade": true,
|
|
"grade_id": "cell-4f5bc81f40ddb3fa",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Conjugate Gradient algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_CG():\n",
|
|
" \"\"\"\n",
|
|
" Check that CG returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # TODO: Do we need to plot the error over epsilon like in task 4?\n",
|
|
" \n",
|
|
" return test_alg(CG, \"Conjugate Gradient\")\n",
|
|
" \n",
|
|
"test_CG()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "65dd782269d01314b02643d0610e3348",
|
|
"grade": false,
|
|
"grade_id": "cell-308595a088486388",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 7\n",
|
|
"\n",
|
|
"Apply all three methods to the following system\n",
|
|
"\n",
|
|
"\\begin{align*}\n",
|
|
"\\begin{pmatrix}\n",
|
|
"0.2& 0.1& 1.0& 1.0& 0.0 \\\\ \n",
|
|
"0.1& 4.0& -1.0& 1.0& -1.0 \\\\\n",
|
|
"1.0& -1.0& 60.0& 0.0& -2.0 \\\\\n",
|
|
"1.0& 1.0& 0.0& 8.0& 4.0 \\\\\n",
|
|
"0.0& -1.0& -2.0& 4.0& 700.0\n",
|
|
"\\end{pmatrix} \\mathbf{x}^*\n",
|
|
"=\n",
|
|
"\\begin{pmatrix}\n",
|
|
"1 \\\\\n",
|
|
"2 \\\\\n",
|
|
"3 \\\\\n",
|
|
"4 \\\\\n",
|
|
"5\n",
|
|
"\\end{pmatrix}.\n",
|
|
"\\end{align*}\n",
|
|
" \n",
|
|
"Plot the number of needed iterations for each method as a function of $\\varepsilon$, using $\\varepsilon = 10^{-1}, 10^{-2}, ..., 10^{-8}$.\n",
|
|
"\n",
|
|
"Explain the observed behavior with the help of the condition number (which you can calculate using $\\text{numpy.linalg.cond()}$). "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "3d94dc67b052b82dae655b7ff562c7d2",
|
|
"grade": true,
|
|
"grade_id": "cell-490cd96dc58cbff8",
|
|
"locked": false,
|
|
"points": 4,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA93UlEQVR4nO3dd3hUVf7H8fdJJyEJhCQk9NBJA0LoAhFEUBEUUURUsGFZbPuzgLq2tcuuLpZVdkVEEVRAhQUVFTGggNIJNYAgoaQBCenJzPn9ccOQCglk5s4k39fzzJOZOyWfGcL5zj3n3nOU1hohhBCiLDezAwghhHA+UhyEEEJUIsVBCCFEJVIchBBCVCLFQQghRCVSHIQQQlQixUEIIUQlUhyEEEJU4mF2gIqUUgnA34EdwAKt9arzPSc4OFi3a9fOrrmEEKK+2bhxY4bWOqSq+xxSHJRSs4FRQJrWOrrM9pHAvwB34L9a61cADeQAPkBKTV6/Xbt2bNiwoc5zCyFEfaaUOlTdfY7qVpoDjCy7QSnlDrwDXAFEAhOUUpHAaq31FcDjwHMOyieEEKIMhxQHrXUicKLC5j7APq31Aa11EbAAGKO1tpbefxLwdkQ+IYQQ5Zk55tASOFzmdgrQVyk1FhgBNAHeru7JSqkpwBSANm3a2C+lEEI0QE43IK21XgwsrsHjZimljgFXe3l59bJ/MiGEaDjMPJT1CNC6zO1WpdtqTGu9VGs9JTAwsE6DCSFEQ2dmcfgd6KSUilBKeQE3AktMzCOEEKKUQ4qDUmo+sBboopRKUUrdobUuAaYC3wG7gM+11jtq+bpXK6VmZWVl1X1oIYRwRkV5cPIQpGyA3cvh1OHzP+cCqPqwElx8fLyW8xyEEC5JayjIgtx045KTVvl62W1FOeWfP/otiLv1gn61Umqj1jq+qvucbkC6NpRSVwNXd+zY0ewoQghxltUCeZnVN/YVr1uKqngRBb7NwC8EGodAy15nr/uFnr0e1N4ub0H2HIQQoiZKCmve2Odlgu2UrTLcPKtu4Ku67tsM3O37/b3e7jkIIcRFKzwNWUcgO8X4efpYFQ1/OhRWM7bp6Xe2UW8aAa37VN/w+zQBpRz69i6USxcH6VYSQpxTcT5kH4WsFMg+Ur4InLldVaPfqKnRqDcOhbCYc3zLDwEvP8e/LweQbiUhhGuyFBsNf5WNfmkxyMus/DzfYAhsCQGtSn+2hMBWpT9bQuMw8PBy/PsxgXQrCSFci9UCOanVN/pZR4z7qfDl1jvQaOgDWxoDuBWLQEBL8PQx5S25GpcuDtKtJIQL0tr4Rl9VV8+ZbaePgbWk/PM8fc9+u+94WZlv/WUKgLe/Oe+pHpJuJSFE3SvKhYxkyNhr/Mw6fLbhzz4KJQXlH+/uBQEtKnT1VPjW36ipywzmugrpVhJC1D2tjSN6MvZWuJQWgzOUG/iHGw18eHfociUEti7f3+8bDG6yarEzkeIghDg3SwmcOgTpe842/hml1wvKHOnj6QfBnaBNfwieBCGdIbizcZKWhyzN4mpcujjImIMQdagwBzKTjca/bCE4sb/8GbyNmxuNfvQ44+eZIuDfQr791yMy5iBEQ1KuK2hPmUKQbAwKn6HcISjCaPSDO0Fwl7PXGzUxLb6oWzLmIERDYymBkwerGA+opiuo3cDSItDZKARBEdIV1MBJcRDClZ3pCkqvUAAy94O1+OzjGocZjX/0OAjpcrYQBLSUI4BElaQ4COEqigsg5Tc48DMc2WgUgewyiyfauoK6QOeRpeMBXaBZR+kKErXm0sVBBqRFvWa1wLEtRjH442f4c51xfoByh7BoaHdJaTdQ2aOCGsa0D8L+ZEBaCGehtTE4/MfPRkE4uObspHChURAxGNoPgbYDwSfA3KyiXpABaSGc1anDZ4vBH4mQc9zY3qQtRI2BiCFGUWgcam5O0eBIcRDCkXIzjCJwpiCc/MPY7hdiFIGIIcbeQdN2psYUQoqDEPZUmAOHfj1bDFK3G9u9/I0xg753GwUhtJscNSScihQHIepSSSGk/H52EPnIRmN2UXdvY4WwoU9BRAK06Gn3JSCFuBjy1ynExbBa4Pi2s8Xg0FooyTcmm2vREwY8YHQTte4Lno3MTitEjbl0cZBDWYXDaW1MNfHHz3BglXFEUcEp476QrhB369kjiuTcAuHC5FBWIc4n60j5I4pOHzW2B7aB9oONbqKIweDf3MyUQtSaHMoqRG3knYCDq892FWXuM7b7NqtwRFGEDCKLekuKgxBgHFG0Z7lREI5vBzR4NTa6h+JvLz2iKFKmpBYNhhQHIX75F3z/tLFUZas+cOkTRjFoGQfunmanE8IUUhxEw6U1rHwBVs+AqLEw5m3w8jM7lRBOQYqDaJisVvh2Gvz2vnGE0ag3wc3d7FRCOA0pDqLhsZTAkvth66fQfypc/oIMLAtRgRQH0bCUFMKiO2DXUrj0SRj8qBQGIarglIdeKKX8lFIblFKjzM4i6pGiXJh/o1EYRr4CQx6TwiBENRxSHJRSs5VSaUqppArbRyql9iil9imlppW563Hgc0dkEw1E/in4eKxxVvOYd6DfvWYnEsKpOWrPYQ4wsuwGpZQ78A5wBRAJTFBKRSqlhgM7gTQHZRP1XU46fDTKmARv3IfQ82azEwnh9Bwy5qC1TlRKtauwuQ+wT2t9AEAptQAYAzQG/DAKRr5SarnW2lrxNZVSU4ApAG3atLFjeuHSso7A3DGQlQITFkCny8xOJIRLMHNAuiVwuMztFKCv1noqgFJqMpBRVWEA0FrPAmaBMbeSfaMKl5S5H+ZeY0yMd8tiaDvA7ERCuAynPVpJaz3nfI+RWVlFtVJ3wsfXGGspTFoKLXqYnUgIl2Lm0UpHgNZlbrcq3VZjWuulWuspgYGBdRpMuLiUjTDnSmNNhdu+kcIgxAUwszj8DnRSSkUopbyAG4EltXkBpdTVSqlZWVlZdgkoXNAfq2HuaPAJhNu/hZAuZicSwiU56lDW+cBaoItSKkUpdYfWugSYCnwH7AI+11rvqM3ryp6DKGfPt/DJdRDYCm77Fpq2MzuREC7LUUcrTahm+3JguSMyiHpu+0L48m4Ii4GJi8CvmdmJhHBpTnmGdE1Jt5IAYMOHsOhOY53mW5dIYRCiDrh0cZBuJcEvM+F/D0Gn4XDzIvAJMDuREPWCSxcH2XNowM6sxfD93yDqWhg/DzwbmZ1KiHrDpYuD7Dk0UFYrfPM4JL5urMVw3Qfg4WV2KiHqFac9CU6IKllKYOkDsGWerMUghB1JcRCuo6TQGHjetQQSnpApt4WwI5cuDjJ9RgNSlAuf3Qz7VxprMciU20LYlYw5COdXkHV2LYbRb0thEMIBXHrPQTQAuRnwyVhjIr1xs40jk4QQdufSew5yKGs9l30UPrwC0vcaazFIYRDCYVy6OEi3Uj124gDMHgHZx4y1GGSRHiEcSrqVhPNJ3QkfXwuWIpi0BFrGmZ1IiAbHpfccRD10pHQtBjDWYpDCIIQppDgI5/HHavhoNHgHGGsxhHY1O5EQDZYUB+Ec9n4H88YZazHc/i0ERZidSIgGzaWLgxytVE8kLYIFN0FIV5i8HAJamJ1IiAbPpYuDHK1UD2ycAwvvMNZimLRU1mIQwkm4dHEQLu7Xt2Dpg9DxMpi4UNZiEMKJyKGswvG0hp9egsTXIPIaGPsfmXJbCCcjxUE4ltUK302H9e9Bz1vg6n+Bm7vZqYQQFUhxEI5Tdi2Gfn+BES/KlNtCOCkpDsIxSopg8Z2w82tImA5DHpfCIIQTc+niIOs5uIiiPPj8Ftj3A4x4Cfr/xexEQojzcOmjleRQVhdQkAWfXAf7foTRb0lhEMJFuPSeg3ByuZnwybWQusNYiyF6rNmJhBA1JMVB2Ef2MZg7Bk4dghvnQ+fLzU4khKgFKQ7CPpY/AlkpcPMiaHeJ2WmEELXk0mMOwkkd2QS7/weXPCSFQQgXJcVB1L2fXoRGQdD3HrOTCCEukBQHUbcOrTUOWb3kIZkrSQgX5nTFQSnVTSn1nlJqoVLqXrPziFrQGla+AH6h0Psus9MIIS6CQ4qDUmq2UipNKZVUYftIpdQepdQ+pdQ0AK31Lq31PcANwEBH5BN15I+f4dAaGPwIePmanUYIcREctecwBxhZdoNSyh14B7gCiAQmKKUiS+8bDSwDljson7hYZ/YaAlpBr8lmpxFCXCSHFAetdSJwosLmPsA+rfUBrXURsAAYU/r4JVrrK4CJjsgn6kDyCkj53dhr8PA2O40Q4iKZeZ5DS+BwmdspQF+lVAIwFvDmHHsOSqkpwBSANm3a2C2kqIEzew1N20HPm81OI4SoA053EpzWehWwqgaPmwXMAoiPj9f2TSXOaddSOL4NrnkP3D3NTiOEqANmHq10BGhd5nar0m01ppS6Wik1Kysrq06DiVqwWoxV3YI7Q+wNZqcRQtQRM4vD70AnpVSEUsoLuBFYUpsXkFlZnUDSYkjfZazRICu6CVFvOOpQ1vnAWqCLUipFKXWH1roEmAp8B+wCPtda76jl68qeg5ksJbDqJWgebawFLYSoNxwy5qC1nlDN9uVcxOGqWuulwNL4+Hg548oMW+fDiQPGrKtuTnc+pRDiIrj0/2jZczBRSSH8/Cq0iIMuV5idRghRx1y6OMiYg4k2zYWswzD0KVkLWoh6yKWLgzBJcT4kzoA2A6DDULPTCCHswKWLg3QrmeT3DyDnOAx9UvYahKinXLo4SLeSCQpzYM0b0D5BFvIRoh5z6eIgTPDb+5CXAZc+ZXYSIYQduXRxkG4lB8s/Bb/8CzqPhNa9zU4jhLAjly4O0q3kYOvehYIsuPQJs5MIIezMpYuDcKDcTFj7LkSOgfDuZqcRQtiZFAdRM7/+C4pyIEH2GoRoCFy6OMiYg4OcToX1s4xZV0O7mp1GCOEALl0cZMzBQdb8EyxFMORxs5MIIRzEpYuDcICsFNgwG3rcBM06mJ1GCOEgUhzEuSXOMJYBHfKY2UmEEA4kxUFU78QfsPlj6DUZmsg63UI0JC5dHGRA2s5+fg3cPGDQ/5mdRAjhYDUuDkqpF6rYZuq6kDIgbUfpe2HbAuh9JwSEm51GCOFgtdlzaKmUsq3oppQKBX6o+0jCKax6GTwawSUPm51ECGGC2iwTejfwnVJqP6CBDwE5trE+Op4EOxbDoEfAL9jsNEIIE5y3OCil5gKbgM3AX4BPgRLgGq31PvvGE6b46SXwDoQBU81OIoQwSU26leYACrgN+ARoB5wEblZKjbNbMmGOIxthzzIYcD80amp2GiGESc6756C1XgmsPHNbKeUBdAO6A32BhXZLJxxv5YvQKAj63WN2EiGEiWoz5gCA1roE2F56+aTOE9WCUupq4OqOHTuaGaP+OLQW9v8Iw58Hb3+z0wghTOTS5znIoax1SGtY+QL4hULvu8xOI4QwmUsXB1GH/vgZDq2BwY+Al6/ZaYQQJpPiIM7uNQS0MqbKEEI0eFIcBCSvgJTfYcij4OFtdhohhBOQ4tDQWa3GXkPTdtBjotlphBBOotZHK4l6ZvdSOL4Nrn0f3D3NTiOEcBKy59CQWS3G2dDBnSHmerPTCCGciOw5NGRJiyB9N1w/B9xMnWBXCOFknLI4KKWuAa4CAoAPtNYrzE1UD1lKjJlXm8dAtzFmpxFCOBmHdSsppWYrpdKUUkkVto9USu1RSu1TSk0D0Fp/pbW+C7gHGO+ojA3K1vlw4gBc+gS4Se+iEKI8R7YKc4CRZTeULhb0DnAFEAlMUEpFlnnIU6X3i7pUUgg/vwot4qDLFWanEUI4IYcVB611InCiwuY+wD6t9QGtdRGwABijDK8C32itNzkqY4OxaS5kHYahT4FSZqcRQjghs/sTWgKHy9xOKd12P3AZME4pVeX0oEqpKUqpDUqpDenp6fZPWl8U50PiDGgzADoMNTuNEMJJOeWAtNZ6JjDzPI+ZBcwCiI+P147IVS/8/gHkHIdxs2WvQQhRLbP3HI4ArcvcblW6rUaUUlcrpWZlZWXVebB6qTAH1vwT2l8K7QaanUYI4cTMLg6/A52UUhFKKS/gRmBJTZ8sU3bX0vr3IC/TGGsQQohzcOShrPOBtUAXpVSKUuqO0oWDpgLfAbuAz7XWO2rxmrLnUFP5p+DXmdD5CmgVb3YaIYSTU1q7fnd9fHy83rBhg9kxnNtPLxmHr969GsJjzU4jhHACSqmNWusqvy2a3a10UWTPoYZyM2HtuxB5jRQGIUSNuHRxkDGHGvr1X1CUAwnTzU4ihHARTnkoq6hDp1Nh/SyIvQFCu5qdRghRBatVk1tUQm6hhZzCEnIKS8gt89O4brFtK3v/lMHtGdQppM4zuXRxUEpdDVzdsWNHs6M4rzX/BEsRDHnc7CRC1CuFJRZyCs426LlF5Rvz06X3VdxuXC9fBPKKLDX6nW4KGnt70NjbA7/SS1GJ1S7vz6WLg9Z6KbA0Pj7+LrOzOKWsFNgwG3pOhGYdzE4jhKlKLFajUS4621DbGvbCEvKKyn87zy1t8I1GvnwRyC0sodhSs4N5vD3c8Pcpbcy9jIY9uLEXbZv5Gtu9jPvONvju5QpA2Z8+nm4oB5286tLFQZxH4uvGz8GPmZtDiAtQbLHaGuq8Iku5b9+Vv5FbbI15TqGFvDMNfJmumpp+w3ZTGI24z9lv54293Qlu7G002j5lGm0v9/INuE+Z615GQ+/h7ppDuy5dHKRb6RxOHIDNn0D87dCk9fkfL0Qd0FqTW2QhK7+YrLxiTuUXkZ1fvo88r6j8N/bcIku5LpczhaBWjbn3mW/g7rbGuZWvL429yzfevl7u5b6J+1Zo3P283Wnk6e6wb+cXIy0vje3p24kKjiLML6zOX9+li4N0K53Dz6+BmwcM+j+zkwgXo7Umv9ho4E/lFRsNfeklu/RnVdtPlf4ssZ67u+VMY16239zPy50gP9/Sbe62rpYz39p9vSp3u5zZ5siuFrMUlBSw68QutqVvMy4Z2zieexyAp/o+xfiudb/sjUsXB1GN9L2w7TPo/xfwr/tvFMI1FJQ28FU15ll5ReVvl2ncs/KLz9mfrhQE+HjSxNeTwEbGpWXTRgQ28qRJo7PbzlwCGnna+twbe3vg7VH/G/OLobXmUPYhtmVssxWD5JPJlOgSAFo2bknPkJ7ERMYQGxJL1yD7HIUoxaE+WvUyePrCwIfMTiLqQEGxhRO5ReUuJ0sb91N5Zxv0so18Vn7xebtlAnw8CCzTwIcF+hDYyKtS4162CAQ08sTf2wM3N2nc60pWYRbbM7bb9gi2p28nuygbAD9PP6KDo7kt+jZigmOICYkhuFGwQ3K5dHGQMYcqHE+CHYth0CPg55g/IlFzVqsmu6C4UmOfmVvEyTO388rfd67DHP29PQgo05B3DG1criGv2Lifufj7eOIuDbzDFVuL2XtyL9vSjSKwLWMbh7IPAeCm3OjQpAPD2w4nNiSW2OBYIgIjcHdzNyWrzK1U38y/CQ6ugYe2QqOmZqep9wpLLJzMLdPY5xVxIqewykbe+MZfjKWaPvlGnu4E+XnV6NLU14sAHw+XPRKmIdBaczz3uK17aHvGdnZm7qTQUghAcKNgYoNjiQmJITY4lqjgKPw8/Rya8VxzK7n0noOo4MhG2LMMLn1KCsMF0FpzurCEEzlnGvmqG/iyl5zCkipfSylo0sjT1phHBPvRq23T0tveBPl5Gj99vQhq7EWQrxeNvMz5hijqRl5xHjsyd9jGCbZnbCc931il0tvdm25B3RjfZTwxITF0D+5OmF+YU4+9SHGoT1a+CI2CoF+VK6s2WBarJjO3kLTsQtJOF5CabVxPPV1AWnYBaacLSc0u4ERuUbUDsV4ebjQr8829bTNfmvp60czPi6Z+XuXuC/Lzoomvl3Tb1GNWbeXAqQPGWEHpnsG+U/uwamOcp21AW/qG97V1D3Vu2hlPd0+TU9eOFIf64tCvsP9HGP538PY3O41DWK2aE3lFpGYXGI19mYY+NbuQ9NJCkJ5TWGVXTpCfF6H+3jQP8KFLc3+aNfautrH39XKNY9+FfWTmZ5YbNN6RsYOc4hwA/L38iQ2OZVibYcagcXAMTXyamBu4Drh0cZAB6VJaw8oXoHFz6H2n2WkumtWqOZlXRGrZb/el11OzC0k7XUhadgHppwurPKa+qa8nzQN8CA3woVNzf5oHGAUg1N+H0NLrIY298fKQ/npRWZGliN0ndtsKwbb0bRzJMVYvdlfudG7amavaX0VsSCwxwTG0DWiLm6p/f0suXRzkJLhSB1bBoV/gitfBy9fsNNXSWnMyr7j0m32BrZFPLfOt/0w3T1WNfhNfT5qXNvAdQ4JtjX7zAG9C/M/89MbbQ/ruzZaRn8G6Y+tYd3QdG1M3UmApMDtSjWUVZlFsLQYgzC+MmOAYbuxyI7EhsXRr1o1GHo1MTugYLl0cBGf3GgJaQa9JpsUoLLFwPKuAY1kFHMvK51hWAcezCsr07Rv9/VX16Qc28rQ19O1DmhkNfml3T2iAN6H+PoT4e+PjKY2+s8orzmND6gajIBxbR/LJZAACvALoE9bHpbpZArwCbEcRhfqGmh3HNFIcXF3yCjiyAa6eCR7edvkVhSUWUrMKbY1+2QJwLCufY6cKyMwtqvS8AB+P0m/2PvSN8CO09Fu+0cVT2r0jjb5LKrGWsCNzB+uOrmPtsbVsTd9KibUELzcvejbvyYNxD9K/RX+6Nu1q2nH64uJIcXBlVqux19A0AnrcdEEvUVRiJTW7QoN/qmwRKCAjp7DS8wJ8PAgPbER4Ex9iWjYhPNCn9GJsCwvwwc9b/rzqC601B7MP2rqKfj/+O6eLTwPQLagbt0TeQr/wfsSFxuHj4WNyWlEX5H+vK9u9FI5vg2vfhyoOkyu2nG34j57Kr9Ttc6bhr3gepL+Ph62hj2oRYGvwbY1/oDT8DUFGfgbrj623dRWdmeithV8LLm93Of1a9KNvWF+a+sg5NfWR/A93UcXFxfDDC5QEduAHPZBjifs5eqqgtAAYjX96VQ2/twdhgT6EN2lEt7AAwpv40CKwEWGBPrRo4kNYYCMaS8PfIOUV57EpbZOtq2jvyb2A0QffN7wvd8XcRf/w/rTybyWH9TYALt0K1PdDWbXWpJ0uZOfRbHYey2bXsWwOn8zneFY+A3J/5A3PvTxY9ADLF2wDwM/LnfAmxjf7LmH+tm/54U0a0SLQh7BAH/x9XOtEHGE/FqvFGDco3TPYkraFYmsxnm6exIXG8WDcg/QL70e3oG4ybtAAydxKTqLEYuVARi47jxpFYOexbHYezS430Ns6qBHtmvnRwt+DaQcmgacvW69aQngTP8Kb+BAgDb84B601f57+k7VH17Lu2Dp+O/4bp4uMcYOuQV3pF96P/uH96dm8Z4M5XLOhk7mVnMzpgmJ2Hz9tFIHSvYLdx0/bplj2cnejc1hjhnULJTI8gMgWgXQN9z/b+G+aCzsPw7ULSOgi6zWI6p0oOGEbN1h7dC3Hco8BEO4XzvC2w+kX3o++4X0J8gkyOalwNlIc7EhrzbGsgnJFYOexbA5l5tke09TXk8gWAUzq35bIFgFEhgfSPsQPz+pm2ywpNFZ5a9kLOo900DsRriK/JJ9NqZtsXUW7T+wGjCke+ob15Y7oO+jXoh9t/NvIuIE4JykOdaTYYmVfWk6lQnAqr9j2mHbNfIlqEcD1vVoR2SKAbuEBhAX41O4/6aa5kHUYRs80pv4UDZrFamHXiV22rqLNaZspthbj4eZBz9Ce3N/zfvqH9yeyWaSMG4hakeJwAbLyi9ldZlxg57FsklNzKLIY3ULeHm50DfPniugwIsONItA1PODijwIqzofEGdB2ILS/tA7eiSi2FnM85zgpOSm2KRNcwfHc46w9upb1x9fbxg26NO3CTV1vol8L43wDX0/nnUpFOD8pDuegtSblZH65AeKdx7JJOZlve0wzPy8iWwRw2yXtjPGB8AAigv3sswjL7x9AznEYN1v2GmohrziPw6cPk3I6hcOnD5e7HMs9hkVXv9KaM2vu25xhbYbRP7w/fcL7OGz5SNEwSHEoVVRiJTnt9NkuodKjhrILjMVclIKIYD+6t27ChD5tiGwRQFR4ACH+3o7pu83cD6v/AR2GQruB9v99LkRrzYmCE7YGv2IRyCzILPf4QO9AWjduTUxwDFdEXEFr/9a08m/lUkfoBHgF0Nq/tYwbCLtp0MXh573pfL3lCDuPZrMvLcc2E2gjT3e6hPkzqnuL0qOFAuga5o+vl0kf16k/Ye4Yo0KNfNWcDCYrsZZwPPd4tQUgr+TsIL9C0dyvOa39WzOk9RBb49/avzWt/VsT4BVg4jsRwjU4XXFQSrUHngQCtdbj7Pm7klNPszo5g8jwAC7tGmorBO2a+TnPKl7Zx+Cj0VCYDZOWQkhnsxPZTX5JfqVG/8ztozlHKdFnl+T0dPO0NfjxYfG2hr+VfytaNm6Jt7t9JiEUoqFwyElwSqnZwCggTWsdXWb7SOBfgDvwX631K2XuW1jT4nChJ8FZrNp5ikBVcjPgwysh+wjc8hW07m12oouitSarMKtSv/+ZIpCWn1bu8f5e/rZGv+Il1De0Xi6wUheKi4tJSUmhoMB11lAQ9uXj40OrVq3w9Cx/oqwznAQ3B3gbmFsmlDvwDjAcSAF+V0ot0VrvdFAm5y4M+Sdh7jVGl9LNC12qMBRZithzYg97T+6tVADOzOR5RmijUFr5t2JAywGVCkCgd6BJ78C1paSk4O/vT7t27WRMQqC1JjMzk5SUFCIiImr8PIcUB611olKqXYXNfYB9WusDAEqpBcAYoEbFQSk1BZgC0KZNm7oL6wwKsuGT6yBjD0yYD+0uMTtRtSxWC/uz9rMjYwdJGUkkZSax9+ReSqxGF5CHmwctG7eklX8ruod0L9f4t/JvJdM720FBQYEUBmGjlKJZs2akp6fX6nlmjjm0BA6XuZ0C9FVKNQNeBHoqpaZrrV+u6sla61nALDC6lewd1mGK8uDT8XB0C4z/GDpeZnYiG601KadTSMpMMgpBRhK7Tuwiv8Q4tLexZ2OigqOYFDmJ6OBougZ1JdwvXE6+MoEUBlHWhfw9ON2AtNY6E7inJo+td7OyFhfAgpvg8Dq47r/Q9SpT46Tnpdv2BnZk7CApM4mswiwAvN296RrUlbGdxhLVLIro4Oh6u9C6EA2RmcXhCNC6zO1WpdtqTGu9FFgaHx9/V10GM4WlGL6YDAd+gjHvQvR1Dv312UXZ7MjYwY5Mo3toe8Z20vKMAWJ35U7HJh25rM1lRAVHEd0smo5NO+LpJrPAiuqlpqby8MMPs27dOpo2bYqXlxePPfYY1157rcMyrFu3jgcffJDCwkIKCwsZP348zz77bLWP37BhA3PnzmXmzJmV7mvXrh0bNmwgOLj6kw1r8hhXYWZx+B3opJSKwCgKNwK1Wuuy3uw5WC2w+C7Y+w1cOQN6TrTrr8svyWfPiT1sz9hOUkYSOzJ3cCj7kO3+tgFtiW8eT3RwtK17yJVOEBPm01pzzTXXMGnSJD799FMADh06xJIlSxyaY9KkSXz++ed0794di8XCnj17zvn4+Ph44uOrPHinwXFIcVBKzQcSgGClVArwjNb6A6XUVOA7jENZZ2utd9TmdevFnoPVCl9PhR1fwvC/Q5+6fSvF1mL2n9rP9ozttkHjfaf22aaMCPUNJbpZNGM6jCEqOIqoZlFylJC4aCtXrsTLy4t77jnbQ9y2bVvuv/9+Dh48yC233EJubi4Ab7/9NgMGDGDVqlXMmDGD//3vfwBMnTqV+Ph4Jk+ezLRp01iyZAkeHh5cfvnlzJgxgy+++ILnnnsOd3d3AgMDSUxMrJQjLS2N8PBwANzd3YmMjAQgNzeX+++/n6SkJIqLi3n22WcZM2ZMuQyZmZlMmDCBI0eO0L9/f8oe9v/JJ58wc+ZMioqK6Nu3L++++y7u7vVrbM1RRytNqGb7cmC5IzI4Ja1h+SOw9VNIeAIGPnBRL2fVVg5lH7LtDSRlJLH7xG4KLYWAMeVCdHA0g1sNtu0VhPqG1sU7EU7quaU72Hk0u05fM7JFAM9cHXXOx+zYsYO4uLgq7wsNDeX777/Hx8eH5ORkJkyYwLnOU8rMzOTLL79k9+7dKKU4deoUAM8//zzfffcdLVu2tG2r6OGHH6ZLly4kJCQwcuRIJk2ahI+PDy+++CJDhw5l9uzZnDp1ij59+nDZZeUP/njuuee45JJLePrpp1m2bBkffPABALt27eKzzz7jl19+wdPTk/vuu4958+Zx6623nvMzcTVONyBdGy7draQ1rHgKNnwAAx+EIY/V8uma1LxU21FDSZlJ7MzYaTuPoJFHI7oFdeOGLjcQ3SyamOAYWftXmOYvf/kLa9aswcvLix9++IGpU6eyZcsW3N3d2bt37zmfGxgYiI+PD3fccQejRo1i1KhRAAwcOJDJkydzww03MHbs2Cqf+/TTTzNx4kRWrFjBp59+yvz581m1ahUrVqxgyZIlzJgxAzAO//3zzz/LPTcxMZHFixcDcNVVV9G0aVMAfvzxRzZu3Ejv3sa5R/n5+YSG1r8vWS5dHFy6W2nVy7D2begzBS57rkazrG5J28K6Y+tsBeHMhHIeyoNOTTsxMmKkbY+gfWB7PNxc+p9X1IHzfcO3l6ioKBYtWmS7/c4775CRkUF8fDxvvPEGzZs3Z+vWrVitVnx8jHNdPDw8sFqttuecOcPbw8OD3377jR9//JGFCxfy9ttvs3LlSt577z3Wr1/PsmXL6NWrFxs3buSRRx5h8+bNtGjRguXLjU6JDh06cO+993LXXXcREhJCZmYmWmsWLVpEly5dyuVOTU0973vTWjNp0iRefrnKo+zrDZc+7lApdbVSalZWVpbZUWpnzRvw86vQ82ZjIr3zFIaCkgL+vvbv3PLNLby75V3+PP0nA1oMYFqfacy7ch7rJq7j86s/5+n+TzO201g6N+0shUGYaujQoRQUFPDvf//bti0vz5gcMSsri/DwcNzc3Pj444+xWIzxr7Zt27Jz504KCws5deoUP/74IwA5OTlkZWVx5ZVX8sYbb7B161YA9u/fT9++fXn++ecJCQnh8OHDfPjhh2zZssVWGJYtW2YbK0hOTsbd3Z0mTZowYsQI3nrrLdt9mzdvrvQeBg8ebBtM/+abbzh58iQAw4YNY+HChaSlGUfznThxgkOHDlV6vqtz6RbEJfcc1r8PPzwL0ePg6pngdu76vO/kPh5NfJR9p/YxOWoyU2Kn4O/l75isQlwgpRRfffUVDz/8MK+99hohISH4+fnx6quvEhcXx3XXXcfcuXMZOXIkfn5+ALRu3ZobbriB6OhoIiIi6NmzJwCnT59mzJgxFBQUoLXmn//8JwCPPvooycnJaK0ZNmwY3bt3r5Tj448/5uGHH8bX1xcPDw/mzZuHu7s7f/vb33jooYeIjY3FarUSERFhGwg/45lnnmHChAlERUUxYMAA20wMkZGRvPDCC1x++eVYrVY8PT155513aNu2rT0/UodzyMR79nahE+853Ka5sOR+6DoKrp8D7tWfJ6C15ou9X/Da76/h5+nHS5e8xMCWso6DOL9du3bRrVs3s2MIJ1PV34UzTLwntn0BSx4wpsMYN/uchSGrMIvn1j7H94e+p394f14a9JKs8iWEcCiXLg4uc7TSrqXw5d3GBHrjPwGP6tca2JS6icdXP05GXgZ/7fVXJkVNkikphBAO59KtjtZ6qdZ6SmCgE5+0lfw9fHEbtOxlzLDqWfWZxharhfe2vsdt392Gh/Lg4ys/5rbo26QwCCFM4dJ7Dk7vj0T47GZoHgkTvwDvqgeSj+ceZ/rq6WxI3cCVEVfyt35/o7FXYweHFUKIs6Q42Muf6+HTG6FpBNz8JTRqUuXDfvrzJ/72698oshTxwsAXGN1htJyoJoQwnUsXB6cdczi6GeaNA/8wuPUr8GtW6SGFlkL+seEfzN89n25B3Xht8Gu0C2zn8KhCCFEVl+7Qdsoxh9Sd8PG14NMEJi0xCkQFB04d4KZlNzF/93xuibyFT678RAqDqFdefPFFoqKiiI2NpUePHqxfvx6AN99803YynFlOnTrFu+++W+397u7u9OjRg6ioKLp3784//vGPcmdum+Gll15y/C/VWrv8pVevXtoppCdr/VpHrWd00Tpzf6W7rVarXrR3ke79SW89aP4g/fPhn00IKeq7nTt3mvr7f/31V92vXz9dUFCgtdY6PT1dHzlyRGutddu2bXV6erqZ8fQff/yho6Kiqr3fz8/Pdj01NVUPGzZMP/30046IVq2ymS5UVX8XwAZdTbvq0nsOTuXkIZg7GrQVbl0CQe3L3Z1dlM1jiY/xzK/PEBscy8LRCxncarBJYYWwn2PHjhEcHIy3t3HIdnBwMC1atGDmzJkcPXqUSy+9lEsvvRSAFStW0L9/f+Li4rj++uvJyckBYOPGjQwZMoRevXoxYsQIjh07BkBCQgIPPvggPXr0IDo6mt9++w0wpuC+/fbb6dOnDz179uTrr78GjNlh+/TpQ48ePYiNjSU5OZlp06axf/9+evTowaOPPnrO9xIaGsqsWbN4++230VpjsVh49NFH6d27N7Gxsbz//vu29zx48GBbrtWrVwPw7bffEhcXR/fu3Rk2bNg5s86ZM4exY8cycuRIOnXqxGOPGZNxTps2jfz8fHr06MHEifZd66UsOUO6LmQfhdkjoSALJv8PwmLK3b0lbQvTVk/jeO5xpvacym1Rt8m6ysJuyp0J+800OL69bn9BWAxc8Uq1d+fk5HDJJZeQl5fHZZddxvjx4xkyZAhQfqW0jIwMxo4dyzfffGObWqOwsJDp06czZMgQvv76a0JCQvjss8/47rvvmD17NgkJCXTq1In//Oc/JCYmct9995GUlMQTTzxBZGQkN998s20K7s2bNzNt2jT69evHxIkTKSoqwmKxkJqayqhRo0hKSqoyf+PGjW1F6owmTZqwZ88evv76a9LS0njqqacoLCxk4MCBfPHFFyxevJiCggKefPJJLBYLeXl5FBQUEBcXR2JiIhEREZw4cYKgoKBqs37xxRc8//zzbN68GW9vb7p06cKaNWto3bp1lZlqq0GdIe0UA9I5afDRaMg7AZO+LlcYLFYLH+74kLc3v02YXxgfXfER3UMqz/8iRH3SuHFjNm7cyOrVq/npp58YP348r7zyCpMnTy73uHXr1rFz504GDjSmhSkqKqJ///7s2bOHpKQkhg8fDoDFYrEt2AMwYYKxPMzgwYPJzs7m1KlT1U7B3b9/f1588UVSUlIYO3YsnTp1uqj3tmLFCrZt28bChQsBYxLB5ORkevfuze23305xcTHXXHMNPXr0YNWqVQwePJiIiAgAgoKCbK9R3XThw4YN48wYamRkJIcOHaJ169YVYziESxcHbfbEe3knYO41kH0Ebl5snOhWKi0vjSdWP8H64+sZ2W4kT/d/WibME453jm/49uTu7k5CQgIJCQnExMTw0UcfVSoOWmuGDx/O/Pnzy23fvn07UVFRrF27tsrXrniot1Kq2im4u3XrRt++fVm2bBlXXnkl77//Pu3bl+/yPZ8DBw7g7u5OaGgoWmveeustRowYUelxiYmJLFu2jMmTJ/PXv/7Vtv5DRdVlXb9+va0rDozPsKSkpFZZ65KMOVyogiz4ZCxk7oMbP4W2/W13JaYkMm7JOLZlbOP5Ac/z2uDXpDCIBmPPnj0kJyfbbm/ZssU2Y6m/vz+nTxsLUvXr149ffvmFffv2AUZf/N69e+nSpQvp6em24lBcXMyOHWdXEP7ss88AWLNmDYGBgQQGBlY7BfeBAwdo3749DzzwAGPGjGHbtm3lMpxPeno699xzD1OnTkUpxYgRI/j3v/9NcXExAHv37iU3N5dDhw7RvHlz7rrrLu688042bdpEv379SExM5I8//gCMqb2BGk0XXpGnp6ftdzqKS+85mKYoF+bdYPTljp8HHYzBtSJLEW9sfINPdn1C56adeX3w67RvUrtvKUK4upycHO6//35OnTqFh4cHHTt2ZNasWQBMmTKFkSNH0qJFC3766SfmzJnDhAkTKCw0lrJ94YUX6Ny5MwsXLuSBBx4gKyuLkpISHnroIaKijIWLfHx86NmzJ8XFxcyePRug2im4P//8cz7++GM8PT0JCwvjiSeeICgoiIEDBxIdHc0VV1zB66+/Xi7/mcHf4uJiPDw8uOWWW/jrX/8KwJ133snBgweJi4tDa01ISAhfffUVq1at4vXXX8fT05PGjRszd+5cQkJCmDVrFmPHjsVqtdqWR63JdOEVTZkyhdjYWOLi4pg3b16d/ntVRwaka6u4AD69AQ6uNmZXjboWgINZB3ks8TF2ndjFTV1v4q/xf8XbvfoJ9oSwl/o8ZXdCQgIzZswgPr7KMVRxDg1qQNrhSorg81uNOZOufQ+irkVrzdf7v+al9S/h7e7NW0PfIqF1gtlJhRDiokhxqClLCSy6A5K/g1FvQPcbySnK4e/r/s7yP5bTO6w3L1/yMs39mpudVIh6a9WqVWZHaDBcujg47FBWqxW+vg92LYERL0P87WxP385jiY9xLPcYU3tM5c6YO+XcBSFEveHSRytpR8ytpDUsexi2fQZDn8La7x5mJ83m1m9uxaItfDjyQ+7ufrcUBiFEveLSew52pzV8Ox02zoFB/0dG79t44vt7WHtsLcPbDueZ/s8Q6O1Ek/4JIUQdkeJwLiv/Duv/DX3v5ZcuQ3liyXXkFufydP+nGddpnKy7IISot1y6W8muEl+H1f+gOO5W/hEczD0/3kuQTxALrlrA9Z2vl8IgxDkcP36cG2+8kQ4dOtCrVy+uvPJK9u7de0Gvdeedd7Jz5846zXe+aburk5OTw7333kuHDh2Ii4ujV69e/Oc//7moLHPmzGHq1KkAvPfee8ydO/eCXufgwYN8+umnF5WlLCkOVVn7Lqx8gT+jR3OLWzpzdn7E+C7jmX/VfDo2dbKFhYRwMlprrr32WhISEti/fz8bN27k5ZdfJjU19YJe77///S+RkZF1mvFCi8Odd95J06ZNSU5OZtOmTXz77be2M5/LutBpL+655x5uvfXWC3quFAd72zAbvpvO0s6XcH3hXg6fPsybCW/yVL+n8PHwMTudEE7vp59+wtPTk3vuuce2rXv37gwaNAitNY8++ijR0dHExMTYpsJYtWoVCQkJjBs3jq5duzJx4kTb9BIJCQmcOcm1ceOza6svXLjQNl/T/v376devHzExMTz11FO2x+Xk5DBs2DDi4uKIiYmxTY9d1bTdr7/+um0q7meeeabS+9q/fz+//fYbL7zwAm5uRtMZEhLC448/bnsPgwYNYvTo0bZids0119CrVy+ioqJsZ4kDfPjhh3Tu3Jk+ffrwyy+/2LY/++yztgn59u/fz8iRI+nVqxeDBg1i9+7dAEyePJkHHniAAQMG0L59e9skgNOmTWP16tX06NGDN954o5b/apXJmENZWxeQu+z/eKl9DEuK/yQuNI5XBr1CeOPw8z9XCCf06m+vsvvE7jp9za5BXXm8z+PV3p+UlESvXr2qvG/x4sVs2bKFrVu3kpGRQe/evRk82FjXZPPmzezYsYMWLVowcOBAfvnlFy655JIaZXrwwQd58MEHmTBhAu+9955tu4+PD19++SUBAQFkZGTQr18/Ro8ezSuvvEJSUhJbtmwBjJlSk5OT+e2339BaM3r0aBITE23ZwFgbonv37rbCUJVNmzaRlJRkm4l19uzZBAUFkZ+fT+/evbnuuusoKirimWeeYePGjQQGBnLppZfSs2fPSq81ZcoU3nvvPTp16sT69eu57777WLlyJWCsH7FmzRp2797N6NGjGTduHK+88gozZsw471QcNSXF4YwdX7Jj+QM81q4dKZzmvu73cVfsXXi4yUckRF1Zs2YNEyZMwN3dnebNmzNkyBB+//13AgIC6NOnD61atQKgR48eHDx4sMbFYe3atXz11VcA3HTTTTzyyCOA0cX1xBNPkJiYiJubG0eOHKmye2vFihWsWLHC1kjn5OSQnJxcrjhU9OKLL/LFF1+QlpbG0aNHAejTp4+tMADMnDmTL7/8EoDDhw+TnJzM8ePHSUhIICQkBIDx48dXGo/Jycnh119/5frrr7dtOzP/FBh7JG5ubkRGRl5wd935OF3Lp5TyA94FioBVWmu7zzJl3b2cj1c8yJvhzWnWKIgPBr9KfJjM3SJc37m+4dtLVFSUraujNmoyXXXZA0EKCgrO+5rz5s0jPT2djRs34unpSbt27ap8ntaa6dOnc/fdd1f7WpGRkWzduhWr1YqbmxtPPvkkTz75ZLmuLj8/P9v1VatW8cMPP7B27Vp8fX1JSEioUWYAq9VKkyZNbHs2FZX9rOw1P55DxhyUUrOVUmlKqaQK20cqpfYopfYppaaVbh4LLNRa3wWMtne2zF1L+Muqh5gRFMjgVoNYNGaxFAYhLsLQoUMpLCws18e+bds2Vq9ezaBBg/jss8+wWCykp6eTmJhInz59avzazZs3Z9euXVitVts3cjCm/160aBEACxYssG3PysoiNDQUT09PfvrpJw4dOgRQadruESNGMHv2bNtqa0eOHCEtLa3c7+7YsSPx8fE89dRTWCwWwChQ1TXOWVlZNG3aFF9fX3bv3s26desA6Nu3Lz///DOZmZkUFxfzxRdfVHpuQEAAERERtvu01mzduvWcn01tpiKvCUcNSM8BRpbdoJRyB94BrgAigQlKqUigFXC49GEWe4Zat2U249ZO5zcfb57s+RBvDntHTmoT4iIppfjyyy/54Ycf6NChA1FRUUyfPp2wsDCuvfZaYmNj6d69O0OHDuW1114jLCysRq8J8MorrzBq1CgGDBhQbnW4N998k3/+85/Exsayb98+22pqEydOZMOGDcTExDB37ly6du0KQLNmzWzTdj/66KNcfvnl3HTTTfTv35+YmBjGjRtXZUP73//+l8zMTFuhGD58OK+99lqVmUeOHElJSQndunWzLVcKEB4ezrPPPkv//v0ZOHBgtTPozps3jw8++IDu3bsTFRVlG0yvTmxsLO7u7nTv3r1OBqQdNmW3Uqod8D+tdXTp7f7As1rrEaW3p5c+NAU4qbX+n1Jqgdb6xmpebwowBaBNmza9znwjqI1FKx7m45QfeW3ELDq37Ffr5wvhjOrblN0xMTEsWbKkXF9+RXl5eTRq1AilFAsWLGD+/PnnbUwbGleasrslZ/cQwCgKfYGZwNtKqauApdU9WWs9C5gFxnoOFxJg7PB/MqowG28f2VsQwhkNHz6cmJiYcxYGgI0bNzJ16lS01jRp0sS2CJC4cE43IK21zgVuq8ljL3ZWVqWUFAYhnNj3339fo8cNGjTovH3yonbMPAnuCNC6zO1WpdtqzCGzsgrhgurDCo+i7lzI34OZxeF3oJNSKkIp5QXcCCypzQsopa5WSs3KysqyS0AhXJGPjw+ZmZlSIARgFIbMzEx8fGo3w4NDBqSVUvOBBCAYSAWe0Vp/oJS6EngTcAdma61fvJDXd+ga0kI4ueLiYlJSUmp8TL2o/3x8fGjVqhWenp7ltp9rQNphRyvZQ5kxh7uSk5PNjiOEEC7lXMXBpSfekzEHIYSwD5cuDkIIIezDpYuDDEgLIYR9uPSYwxlKqXTgEBAI1LZSBAMZZZ5b9jVqcp1abqvN/bXJe65857vt6KwVX08+W3M+2+p+r3y2Deezbau1DqnyVbXW9eYCzLqA52wo+9yyr1GT67XdZq+858p3vtuOziqfrXN8ttX9XvlsG95nW9XFpbuVqlDtdBu1eO7SKrad63ptt9Xm/po893y5a3K7Jlnks636dV35s63u98pnW/65DeGzraRedCtdDKXUBl3NoVzOyJXyulJWcK28rpQVXCuvK2UF++Wtb3sOF2LW+R/iVFwprytlBdfK60pZwbXyulJWsFPeBr/nIIQQojLZcxBCCFGJFAchhBCVSHEQQghRiRSHc1BKtVdKfaCUWmh2lppQSl2jlPqPUuozpdTlZuc5F6VUN6XUe0qphUqpe83Ocz5KKT+l1Aal1Cizs5yPUipBKbW69PNNMDvPuSil3JRSLyql3lJKTTI7z/kopQaVfq7/VUr9anaec1FKtVFKfaWUmq2Umlbb59fb4lD6gaQppZIqbB+plNqjlNp3vg9Ma31Aa32HfZPactVF3q+01ncB9wDjnTzrLq31PcANwEBnzlrqceBz+6Qsl6su8mogB/DBWH7XmbOOwVjoq9ieWUtz1cXf7erSv9v/AR85c1YgBliotb4d6FnrELU9G89VLsBgIA5IKrPNHdgPtAe8gK1AZOmH+L8Kl9Ayz1voYnn/AcQ5e1ZgNPANcJMzZwWGYyxGNRkY5ex/B4Bb6fOaA/OcPOs04G5H/D+r4/9jnwP+zpwVaAb8BKwEbqt1Bnv+Y5h9AdpV+HD7A9+VuT0dmF6D17F7caiLvIACXgUuc/asFV5rmTNnBV7EWJRqBfD1mcbXWfOWeZyXAxrci/1sbwZuKL3+mT2z1tVnC7QB/uPsWYFHgMGl12v9d+BBw9ISOFzmdgrQt7oHK6WaYTQMPZVS07XWL9s5X0W1ygvcD1wGBCqlOmqt37NnuApq+9kmAGMBb2C5PYNVoVZZtdZPAiilJgMZWmurXdNVVtvPdiwwAmgCvG3XZJXV9m92MfCWUmoQkGjPYNWobV6AO4AP7ZaoerXN+i3wrFLqJuBgbX9ZQysOtaK1zsTov3cJWuuZwEyzc9SE1noVsMrkGLWitZ5jdoaa0Fovxmh0nZ7WOg+jsXUZWutnzM5QE1rrJGDchT6/3g5IV+MI0LrM7Val25yVK+WVrPbjSnldKSu4Vl6HZm1oxeF3oJNSKkIp5YUxyLjE5Ezn4kp5Jav9uFJeV8oKrpXXsVntPahi1gWYDxzj7CFyd5RuvxLYizHq/6TZOV0xr2SVvK6W1dXyOkNWmXhPCCFEJQ2tW0kIIUQNSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIeykdGGWLaWX9Uop+f8mXIacIS2EnSilkjHm0z9mdhYhaku+yQhhP8uBbUqpN80OIkRtyXoOQtiBUmoAxsp84VrrErPzCFFbsucghH1cD+zVWpcoQ4DZgYSoDRlzEMIOlFJ9gA8ADeQD92mtN5qbSoiak+IghBCiEulWEkIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJX8PwHirA85xCEuAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"A = np.array([[ .2, .1, 1.0, 1.0, 0.0],\n",
|
|
" [ .1, 4.0, - 1.0, 1.0, - 1.0],\n",
|
|
" [ 1.0, -1.0, 60.0, .0, .0],\n",
|
|
" [ 1.0, 1.0, .0, 8.0, 4.0],\n",
|
|
" [ .0, -1.0, - 2.0, 4.0, 700.0]])\n",
|
|
"b = np.array( [ 1 , 2 , 3 , 4 , 5 ] )\n",
|
|
"x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
"eps_list = np.logspace(-8, -1, 8)\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"\n",
|
|
"for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\")]:\n",
|
|
" #diff_list = []\n",
|
|
" k_list = []\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" #diff_list.append(diff(x_exact, x))\n",
|
|
" k_list.append(k)\n",
|
|
" #ax.plot(eps_list, diff_list, label=alg_name)\n",
|
|
" ax.plot(eps_list, k_list, label=alg_name)\n",
|
|
" \n",
|
|
"ax.set_xscale(\"log\")\n",
|
|
"ax.invert_xaxis()\n",
|
|
"ax.set_yscale(\"log\")\n",
|
|
"ax.set_xlabel(\"$\\epsilon$\")\n",
|
|
"#ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n",
|
|
"ax.set_ylabel(\"$k$\")\n",
|
|
"ax.legend()\n",
|
|
"\n",
|
|
"fig.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"The condition number for A is K(A) = 12269.877667702964 >> 1,\n",
|
|
"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\n",
|
|
"This explains the great difference in order of required iterations k as observed above.\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(\"The condition number for A is K(A) =\", linalg.cond(A), \">> 1,\")\n",
|
|
"print(\"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\")\n",
|
|
"print(\"This explains the great difference in order of required iterations k as observed above.\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "4d7421993a0fef53e4f0f83ffae975b1",
|
|
"grade": false,
|
|
"grade_id": "cell-c1b8b533da043a26",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"### Task 8\n",
|
|
"\n",
|
|
"Try to get a better convergence behavior by pre-conditioning your matrix $A$. Instead of $A$ use\n",
|
|
"\n",
|
|
"$$ \\tilde{A} = C A C,$$\n",
|
|
"\n",
|
|
"where $C = \\sqrt{D^{-1}}$. If you do so, you will need to replace $\\mathbf{b}$ by \n",
|
|
"\n",
|
|
"$$\\mathbf{\\tilde{b}} = C \\mathbf{b}$$\n",
|
|
"\n",
|
|
"and the vector $\\mathbf{\\tilde{x}}$ returned by your function will have to be transformed back via\n",
|
|
"\n",
|
|
"$$\\mathbf{x} = C \\mathbf{\\tilde{x}}.$$ \n",
|
|
"\n",
|
|
"What is the effect of $C$ on the condition number and why?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "213ac0f9b57d67bb90bbd75dd3808955",
|
|
"grade": true,
|
|
"grade_id": "cell-de3d1bd5ccfa8dfb",
|
|
"locked": false,
|
|
"points": 3,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"The number of iterations is brought down by a lot due to the conditioning, bringing down the runtime.\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAF3CAYAAAB65myjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACZ/ElEQVR4nOzdd3xUdb7/8dd3JpPeG0kICYHQQgkkIYAUsYBYsDdQV10VXXV19a6/te3q7t2967re7XvXRUXXXTs2bAgWxIJCQpNeAwQCCek9U76/P85kMgmhJzkzyef5eOQxM2fOnHwGY755z7cprTVCCCGEEEIIIYQ/sphdgBBCCCGEEEIIcaok1AohhBBCCCGE8FsSaoUQQgghhBBC+C0JtUIIIYQQQggh/JaEWiGEEEIIIYQQfktCrRBCCCGEEEIIvxVgdgFdIT4+Xg8cONDsMoQQQvQShYWFh7XWCWbX4Y+UUrOB2REREbcNHTrU7HKEEEL0Esdqm1Vv2Kc2Ly9PFxQUmF2GEEKIXkIpVai1zjO7Dn8mbbMQQoiudKy22eeGHyulpiulvlRKPa2Umm52PUIIIYQQQgghfFePhFql1AKlVKlSakOH47OUUluVUjuUUg+6D2ugDggGinuiPiGEEEIIIYQQ/qmnempfAGZ5H1BKWYG/A+cDWcAcpVQW8KXW+nzgZ8Ave6g+IYQQQgghhBB+qEcWitJaL1dKDexwOB/YobXeBaCUehW4RGu9yf18JRB0qt/TbrdTXFxMU1PTqV5C9DLBwcGkpqZis9nMLkUIIfokaZtFR9I2CyG6gpmrH/cH9nk9LgYmKKUuB84DooG/He3FSql5wDyAtLS0I54vLi4mIiKCgQMHopTqwrKFP9JaU15eTnFxMRkZGWaXI4QQfZK0zcKbtM1CiK7icwtFaa3f0lrfrrW+Rmu97Bjnzdda52mt8xISjlzZuampibi4OGk0BQBKKeLi4qR3QAghTCRts/AmbbMQoquYGWr3AwO8Hqe6j3UZaTSFN/l5EEII88nvYuFNfh6EEF3BzFC7ChiilMpQSgUC1wKLTuYCSqnZSqn51dXV3VJgVzh06BBz585l0KBB5ObmMmnSJN5+++0ereHbb79lwoQJjB07lhEjRvD4448f8/yCggLuueeeTp8bOHAghw8fPubrT+QcIYQQwizSNgshRO/SI3NqlVKvANOBeKVUMfCY1vo5pdTdwMeAFVigtd54MtfVWr8HvJeXl3dbV9fcFbTWXHrppdx44428/PLLAOzZs4dFi04qu5+2G2+8kddff53s7GycTidbt2495vl5eXnk5XW6r7EQQgjh16RtFkKI3qdHemq11nO01slaa5vWOlVr/Zz7+Ida66Fa68Fa69/0RC096bPPPiMwMJA77rjDcyw9PZ0f//jHFBUVMXXqVHJycsjJyeGbb74BYNmyZVx00UWe8++++25eeOEFAB588EGysrIYM2YMP/3pTwF44403GDVqFNnZ2UybNq3TOkpLS0lOTgbAarWSlZUFQH19PT/84Q/Jz89n3LhxvPvuu0fUUF5ezsyZMxk5ciS33norWmvPdf/zn/+Qn5/P2LFjuf3223E6nV3xzyaEEEJ0G2mbhRCi9zFz9ePTppSaDczOzMw85nm/fG8jmw7UdOn3zkqJ5LHZI495zsaNG8nJyen0ucTERJYuXUpwcDDbt29nzpw5FBQUHPVa5eXlvP3222zZsgWlFFVVVQD86le/4uOPP6Z///6eYx3dd999DBs2jOnTpzNr1ixuvPFGgoOD+c1vfsPZZ5/NggULqKqqIj8/n3PPPbfda3/5y18yZcoUfvGLX/DBBx/w3HPPAbB582Zee+01vv76a2w2G3feeScvvfQSP/jBD475byKEEEKAtM3SNgshRNfxudWPT4bW+j2t9byoqCizSzkhd911F9nZ2YwfPx673c5tt93G6NGjueqqq9i0adMxXxsVFUVwcDC33HILb731FqGhoQBMnjyZm266iWeeeeaon8b+4he/oKCggJkzZ/Lyyy8za9YsAJYsWcITTzzB2LFjmT59Ok1NTezdu7fda5cvX871118PwIUXXkhMTAwAn376KYWFhYwfP56xY8fy6aefsmvXrtP69xFCiFPmckLtQXA6zK5E+Blpm4UQopu4XFBzwGiju5lf99SeqON9attdRo4cyZtvvul5/Pe//53Dhw+Tl5fHH//4R/r168e6detwuVwEBwcDEBAQgMvl8rymdZn7gIAAVq5cyaeffsrChQv529/+xmeffcbTTz/Nd999xwcffEBubi6FhYX89Kc/Zc2aNaSkpPDhhx8CMHjwYH70ox9x2223kZCQQHl5OVpr3nzzTYYNG9au7kOHDh33vWmtufHGG/ntb3972v9OQghxTE4H1B00Gsaa/e7bDvdrS8DlgLsLIH6I2RWLEyBts7TNQgg/pzXUlULVXqja4/7aC5Xu2+p94GyBe9ZCbPfuRe3XofZEhx+b5eyzz+bhhx/mH//4Bz/60Y8AaGhoAKC6uprU1FQsFgv/+te/PJ/kpqens2nTJpqbm2lsbOTTTz9lypQp1NXV0dDQwAUXXMDkyZMZNGgQADt37mTChAlMmDCBjz76iH379vH888+3q+ODDz7gggsuQCnF9u3bsVqtREdHc9555/HXv/6Vv/71ryilWLNmDePGjWv32mnTpvHyyy/z6KOP8tFHH1FZWQnAOeecwyWXXMJ9991HYmIiFRUV1NbWkp6e3q3/pkKIXsbRYgTSTgOr+3HdIdCu9q+zhUJkivE1cIr7fn8IiTXnffRySqkw4P+AFmCZ1volk0s6ZdI2CyHECdIaGiqgqqh9WPWE2L3g6LDPdGg8RKdB8hgYcZFxP7j7R9X6daj19dWPlVK888473HfffTz55JMkJCQQFhbG7373O3Jycrjiiit48cUXmTVrFmFhYQAMGDCAq6++mlGjRpGRkeFpyGpra7nkkktoampCa80f/vAHAB544AG2b9+O1ppzzjmH7OzsI+r497//zX333UdoaCgBAQG89NJLWK1Wfv7zn/OTn/yEMWPG4HK5yMjI4P3332/32scee4w5c+YwcuRIzjjjDNLS0gDIysri17/+NTNnzsTlcmGz2fj73/8uDacQoo290SukHqWXtb70yNcFRbYF1sQRRlhtDa2tx4OjQfa3PC1KqQXARUCp1nqU1/FZwJ8xdiZ4Vmv9BHA5sFBr/Z5S6jXAb0OttM1CCOGlsbJDYPUKrpV7wF7f/vzgaIhJh4RhMGQmRKcbwTUmHaIGQFC4KW9Dea+Y56/y8vJ0x4UcNm/ezIgRI0yqSPgq+bkQoou01HcIql6Btdrdy9pYceTrgqPbh9PO7gdH9vjb6UgpVai17tX7pyilpgF1wIutoVYpZQW2ATOAYow95ecAlwAfaa3XKqVe1lrPPd71pW0WJ0p+LoToRk01R/aueve4Nle3Pz8wwgio3mE1Oq3tqwd6XY/mWG2zX/fUCiGE6AZN1ceev1qz3zino9A4I5RGpcKA/A6htT9EJkNgWM+/H9EprfVypdTADofzgR1a610ASqlXMQJtMZAKrMXPF5kUQohepaXeK7Tuhcqi9gG2sbL9+bbQtsCaPskdVr2Ca0iMX46E8utQ6+tzaoUQwuc4WoyFGyp2QXVx56G1pbbDixSEJxrhNHZQ+zmsrb2sESlgCzblLYku1R/Y5/W4GJgA/AX4m1LqQuC9o71YKTUPmAd4hsQKIYQ4DfYmo92u3NN+MabWHteGw+3Ptwa19bD2z+3Q45pufADth6H1ePw61Pr6nFohhDCFvcn4pLZi15Ff1fvaL7qkLBCeZATThGEw+Gx3b2v/ttAangQBgaa9HWE+rXU9cPMJnDcfmA/G8OPurksIIfye1saWdBU7oXynEVq9hwfXHWx/vsUG0QOMoDr8QndgHdjW4xqWAJa+N6DGr0OtEEL0Wc11ULnbK7Dubrut2Q945YngKIgdDKnjYcw1Rm9rbIaxoEN4P7BKUyA89gMDvB6nuo8JIYQ4Va2rCJfvaAuvnvu72i/GpKzGNJ7oNMg812tOq/s2IgksVvPei4+Sv2SEEMJXNVYZQbVyd4fgusvY5sZbWIIRVjOmQkyGO7i6w2uobHMjTtgqYIhSKgMjzF4LHHdRKG8yNUgI0Wc1VbsD684jw6v3WhTKaoTV2MGQPtm4jXN/RabKh82nQP7FhBDCLK2f3HY2TLhi15GrB0ekGCF1yAyv0DrICLE+sGKw8C9KqVeA6UC8UqoYeExr/ZxS6m7gY4wtfRZorTeezHVlapAQoldrqTfaaE9g9brfbn6rMnpc4wbDqCvdoTXTCLAx6WC1mfYWeiO/DrW+/mnwb37zG15++WWsVisWi4V//vOfTJgwgT/96U/MmzeP0NBQ02qrqqri5Zdf5s477+z0eavVyujRo7Hb7QQEBPCDH/yA++67D4uJY/T/53/+h4cffti07y/EKdHa6FXtNLgWdVhKXxlDgmMzIOuSDsF1IASa9ztD9D5a6zlHOf4h8GEPl9NjpG3uWtI2i17J0eweHbXzyPBae6D9ueFJRmAdfoFXj2um0W7bQkwpvy+SfWq7yYoVK7j//vtZtmwZQUFBHD58mJaWFlJSUhg4cCAFBQXEx8ebVl9RUREXXXQRGzZs6PT58PBw6urqACgtLWXu3LlMnjyZX/7ylz1Z5lFrOlVm/1yIXsrlMuaxegfWyt1tw4XtDW3neoYcdehpjR1kHA8IMu99CI++sE9td/H6wPm27du3t3vO7N/B0jZ3PWmbhd9yOoxFmSp2GaHVM2R4h7E7gPeiiqFx7YcIt96PHQRBEea9hz5G9qk1QUlJCfHx8QQFGX+gtjaSf/nLXzhw4ABnnXUW8fHxfP755yxZsoTHHnuM5uZmBg8ezPPPP094eDiFhYXcf//91NXVER8fzwsvvEBycjLTp08nOzubL774AofDwYIFC8jPz6e+vp4f//jHbNiwAbvdzuOPP84ll1zCxo0bufnmm2lpacHlcvHmm2/y85//nJ07dzJ27FhmzJjB73//+6O+l8TERObPn8/48eN5/PHHcblcPPjggyxbtozm5mbuuusubr/9dkpKSrjmmmuoqanB4XDwj3/8g6lTp7J48WIefvhhnE4n8fHxfPrpp0et9YUXXmDRokU0NDSwc+dOLrvsMp588kkefPBBGhsbGTt2LCNHjuSll17qkf+OQng4HW1b4XSc31pZBM7mtnOtge6gmgEZ09rmtsYOMnpiZciR6MV8efixtM3SNos+xuWCmuIOc1zd9yuLwOVoOzco0giqqfmQPadtqHDcIGPvVuHT+kZP7UcPwsHvu/abJo2G85846tN1dXVMmTKFhoYGzj33XK655hrOPPNMgHafBh8+fJjLL7+cjz76iLCwMH73u9/R3NzMQw89xJlnnsm7775LQkICr732Gh9//DELFixg+vTpDBkyhGeeeYbly5dz5513smHDBh5++GGysrK4/vrrqaqqIj8/nzVr1vDggw8yceJErrvuOlpaWnA6nRw6dOiEPw1uFR0dzdatW3n33XcpLS3l0Ucfpbm5mcmTJ/PGG2/w1ltv0dTUxCOPPILT6aShoYGmpiZycnJYvnw5GRkZVFRUEBsbe9Ra33jjDX71q1+xZs0agoKCGDZsGF999RUDBgyQT4NFz2mqNn5nlKxr+yrf0b7xCwhpH1a9vyJTZGVCPyc9tadP2mZpm0+UtM3itLRO82nX2+r+qtwNjqa2c22hbUE11j1MuLXnNSy+V+7f2ptIT60JWj/N/fLLL/n888+55ppreOKJJ7jpppvanfftt9+yadMmJk+eDEBLSwuTJk1i69atbNiwgRkzZgDgdDpJTk72vG7OHGMq1LRp06ipqaGqqoolS5awaNEinnrqKQCamprYu3cvkyZN4je/+Q3FxcVcfvnlDBky5LTe25IlS1i/fj0LFy4EoLq6mu3btzN+/Hh++MMfYrfbufTSSxk7dizLli1j2rRpZGRkABAbG+u5Rme1ApxzzjlERUUBkJWVxZ49exgwYEDHMoToGvXlcHBd+wBbsavt+YgUSM6GYRcYDV/rUOGIJGn8hPAz0jZL2yz8XEs9lKyHA2uMr7LNxsipFq8PVlpHS8VlwpBz24fXiGRpu3upvhFqj/GpbXeyWq1Mnz6d6dOnM3r0aP71r38d0XBqrZkxYwavvPJKu+Pff/89I0eOZMWKFZ1eW3X4H1IphdaaN998k2HDhrV7bsSIEUyYMIEPPviACy64gH/+858MGjTopN7Lrl27sFqtJCYmorXmr3/9K+edd94R5y1fvpwPPviAm266ifvvv5+YmM6Haxyt1u+++84zLAyMf0OHw9Hx5UKcmtqD7cNryTpjSHGr6HQjwI69DpLHQvIYCE80rVwh/NEJL+IobbO0zUIci70RDm5oC7AH1sDhrW1zXSNSoN9ISJ/SNr81LtNYcVhGS/U55i2X1wWUUrOVUvOrq6uPf3IP27p1K94LZKxdu5b09HQAIiIiqK2tBWDixIl8/fXX7NixA4D6+nq2bdvGsGHDKCsr8zScdrudjRvbdlV47bXXAPjqq6+IiooiKiqK8847j7/+9a+0Dilfs2YNYDR6gwYN4p577uGSSy5h/fr17Wo4nrKyMu644w7uvvtulFKcd955/OMf/8ButwOwbds26uvr2bNnD/369eO2227j1ltvZfXq1UycOJHly5eze/duACoqjC1KjlbrsdhsNs/3FOKYtIaqvbD5Pfjs1/DSVfDUUPjfYfDy1fD5/0DZVhgwAWb8N/xgEfysCH6yHq75N0z7qfHprgRaIU6a1vo9rfW81l49XyJts7TNwkc5mmH/alj1HLx7N/xjCvxPf3juXPjoAdixFKLTYNr/gzmvwX9thf/aDNcvND4gy78NMs8xFlyUQNsn+XVPrS8vRlFXV8ePf/xjqqqqCAgIIDMzk/nz5wMwb948Zs2aRUpKCp9//jkvvPACc+bMobnZWGjm17/+NUOHDmXhwoXcc889VFdX43A4+MlPfsLIkSMBCA4OZty4cdjtdhYsWADAz3/+c37yk58wZswYXC4XGRkZvP/++7z++uv8+9//xmazkZSUxMMPP0xsbCyTJ09m1KhRnH/++UcsRtG68EPrtgE33HAD999/PwC33norRUVF5OTkoLUmISGBd955h2XLlvH73/8em81GeHg4L774IgkJCcyfP5/LL78cl8tFYmIiS5cuPWqtxzJv3jzGjBlDTk6OLEYh2rhcxpyZkrXte2AbK43nlRUShsPgc4xe2ORsSBolqxUK0QdJ2yxts/ABTjuUbm7fA3toI7jcH46ExELKOBh6nnGbMs5Yq0KGDYtj6BsLRfUy06dP56mnniIvT9YwOVm9+eeiT3A6oHx7hyHE66HF3bNhDYTErLbwmjwW+mXJPnHipMlCUadP2mZxonrzz0Wf53IaI6O8A+zB79t2DAiKgpSxbeE1ZZzRIysBVnRCFooSQvgfR4uxAIR3gD24ARyNxvO2UOg3CrKvbQuxCcMhINDcuoXo4054Tq0QondxuYwViNsF2PVte7UHhhsfNuff1hZgYwdJgBVdQkKtH1q2bJnZJQjRteyNxtAj7yHEhza1DUUKijRC6/hb2gJsXKbMmxHCB/ny1KDuJG2z6FO0NnYK8ATYtUbb3TpyKiDEaKtzbmwLsHGZYPHr5XyED5NQK4ToWc21R+4BW7YVtNN4PiTWGIp0xt1tATZ6oDSEQgghhBlaF1/07oEtWWvs6Q5gDTL2iM6+ti3Axg8Fq8QM0XPkp00I0X0aKoyhR94BtnxH2/MRyUZoHTG7LcBG9pehSEIIIYQZtIaaA+0D7IE10GiskI3FZmyjM/LytgCbOAKsNnPrFn2ehFohRNepPQRb3oddnxsBtmpv23PRaUZozb7WmFOTNAYi+plWqhBCCNHn1R4yel0PrDG21DmwBupLjeeU1Vh8cfiFbQG230gICDrmJYUwg1+HWlmMQggfULXP2A928yLY+y2gITodUsfD+FvdW+iMgdBYsysVQvQAaZuF8FH15VDiNQf2wBqo2e9+UhmLLWae2xZgk0bJ7gHCb/j1JDVf3uAd4ODBg1x77bUMHjyY3NxcLrjgArZt23ZK17r11lvZtGlTl9ZXVVXF//3f/5306+rq6vjRj37E4MGDycnJITc3l2eeeea0annhhRe4++67AXj66ad58cUXT+k6RUVFvPzyy6dVizgB5Tvhyz/A/Onwp1Hw8UPGXNnpD8KPVsC96+DKBTD5Xhg0XQKtEH2ItM2nR9pm0aUqi+CrP8E/z4TfD4L/XAGf/dpYyyJ9Mpz3W7h5MTxUDHd9C5f9AybMgwHjJdAKv+LXPbW+TGvNZZddxo033sirr74KwLp16zh06BBDhw496es9++yzXV2ip+G88847T+p1t956K4MGDWL79u1YLBbKyso8m8x7czgcBASc/I/YHXfccdKvadXacM6dO/eUryE6oTWUboJNi4xe2dKNxvGUHDj3cRhxMcQNNrVEIYQ4HmmbpW3uEyr3wKZ3YOPbRm8sGO31Ob+A1HxIHgPBvvmhkxCnyq97an3Z559/js1ma9cIZGdnM3XqVLTWPPDAA4waNYrRo0fz2muvAcZ2ANOnT+fKK69k+PDhXHfddWitAWNT99ZN7MPDwz3XXLhwITfddBMAO3fuZOLEiYwePZpHH33Uc15dXR3nnHMOOTk5jB49mnfffReABx98kJ07dzJ27FgeeOABAH7/+98zfvx4xowZw2OPPXbE+9q5cycrV67k17/+NRb3arQJCQn87Gc/87yHqVOncvHFF5OVlQXApZdeSm5uLiNHjmT+/Pmeaz3//PMMHTqU/Px8vv76a8/xxx9/nKeeesrz/WbNmkVubi5Tp05ly5YtANx0003cc889nHHGGQwaNIiFCxd63tOXX37J2LFj+eMf/3iS/9VEO1rD/kL45HH4ay784wz44ndGQzjrCfjJBpj3OUy5TwKtEMIvSNssbXOvVbkHvv4zzD8L/jwGlv4CUDDjV8boqXmfw9T/goypEmhFr9Qnemp/t/J3bKnY0qXXHB47nJ/l/+yoz2/YsIHc3NxOn3vrrbdYu3Yt69at4/Dhw4wfP55p06YBsGbNGjZu3EhKSgqTJ0/m66+/ZsqUKSdU07333su9997LnDlzePrppz3Hg4ODefvtt4mMjOTw4cNMnDiRiy++mCeeeIINGzawdu1aAJYsWcL27dtZuXIlWmsuvvhili9f7qkNYOPGjWRnZ3sazc6sXr2aDRs2kJGRAcCCBQuIjY2lsbGR8ePHc8UVV9DS0sJjjz1GYWEhUVFRnHXWWYwbN+6Ia82bN4+nn36aIUOG8N1333HnnXfy2WefAVBSUsJXX33Fli1buPjii7nyyit54okneOqpp3j//fdP6N9MdOBywr6VxvzYze9B9T5joYiMaTDpLhh+kSzuJIToEtI2S9ssTlPVXtj4jtEru7/QOJYyzgiyWZdAzEATixOiZ/WJUOtrvvrqK+bMmYPVaqVfv36ceeaZrFq1isjISPLz80lNTQVg7NixFBUVnXDDuWLFCt555x0A5s6dy09/+lPAGG718MMPs3z5ciwWC/v37+fQoUNHvH7JkiUsWbLE04DV1dWxffv2dg1nR7/5zW944403KC0t5cCBAwDk5+d7Gk2Av/zlL7z99tsA7Nu3j+3bt3Pw4EGmT59OQkICANdcc80Rc5rq6ur45ptvuOqqqzzHmpubPfcvvfRSLBYLWVlZnb4fcYKcdij6ygiyWz6AukNgDYTBZ8P0h2DY+TInVgjR60nbLG2zX6jaC5veNYYWtwbZ5LFw7i+NIBubccyXC9Fb9YlQe6xPbbvLyJEjPcNuTkZQUNsy6VarFYfDccQ5ymsPz6ampuNe86WXXqKsrIzCwkJsNhsDBw7s9HVaax566CFuv/32o14rKyuLdevW4XK5sFgsPPLIIzzyyCPthl2FhYV57i9btoxPPvmEFStWEBoayvTp00+oZgCXy0V0dLTn0+qOvP+tWoeCiRPkaIadnxtBduuH0FgJtlAYMsOYHztkJgRHml2lEMIPnejqx9I2S9ssTlDVPq8gawx3JznbWNMi61IJskIgc2q7zdlnn01zc3O7eSrr16/nyy+/ZOrUqbz22ms4nU7KyspYvnw5+fn5J3ztfv36sXnzZlwul+dTVoCJEyfy5ptvAngWwACorq4mMTERm83G559/zp49ewCIiIigtrbWc955553HggULqKurA2D//v2Ulpa2+96ZmZnk5eXx6KOP4nQ6AaPxPlrDVV1dTUxMDKGhoWzZsoVvv/0WgAkTJvDFF19QXl6O3W7njTfeOOK1kZGRZGRkeJ7TWrNu3bpj/tt0fE/CS0u90SguvAWeHAyvXAOb3zcC7DX/gQd2wtUvwugrJdAKIU6ZL69+LG1z2/eWttnHVRfDir/Ds+cauwwseQRcdjjnMbhnDdy+3FjTQgKtEICE2m6jlOLtt9/mk08+YfDgwYwcOZKHHnqIpKQkLrvsMsaMGUN2djZnn302Tz75JElJSSd0TYAnnniCiy66iDPOOIPk5GTP83/605/4wx/+wJgxY9ixYwetf1Bcd911FBQUMHr0aF588UWGDx8OQFxcHJMnT2bUqFE88MADzJw5k7lz5zJp0iRGjx7NlVde2Wkj9Oyzz1JeXu5pRGfMmMGTTz7Zac2zZs3C4XAwYsQIHnzwQSZOnAhAcnIyjz/+OJMmTWLy5MmMGDGi09e/9NJLPPfcc2RnZzNy5EjPQhpHM2bMGKxWK9nZ2bIYBUBTNax/HV69zgiyr/8Adn4GIy+F6xbCAzvg8vkwYjYEhppdrRBCdCtpmw3SNvsoT5CdAX8cCR8/bIysOucx+PFqI8hOvR9iB5ldqRA+R/WGoSF5eXm6dfXBVps3bz7qL2N/NHr0aBYtWtRuPkxHDQ0NhISEoJTi1Vdf5ZVXXjluQ9PX9Lafi07Vl8PWD4ztd3YtMz7ZDU8ygmvWxZB2Blj7xMwDIU6ZUqpQa51ndh3+TNpmg7TNx9fbfi5OSvX+tqHFxSuNY0mjYeRlxtBi2V1ACI9jtc3yl60fmDFjBqNHjz5mowlQWFjI3Xffjdaa6OjoTvenE71UTQlsed+YI1v0NWgnRKfBhNuNhSP658ExVsUUQghxcqRtFqes5kBbkN33nXGs32g4++dGmJUgK8RJ8+tQe6KLUfi7pUuXntB5U6dOPe68FtGLVO4xtt3ZvMjYhgcNcUNgyk+MxZ6Ss8Fr4RIhhBBdR9pmcVJqDhgjqDa+DfuMOcxGkH0Usi6D+N79t6wQ3c2vQ63W+j3gvby8vNvMrkWIHnF4u/Hp7uZFUOL+I6nfaDjrYSPIJg43tz4hhBBCGGpKjDZ70zuwd4VxrN8oOOtRY22L+CFmVidEr+LXoVaIXk9rOLTB+HR383tQttk43j/P2Fx9xGxZMEIIIYTwFTUlxgfPG99xB1kNiSMlyArRzSTUCuFrtDY2VN+8yAizlbtBWYwFns5/EoZfCFGpZlcphBCd6itTg4TwqD3YNrTYE2SzjFFUWZdCwlCzKxSi15NQK4QvcDlh77dGkN38HtTsB0sAZEyDyffC8IsgPMHsKoUQ4rhkapDoE2oPuXtk34Y93wAaEkbA9IeMHtmEYWZXKESfIsuhdqODBw9y7bXXMnjwYHJzc7ngggvYtm3bKV3r1ltvZdOmTV1aX1VVFf/3f/930q+rq6vjRz/6EYMHDyYnJ4fc3FyeeeaZ06rlhRde4O677wbg6aef5sUXXzyl6xQVFfHyyy+fVi09yt4IH/4/+N9h8MIFUPC8scDTpU8be8je8Dbk3SyBVgghuoi0zSeuz7bNR1N7CFY+A89faLTbH/4UGipg+oNw53dw17cw/WcSaIUwgfTUdhOtNZdddhk33ngjr776KgDr1q3j0KFDDB168sNQnn322a4u0dNw3nnnnSf1ultvvZVBgwaxfft2LBYLZWVlnW5R4HA4CAg4+R+xO+6446Rf06q14Zw7d+4pX6PHuFzw1jyjZzbrEmMP2SEzISjC7MqEEKJXkrZZ2uaTVlfqXuzpXSj6CtAQP8wIslmXygKNQvgI6antJp9//jk2m61dI5Cdnc3UqVPRWvPAAw8watQoRo8ezWuvvQbAsmXLmD59OldeeSXDhw/nuuuuQ2sNwPTp02ndxD48PNxzzYULF3LTTTcBsHPnTiZOnMjo0aN59NFHPefV1dVxzjnnkJOTw+jRoz2bvj/44IPs3LmTsWPH8sADDwDw+9//nvHjxzNmzBgee+yxI97Xzp07WblyJb/+9a+xuPc9TUhI4Gc/+5nnPUydOpWLL76YrKwsAC699FJyc3MZOXIk8+fP91zr+eefZ+jQoeTn5/P11197jj/++OM89dRTnu83a9YscnNzmTp1Klu2bAHgpptu4p577uGMM85g0KBBLFy40POevvzyS8aOHcsf//jHk/yv1sM+/aUxdGnmr+Hqf8GoKyTQCiFEN5K2WdrmE+JyweoX4YWL2npk60rhzJ/Bnd/C3SuNUCuBVgif0Sd6ag/+z//QvHlLl14zaMRwkh5++KjPb9iwgdzc3E6fe+utt1i7di3r1q3j8OHDjB8/nmnTpgGwZs0aNm7cSEpKCpMnT+brr79mypQpJ1TTvffey7333sucOXN4+umnPceDg4N5++23iYyM5PDhw0ycOJGLL76YJ554gg0bNrB27VoAlixZwvbt21m5ciVaay6++GKWL1/uqQ1g48aNZGdnexrNzqxevZoNGzZ4NqRfsGABsbGxNDY2Mn78eK644gpaWlp47LHHKCwsJCoqirPOOotx48Ydca158+bx9NNPM2TIEL777jvuvPNOPvvsMwBKSkr46quv2LJlCxdffDFXXnklTzzxBE899RTvv//+Cf2bmabwBfj6T5B3C0y6y+xqhBCix0nbLG2zz1r2W1j+JMQPhWkPuHtkR8je70L4sD4Ran3NV199xZw5c7BarfTr148zzzyTVatWERkZSX5+Pqmpxsq2Y8eOpaio6IQbzhUrVvDOO+8AMHfuXH76058CxnCrhx9+mOXLl2OxWNi/fz+HDh064vVLlixhyZIlngasrq6O7du3t2s4O/rNb37DG2+8QWlpKQcOHAAgPz/f02gC/OUvf+Htt98GYN++fWzfvp2DBw8yffp0EhKMuaLXXHPNEXOa6urq+Oabb7jqqqs8x5qbmz33L730UiwWC1lZWZ2+H5+18zN4/37IPNdYzVgaSSGEMJ20zX28bW614U0j0I67AS7+q7TRQviJPhFqj/WpbXcZOXKkZ9jNyQgKCvLct1qtOByOI85RXr9gm5qajnvNl156ibKyMgoLC7HZbAwcOLDT12mteeihh7j99tuPeq2srCzWrVuHy+XCYrHwyCOP8Mgjj7QbdhUWFua5v2zZMj755BNWrFhBaGgo06dPP6GaAVwuF9HR0Z5Pqzvy/rdqHQrm8w5tgtdvND7xvfJ5sPaJ/wWFEOII0jZL2+xz9q+Gd+40ttC78A8SaIXwIzKntpucffbZNDc3t5unsn79er788kumTp3Ka6+9htPppKysjOXLl5Ofn3/C1+7Xrx+bN2/G5XJ5PmUFmDhxIm+++SaAZwEMgOrqahITE7HZbHz++efs2bMHgIiICGpraz3nnXfeeSxYsIC6ujoA9u/fT2lpabvvnZmZSV5eHo8++ihOpxMwGu+jNVzV1dXExMQQGhrKli1b+PbbbwGYMGECX3zxBeXl5djtdt54440jXhsZGUlGRobnOa0169atO+a/Tcf35FNqD8HL14AtFOa+BsGRZlckhBB9irTNbd9b2uYOakrg1bkQngjX/BsCAs2uSAhxEiTUdhOlFG+//TaffPIJgwcPZuTIkTz00EMkJSVx2WWXMWbMGLKzszn77LN58sknSUpKOqFrAjzxxBNcdNFFnHHGGSQnJ3ue/9Of/sQf/vAHxowZw44dO4iKigLguuuuo6CggNGjR/Piiy8yfLixsEFcXByTJ09m1KhRPPDAA8ycOZO5c+cyadIkRo8ezZVXXtlpI/Tss89SXl7uaURnzJjBk08+2WnNs2bNwuFwMGLECB588EEmTpwIQHJyMo8//jiTJk1i8uTJjBgxotPXv/TSSzz33HNkZ2czcuRIz0IaRzNmzBisVivZ2dm+tRhFSwO8ci00HIa5r0JUqtkVCSFEnyNts0Ha5g7sjUagba6FOa9CWLzZFQkhTpLyu6EhncjLy9Otqw+22rx581F/Gfuj0aNHs2jRonbzYTpqaGggJCQEpRSvvvoqr7zyynEbmr7GlJ8LlwtevwG2fADXvgzDL+jZ7y+EOGlKqUKtdZ7ZdfgzaZsN0jYfn6k/F1rDm7cac2mljRbCpx2rbfbJCX1KqTDgC+BxrbWfLJXXfWbMmMHo0aOP2WgCFBYWcvfdd6O1Jjo6utP96YQJPvkFbHkfZj0hjaUQotdTSs0GZmdmZppdSreStrmX+PJ/YcNCOPdxaaOF8GM9EmqVUguAi4BSrfUor+OzgD8DVuBZrfUT7qd+BrzeE7X5g6VLl57QeVOnTj3uvBbRwwoWwDd/hfG3wYRT37heCCH8hdb6PeC9vLy828yupTtJ29wLbH4PPvtvGHMNTP6J2dUIIU5DT82pfQGY5X1AKWUF/g6cD2QBc5RSWUqpGcAmoLTjRYTwKzs+gQ9+CkNmGr20soqiEEII4RtK1sNb86B/Hsz+i7TRQvi5Hump1VovV0oN7HA4H9ihtd4FoJR6FbgECAfCMIJuo1LqQ6216xS/b7sl9kXf1qPzxw9thNdvgsQsuHKBbN0jhBBu0jYLb6as7VJXCq/MgZAYYx6tLbjnaxBCdCkz/9LuD+zzelwMTNBa3w2glLoJOHy0QKuUmgfMA0hLSzvi+eDgYMrLy4mLi5PGU6C1pry8nODgHmi4ag/CS1dDULixdU9QRPd/TyGE8APSNgtvPdo2t3I0w2vXQ0M5/HAxRPTrue8thOg2Ptt9pLV+4TjPzwfmg7HCYsfnU1NTKS4upqysrHsKFH4nODiY1NRu3kqnpd7YuqexEm7+EKL6d+/3E0IIPyJts+ioR9rmVlrD+/fBvu/gqhcgZWzPfF8hRLczM9TuBwZ4PU51H+sSNpvtuCsSCtGlXE548zYoWQfXviKNpRBCdCBtszDVir/B2pdg+kMw8jKzqxFCdKGeWiiqM6uAIUqpDKVUIHAtsOhkLqCUmq2Uml9dXd0tBQpxUpb+ArZ+YCwKNWzW8c8XQgghRM/YtgSW/ByyLoVp/8/saoQQXaxHQq1S6hVgBTBMKVWslLpFa+0A7gY+BjYDr2utN57MdbXW72mt50VFRXV90UKcjFXPGp8A598OE243uxohhBBCtCrdAgt/CMlj4NJ/gMXMPh0hRHfoqdWP5xzl+IfAhz1RgxDdZvtS+PABGDoLZv3W7GqEEEII0aqhAl65BgJDjalBgaFmVySE6AZ+/VGVDD8Wpjv4PbxxE/QbBVc8Bxar2RUJIYQQAsBph9d/ADUlxtY9snijEL2WX4daGX4sTFVTAi9fA0GR7q17ws2uSAghhBBgrHT84QNQ9CVc/FdIzTO7IiFEN/LZLX2E8Gkt9cZwpsYqY5+7yBSzKxJCCCFEq1XPQuHzMOU+yL7G7GqEEN3Mr3tqZfixMIXLCW/eagw9vup5Y+EJIYQQQviGnZ/DRz+DoefD2b8wuxohRA/w61Arw4+FKZY8Cls/hPOfhKHnmV2NEEIIIVqV74Q3boSEYXDFM7LSsRB9hPyfLsTJ+G4+fPt/MOFHkH+b2dUIIUSPUUoNUko9p5RaaHYtQnSqscpY68ISAHNegaAIsysSQvQQCbVCnKhtH8Ni93Cm835jdjVCCHHClFILlFKlSqkNHY7PUkptVUrtUEo9eKxraK13aa1v6d5KhThFToexF21lEVz9b4gZaHZFQoge5NehVubUih5Tst5oLJNGwxXPytY9Qgh/8wIwy/uAUsoK/B04H8gC5iilspRSo5VS73f4Suz5koU4CUt/Djs/hYv+AAMnm12NEKKH+XWolTm1okfUHDCGMwVHwRzZukcI4X+01suBig6H84Ed7h7YFuBV4BKt9fda64s6fJX2eNFCnKjCfxlTgybeCTk/MLsaIYQJ/DrUCtHtmuvg5auhuQbmvg6RyWZXJIQQXaU/sM/rcbH7WKeUUnFKqaeBcUqph45yzjylVIFSqqCsrKxrqxWiM0Vfwwf/BYPPgRn/bXY1QgiTyD61QhyNywlv3gKHNhqBNmmU2RUJIYRptNblwB3HOWc+MB8gLy9P90Rdog+rLILXbzDmz165AKzyZ60QfZX83y/E0Xz8MGxbDBf+LwyZYXY1QgjR1fYDA7wep7qPCeH7mmvhlTnGB9BzX4OQaLMrEkKYyK+HH8tCUaLbfPs0fPc0TLobxt9qdjVCCNEdVgFDlFIZSqlA4Fpg0eleVNpm0e1cTnjzNijbClf/C+IGm12REMJkfh1qZaEo0S22fgQfPwTDLoQZvzK7GiGEOG1KqVeAFcAwpVSxUuoWrbUDuBv4GNgMvK613ni630vaZtHtPvtv2PYRnP87GDTd7GqEED5Ahh8L4a1kHSy8BZLGwBXPyNY9QoheQWs95yjHPwQ+7OFyhDh1616Dr/4IebdA/m1mVyOE8BF+3VMrRJeq3m9s3RMSY8zPCQwzuyIhhPA7MvxYdJt9q2DRj2HgVKOXVggh3CTUCgHGghMvX2Ns4XPd6xCRZHZFQgjhl2T4segW1cXw6lyITIGrXwSrzeyKhBA+RIYfC+F0wMIfQukmI9D2G2l2RUIIIYRo1VJvrHTsaIKb3ofQWLMrEkL4GL/uqZUhTuK0aQ2LH4TtS+DCpyDzXLMrEkIIIUQrlwvevgMObTD2ok0YZnZFQggf5NehVoY4idP23dOw6hk448eQ90OzqxFCCL8nHziLLvXF72DzIpjx37JnvBDiqPw61ApxWrZ8AIsfguEXwbmydY8QQnQF+cBZdJkNb8EXT8DY62HSXWZXI4TwYRJqRd90YA28eSukjIPLnwGL/K8ghBBC+IwDa+CdO2HARLjoD6CU2RUJIXyY/CUv+p7qYnj5WgiNgzmvQmCo2RUJIYQQolXtQXhlLoTFwzX/gYAgsysSQvg4CbWib2mqgZeuBnsDzH0dIvqZXZEQQvQqMqdWnBZ7o7F1T1M1zHkFwhPMrkgI4Qck1Iq+w+mAhTdD2Ra4+l/QL8vsioQQoteRObXilGkNi+6B/YVw+XxIGm12RUIIPyH71Iq+QWv46AHY8QnM/jMMPtvsioQQQgjh7as/wvevw9k/hxEXmV2NEMKP+HVPrQxxEidsxd+hYAFMvhdybzK7GiGEEEJ42/IBfPorGHUlTP0vs6sRQvgZvw61MsRJnJDN78OSR2HExXDO42ZXI4QQQghvBzfAm7cZOxJc8jdZ6VgIcdL8OtQKcVz7Vxtb9/TPMebnyNY9QgjRrWQUlTgp9YfhlTkQHAnXvgy2ELMrEkL4IfkLX/ReVXvhlWshLMHYukcaSiGE6HYyikqcMEcLvHY91JcagTYy2eyKhBB+ShaKEr1TUzW8fA3Ym+AHiyA80eyKhBBCCNFKa/jgPti7Aq5cYIyoEkKIUyShVvQ+Tju8cRMc3gbXLYTE4WZXJIQQQghv3/4D1vwHpv0/GHWF2dUIIfychFrRu2gNHz4AOz+D2X+BwWeZXZEQQgghvG3/BJY8AiNmw/SHzK5GCNELyJxa0bt881cofB6m3Ae5N5pdjRBCCCG8lW2DhTdD4ki47J+ygKMQokvIbxLRe2xaBEt/AVmXwtm/MLsaIYTok2T1Y3FUDRXwyjUQEARzXoHAMLMrEkL0EhJqRe9QXAhvzYPUPLjsafnkVwghTCKrH4tOta53UV1srHQcPcDsioQQvYjMqRX+r3KP8clveAJc+4ps3SOEEEL4msUPwe4v4NKnYUC+2dUIIXoZv+7OkiFOgsYqePlqY6+76xYawVYIIYQQvmPVs7DqGZh8L4ydY3Y1QoheyK9DrQxx6uOcdnjjRijfAdf8GxKGmV2REEIIIbzt+gI+/H8wdBac85jZ1QgheikZfiz8k9bwwf2waxlc8ncYdKbZFQkhhBDCW/lOeP0HED8ULn8GLFazKxJC9FJ+3VMr+rCv/wyrX4Sp/wXjrje7GiGEEEJ4a6qGV64FZTFWOg6ONLsiIUQvJj21wv9sfAc+eQxGXg5nPWp2NUIIIYTw5nLCwh9CxS74wbsQm2F2RUKIXk56aoV/2bcK3r4dUvPh0n/I1j1CCOFjZBFHwdJfwI5P4ML/hYFTzK5GCNEHSCIQ/sPlhDd/COH9jKFMtmCzKxJCCNGBLOLYx63+N6z4G0y4A3JvMrsaIUQfIcOPhf/Y9jFU7YWrX4SweLOrEUIIIYS3PSvg/ftg0Fkw8zdmVyOE6EOkp1b4j1XPQEQyDLvQ7EqEEEII4a1yD7x2PcSkw1XPg1X6TYQQPUd+4wj/UL4Tdn4G0x+WhlII0ScppSxANpACNAIbtNal5lYlBNBcB6/MAZcd5rwGITFmVySE6GMkHQj/sOo5sARA7o1mVyKEED1KKTUY+BlwLrAdKAOCgaFKqQbgn8C/tNYu86oUfdqiH0PZFrh+IcRnml2NEKIPklArfF9LA6z9D4y4GCKSzK5GCCF62q+BfwC3a6219xNKqURgLnAD8C8TahN9XflO2PgWTHsABp9tdjVCiD5KQq3wfRsWGpu4j7/V7EqEEKLHaa3nHOO5UuBPPVeNEB0UPm+MpJI2WghhIlkoSvg2rWHlM5CYBelnmF2NEEKYRil1lVIqwn3/50qpt5RSOWbXJfowexOseQmGXygjqYQQpvK5UKuUGqGUeloptVAp9SOz6xEmKy6Ag+th/C2glNnVCCGEmX6uta5VSk0BzgGewxiWLIQ5Ni+CxgrI+6HZlQgh+rgeCbVKqQVKqVKl1IYOx2cppbYqpXYopR4E0Fpv1lrfAVwNTO6J+oQPW/UMBEbAmGvMrkQIIczmdN9eCMzXWn8ABJpYj+jrChZA7GAYOM3sSoQQfVxP9dS+AMzyPqCUsgJ/B84HsoA5Sqks93MXAx8AH/ZQfcIX1R+GjW/D2DkQFGF2NUIIYbb9Sql/AtcAHyqlgvDNEVezlVLzq6urzS5FdKdDm2DvCsi7GSw+92MohOhjeuS3kNZ6OVDR4XA+sENrvUtr3QK8ClziPn+R1vp84LqeqE/4qNUvgrMF8m4xuxIhhPAFVwMfA+dprauAWOABUyvqhNb6Pa31vKioKLNLEd2p8HmwBkH2XLMrEUIIU1c/7g/s83pcDExQSk0HLgeCOEZPrVJqHjAPIC0trduKFCZxOaHgeRg4FRKHm12NEEKYRikV6/VwmdexZqDAjJpEH9dSD+tehZGXQlic2dUIIYTvbemjtV6Gu9E+znnzgfkAeXl5+jinC3+zfQlU74WZ/212JUIIYbZCQAMKSAMq3fejgb1AhmmVib5pw5vQXCMLRAkhfIaZkyD2AwO8Hqe6jwlhbOMTkWxsEyCEEH2Y1jpDaz0I+ASYrbWO11rHARcBS8ytTvRJBQuMrfYGTDC7EiGEAMwNtauAIUqpDKVUIHAtsOhkLiCLUfRS5Tth56eQexNYbWZXI4QQvmKi1tozLUdr/REgG3iLnnVgjfGV90PZak8I4TN6akufV4AVwDClVLFS6hattQO4G2PRi83A61rrjSdzXVmMopcqWACWACPUCiGEaHVAKfWoUmqg++sR4IDZRYk+puB5sIXCmKvNrkQIITx6ZE6t1nrOUY5/iGzbI7y1NMCa/8CI2RCRZHY1QgjhS+YAjwFvux8vdx8Tomc0VcP3C2HUFRAsHQpCCN/hcwtFnQyl1GxgdmZmptmliK6y4U1oqoLxt5pdiRBC+BStdQVwr9l1iD5s/etgr5cFooQQPsevd8uW4ce9jNaw6hlIGAHpk82uRgghfIpSaqhSar5SaolS6rPWL7PrEn2E1sbQ4+Sx0D/H7GqEEKIdv+6pFb3M/kIoWQcXPCWLTwghxJHeAJ4GngWcJtci+pp9K6F0I8z+i9mVCCHEEfw61Mrw415m5TMQGAHZ15pdiRBC+CKH1vofZhch+qiCBRAUacynFUIIHyPDj4VvqC+HjW8ZgTYowuxqhBDCF72nlLpTKZWslIpt/TK7KNEHNFTAxrdhzDUQFG52NUIIcQS/7qkVvciaF8HZAuNvMbsSIYTwVTe6bx/wOqaBQSbUIvqStS+Dsxnybja7EiGE6JSEWmE+l9MY1jRwKiSOMLsaIYTwSVrrDLNrEH2Q1kYbPWAi9BtpdjVCCNEpvx5+LHqJ7Uuhaq/00gohxDEopWxKqXuUUgvdX3crpWxm1yV6ud3LoWKnbOMjhPBpfh1qlVKzlVLzq6urzS5FnI5Vz0BEMgy/yOxKhBDCl/0DyAX+z/2V6z7WI5RSlyqlnlFKvaaUmtlT31eYrGABhMRA1iVmVyKEEEfl16FWForqBSp2wY5PIPcmsEqHgxBCHMN4rfWNWuvP3F83A+NP5IVKqQVKqVKl1IYOx2cppbYqpXYopR481jW01u9orW8D7gCuOeV3IfxH7SHY8j6MvQ5swWZXI4QQR+XXoVb0AqueA0sA5Nx4/HOFEKJvcyqlBrc+UEoN4sT3q30BmOV9QCllBf4OnA9kAXOUUllKqdFKqfc7fCV6vfRR9+tEb7f2P+ByQK4sECWE8G2yUJQwj70R1vzHGHYcmWx2NUII4eseAD5XSu0CFJAOnFDa0FovV0oN7HA4H9ihtd4FoJR6FbhEa/1b4Ij5IEopBTwBfKS1Xt3Z91FKzQPmAaSlpZ1IacJXuZxQ+AJkTIP4TLOrEUKIY/LrUKuUmg3MzsyUX7Z+acOb0FQF4281uxIhhPB5WutPlVJDgGHuQ1u11s2nccn+wD6vx8XAhGOc/2PgXCBKKZWptX66kxrnA/MB8vLy9GnUJsy28zNjEccZvzK7EiGEOC6/Hn4sc2r9mNaw8hlIGAEDp5hdjRBC+Dyl1F1AiNZ6vdZ6PRCqlLqzp76/1vovWutcrfUdnQVa0csULICwRBh2odmVCCHEcfl1qBV+bP9qKFlrbOOjlNnVCCGEP7hNa13V+kBrXQncdhrX2w8M8Hqc6j52WmRngl6guhi2LYacGyAg0OxqhBDiuCTUCnOsegYCw2GMLKAphBAnyOqe1wp4Fno6ncSxChiilMpQSgUC1wKLTrNGGUXVG6x+0RhRJYs4CiH8hIRa0fPqy2HDW5B9LQRHml2NEEL4i8XAa0qpc5RS5wCvuI8dl1LqFWAFMEwpVayUukVr7QDuBj4GNgOva603dlPtwl847VD4LxgyA2LSza5GCCFOiF8vFCX81Jp/g7NZFogSQoiT8zPgduBH7sdLgWdP5IVa6zlHOf4h8GGXVCd6h22Loe4g5P3J7EqEEOKE+XWoldWP/ZDLaSw+kT4FEkeYXY0QQvgNrbVLKfUC8JnWeqvZ9RyNtM1+rmABRKbCkJlmVyKEECfMr4cfy7wdP7TjE6jaYywQJYQQ4oQppS4G1uIecqyUGquUOu05sF1N2mY/VrHL2Mon90awWM2uRgghTphfh1rhh1Y+A+FJMGK22ZUIIYS/eQzIB6oAtNZrgQwT6xG9TeELoKww7gazKxFCiJMioVb0nIpdRk9t7k1gtZldjRBC+Bu71rrjPjnalEqOQbb08VOOZljzHxh+AUQmm12NEEKcFAm1oucULABlMYY1CSGEOFkblVJzMbb2GaKU+ivwjdlFdSTDj/3U5vegoRzyfmh2JUIIcdIk1IqeYW80PgEecRFEpphdjRBC+KMfAyOBZoztfGqAn5hZkOhFCp6HmAzImG52JUIIcdL8evVj4Uc2vAWNlbKNjxBCnCKtdQPwCPCIUsoKhGmtm0wuS/QGZVthz1dw7i/BIv0dQgj/I7+5RM9Y9SwkDIeBU82uRAgh/JJS6mWlVKRSKgz4HtiklHrA7Lo6kjm1fqjgebDYYNz1ZlcihBCnxK9DrTScfmJ/IRxYbfTSKmV2NUII4a+ytNY1wKXARxgrH/vcMrUyp9bPtDTAupch6xIIize7GiGEOCUnHGqVUr/u5Jipm5hJw+knVj4LgeEw5hqzKxFCCH9mU0rZMELtIq21HR9c/Vj4mY1vQ1O1LBAlhPBrJ9NT218pNaf1gVIqEfik60sSvUpDBWx40wi0wZFmVyOEEP7sn0AREAYsV0qlYywWJcSpK1gA8cMg/QyzKxFCiFN2MgtF3Q58rJTaifHJ8PPAz7qlKtF7rPk3OJtlgSghhDhNWuu/AH9pfayU2gucZV5Fwu+VrIP9BTDrdzI9SAjh144bapVSLwKrgTXAXcDLgAO4VGu9o3vLE37N5YJVz0H6ZOiXZXY1Qgjhl5RS1wMva61d3se11hpwKKUGA8la669MKbADpdRsYHZmZqbZpYjjKXgeAkIgW6YHCSH824n01L4AZAM3A2OAgcAq4Hql1Aat9cJuq074tx2fQNUeOPdxsysRQgh/FgesUUoVAoVAGRAMZAJnAoeBB80rrz2t9XvAe3l5ebeZXYs4huZa+P4NGHUFhMSYXY0QQpyW44ZarfVnwGetj5VSAcAIjKA7AZBQKzq36hkI7wfDLzK7EiGE8Fta6z8rpf4GnA1MxviAuRHYDNygtd5rZn3CT61/HVrqZIEoIUSvcDJzagHQWjsw9sf7HvhPl1ckeoeK3bB9KZz5/yAg0OxqhBDCr2mtncBS95cQp0drY+hx0hjon2N2NUIIcdr8ep9a4cMKFoCyQO5NZlcihBBCCG/7C+HQ90YvrSwQJYToBSTUiq5nbzRWPR5+IUSmmF2NEEIIIbwVLIDACBh9pdmVCCFEl5BQK7rexrehsVK28RFCCCF8TWOle//4qyAowuxqhBCiS/h1qFVKzVZKza+urja7FOFt1bPGRu4Z08yuRAgheg2lVD+l1HNKqY/cj7OUUreYXVdH0jb7uHWvgqNJFogSQvQqfh1qtdbvaa3nRUVFmV2KaLV/tTFXZ/ytMk9HCCG61gvAx0DrvI5twE/MKuZopG32YVobQ49Tx0PSaLOrEUKILuPXoVb4oFXPgi1MNnIXQoiuF6+1fh1wgWc3Aqe5JQm/sudrOLxNemmFEL2OhFrRdRoqjHk62ddAsHxCL4QQXaxeKRUHaACl1ERAxviKE1ewwGifR15mdiVCCNGlTnqfWiGOas1/jHk6skCUEEJ0h/uBRcBgpdTXQAIgy9eKE1NXBpsWQf5tYAsxuxohhOhSEmpF13C5oOA5SDsD+o00uxohhOh1tNarlVJnAsMABWzVWttNLkv4i7X/AZcdcm82uxIhhOhyEmpF19j5KVQWwTm/MLsSIYTolZRSVuACYCBG+z1TKYXW+g+mFiZ8n8sFBc/DwKmQMNTsaoQQostJqBVdY+UzEJYIw2ebXYkQQvRW7wFNwPe4F4vyRUqp2cDszMxMs0sRrXZ9BlV74NzHzK5ECCG6hYRacfoqi2D7Epj2AAQEml2NEEL0Vqla6zFmF3E8Wuv3gPfy8vJuM7sW4VbwPITGywfPQoheS1Y/FqevYAEoC+TeZHYlQgjRm32klJppdhHCz9QcgK0fQc4N8sGzEKLXkp5acXrsTbD63zD8Aojqb3Y1QgjRm30LvK2UsgB2jMWitNY60tyyhE9b/W/QLsi50exKhBCi20ioFadn49vQWAHjZZSZEEJ0sz8Ak4Dvtdba7GKEH3A6YPW/YPDZEJthdjVCCNFtZPixOD2rnoX4oZAxzexKhBCit9sHbJBAK07Y9iVQsx/yfmh2JUKIPqrF4cLl6v5mS3pqxak7sAb2F8D5T4JSZlcjhBC93S5gmVLqI6C59aBs6SOOqmABRCTD0FlmVyKE6AXsTheVDS1UNdipqG+hsr6FigbjtrLB3unj2mYHyx84i7S40G6tzSdDrVLqUuBCIBJ4Tmu9xNyKRKdWPgu2MMi+1uxKhBCiL9jt/gp0fwlxdJVFsOMTOPNnYPXJP/eEECZyOF1UNtipamgxAmpDCxX1diob2ofVCq9zapscR71eeFAAMWE2YkIDiQkNZFBCuPu+jbAga7e/nx77LaeUWgBcBJRqrUd5HZ8F/BmwAs9qrZ/QWr8DvKOUigGeAiTU+pqGCtiwELLnQHCU2dUIIUSvp7X+pdk1CD9S+C9jFFXOD8yuRAjRzRxOF9WNdk8wbQ2pnoBa3/pciyeg1hwjoIYGWokJDSQ2LJCYsEAy4kKJ9nocGxroCbCxYYFEh9oICuj+4HosPfnR3QvA34AXWw8opazA34EZQDGwSim1SGu9yX3Ko+7nha9Z+xI4mmD8rWZXIoQQvZpS6k9a658opd4DjpiYpLW+2ISyhC9ztMCaf8PQ82VnAiH8jNOlqW60e/WetgbRtmBa6QmtxnnVjfajXi/EZnWHUSOEpsWGGo9DA4kNsxHjvu8dUINt5gbUU9FjoVZrvVwpNbDD4Xxgh9Z6F4BS6lXgEqXUZuAJ4COt9eqeqlGcIJcLVj0HaZMgadTxzxdCCHE6/u2+fcrUKoT/2PI+1JfJAlFCmEBrTbPD6DmtabRT02SnptFhPG4yjhnPOYzHTW2PW8852nKAQQEW4sICPb2m/WNCiQ31CqZevaitwdUfA+qpMHuSRX+M1RxbFQMTgB8D5wJRSqlMrfXTHV+olJoHzANIS0vrgVKFx87PoHI3nP2o2ZUIIUSvp7UudN8dq7X+s/dzSql7gS96virh0woWQHSasZWPEOKk2Z0udyB1tIXQTsJpTZOjQ3g1zmlxuo55/RCblciQACKDbUSG2EiMCCYzIYDIEJsRWN1BNbZDWA0J7BsB9VSYHWo7pbX+C/CX45wzH5gPkJeXJ9sb9KRVz0JYIoyQEW9CCNGDbsRYg8LbTZ0cM5VSajYwOzMz0+xS+qaybVD0JZzzGFhk50bRN7lcmtrmIwNpW/A8MpB69542tDiPef0AiyIqxAikkSE2IoMD6B8TYhwLtnkCa5TX8633I4IDTJ9/2huZHWr3AwO8Hqe6jwlfVbkHti2GaT+FAFl8UwghuptSag4wF8hQSi3yeioCqDCnqqPTWr8HvJeXl3eb2bX0SYUvgMUG4643uxIhukSzw0lFfQvldS2U1TVTXtdCeV0z5e65pZ31oNY2O446hBeMNdQiggKICnWH0GAbGfFhXqG0NZB2DKfGsRCbFSXbWfoUs0PtKmCIUioDI8xei9FwnxD5NNgEBQuM3wS5N5ldiRBC9BXfACVAPPC/XsdrgfWmVCR8k70R1r0MI2ZDeKLZ1QjRKZd7IaTy+mYO1xlhte1+M4dbg2t9C4frmo+6jUxQgIWY0EBP+EyJDmZ4cES73tPWINouoIbaCA8MwGKRUNqb9OSWPq8A04F4pVQx8JjW+jml1N3Axxhb+izQWm880WvKp8E9zN5krKY47AKISjW7GiGE6BO01nuAPcAks2sRPm7Tu9BYKQtEiR7XZHd6hdFmDte2cLi+fa9qa2itqG/B4TqyG1UpiA0NJC48kLiwIEamRBIfHkRcWCDxEcZtXHgQ8eHGbVig9JaKNj25+vGcoxz/EPiwp+oQp2HTO9BQDvnyGYIQQvQ0pdTlwO+AREC5v7TWOtLUwoTvKFgAcUNg4BSzKxF+zunSVDW0BdHD9e5w6tWrethrKHD9UeaghgZaPSG1f3Qw2alRnsdx4YFGaHU/jg0LxCq9p+IUmT38+LTI8OMetupZo7HMONPsSoQQoi96Epittd5sdiHCBx3cAPu+g/P+x+jyEqKD+mYH5XVH9qCW1Rq33qG1or6FTjpTsSiIDTN6S+PDg0hLC/UKqEeG1dBAv44awo/49U+aDD/uQQfWQvEqmPU7aSyFEMIchyTQiqMqfB6sQZDd6cA40Us12Z2U1TZTVtdMWa0xH9X71rhvBNdGe+e9qRFBAUZvaXgQ6XGh5KTHuANq67DftiG/0SE2mYsqfJJfh1rRg1Y9A7ZQyL7W7EqEEKKvKlBKvQa8AzS3HtRav2VaRcI3NNfButdg1OUQGmt2NeI0tThcHK7rGEzbB9SyumYO1zZT29z5IkoxoTbiw4NIiAhiXFo08eFB7i+vIb/u+arBNtleRvg/CbXi+Bor4fuFRqANiTa7GiGE6KsigQZgptcxDUio7es2LISWWlkgyoc5nC4q6lso9e5FrTMWVGoNqK29rdWN9k6vEREcQEKEEU6zUiJJcIfWhPAg4iMCSQgPJiHCmJsaGCB7FIu+xa9Drcyp7SFrXgJHE4y/1exKhBCiz9Ja32x2DcIHaQ2rnoN+oyB1vNnV9Ckul6aioeWYQ35bj1U0tHS6b2pYoJV4dzAdkhjOpEFxnuBq3AZ6HkuPqhBH59ehVubU9gCXCwqegwETIWm02dUIIUSfpZR6HqNnth2ttXTP9WUHVsPB9XDh/8qaF12kscXJ/qpGSmub2g/57TBntby+BWcnqykFBVg8QXRAbOsc1dZe1UD3bTDxEbKQkhBdRf5PEse26zOo2AVnPWJ2JUII0de973U/GLgMOGBSLcJXFCwAWxiMvtrsSvxGi8PFgapG9lU2sK+ikeLKBvZVNrKvooHiykYO1zUf8RqbVXmCaXJUMKP7R3n1pAZ7elQTIoIIDwqQ/VOF6GESasWxrXoOwhJgxGyzKxFC9EEtDhc1TXZqGu3UNDmobmy9b3ffd3ier2608+SVY0iOCjG77G6htX7T+7FS6hXgK5PKEb6gsQq+fxOyr4Fg2a64lcPp4mBNU7vAWuwOrPsqGzhY09RuKHCARZESHcKA2BDOGZ7IgNgQUmNCSYwMItHd4xoVYpOgKoQP8+tQK3Nqu1nVXti2GKbcDwFBZlcjhPBDTpemrjWMeoVP477DK5waobXj80fbgqKVzaqICrERGWwjMsRGk93VQ+/MJwwBEs0uQpho/WvgaOxzC0S5XJqyumZPz+q+igb2VbaF1pKqJhxew4KVguTIYFJjQpk0OI4BMaGkxoQwIDaUAbGh9IsIIsAqCysJ4c/8OtTKnNpuVrDAuM29ydQyhBDm0VpT3+JsC5tHCZ+dhdOaRvtRt5toZVEQ6QmlAUSF2MiMCG/3uPV5436AJ8BGhdgICrD0md4TpVQt7efUHgR+ZlI5wmxaQ8Hz0D8XkrPNrqZLaa2pqG/xhNQjelyrGmlxtP8Ay5i/GsK4ATHMHuMOrO7wmhIdIqsBC9HL+XWoFd3I0QyrX4RhF0D0ALOrEUKcIq01jXanJ3jWNrUfsnu8Ib01TY5OF0LxFh5khM+I4AAiQ2z0jw4hKznSE0DbgmmAJ4y2Ppa5ZydGGf9II7XWe036/iOAe4F44FOt9T/MqEN42fstlG2GS/5udiWnpKbJbvSwugNrsdec1n2VDTS0tB+lER1qY0BMKMOSIjg3qx8DYowhwgNiQ+gfHUpIoKwMLERfJqFWdG7jO9BQLtv4CGEyl0tT1+Kg1qv3s6bJ4Q6nXr2iTXbjnA6htbbJ0W4YXmeCbZZ24TM+PJBBCWFew3o7htO2XtPwoAAZttcDtNZaKfUBcNLL0CulFgAXAaVa61Fex2cBfwaswLNa6yeO8f03A3copSzAi4CEWrMVLICgKBh5udmVdKqhxUFxpbuHtaJ9YN1X0UBNU/tRHGGBVs9w4DMy44zA6h4inBoTQkSwzaR3IoQ4FQ32BlYeXMmKAyv4ad5PsVm79/9hvw61Mqe2G616FuIyIeNMsysRwq85nC7qmh1tQbOT0Hn0Y8bw3c72NvQWGmj1hM+I4LZQ6n3MO5xGBBvHWoNpUID0cPiJ1Uqp8VrrVSf5uheAv2GEUQCUUlbg78AMoBhYpZRahBFwf9vh9T/UWpcqpS4GfgT8+xTrF12lvhw2vQO5N0NgqCklOF2akupG9pY3UFTe4Amr+yob2V/ZwOG6lnbnBwVYPPNYx6VFMyAm1BNYB8SEEh0qCzEJ4c+01myv2s7X+7/m6/1fU1haiMPlICQghMuHXM6w2GHd+v39OtTKnNpuUrIOilfCrCfAIj0wom9zuTTVjXaqvIbo1jS6e0q9wmitV4+p97G648wpBYgIMobltg7fTYkOYXhwRLshu61h1Hv+aYT7mE16StFas7tmNxmRGb35D+MJwHVKqT1APaAwOnHHHOtFWuvlSqmBHQ7nAzu01rsAlFKvApdorX+L0avb2XUWAYvcPcYvn9Y7Eadn7UvgbIG8m7v127Q4XOyrbHAH13r2lDewt8K4X1zRSIuzbV5rgEXR3x1Qzx3RzxNYW4cIJ4QH9eb/N4Xok6qbq/m25FtPkC1tLAUgMzqT60dcz+T+k8lJzCHQGtjttfh1qBXdZNWzYAuF7DlmVyJEl3O5NFWNdirqmzlc10JFfQvl9S1U1LVQXt/suW8cb6aywX7MOaXeCx1FBBs9oQPjQz2LGbUe8w6o3sfCgwKwWuQPvVNR1lDGtyXfGl8HvqW0sZR3L32XQVGDzC6tu5zXhdfqD+zzelyMEZo7pZSaDlwOBAEfHuO8ecA8gLS0tC4oUxzB5YLC5yHtDEgccdqXa2hxsKe8gT3u0FpU3sDeinqKDjdQUt2I96+/0EAr6XFhDE2MYEZWP9JjwxgYF0paXCjJUSHyu0yIXs6lXWwq38RX+7/i6/1fs/7welzaRURgBJOSJzGl/xQmpUwiKSypx2uTUCvaa6yE9W/AmKshJNrsaoQ4LqdLU9XQFk7L61qoaA2n7sfl9c2e+5UNLRwto0aF2IgLCyQ2LJD0uFBy0mM8j2PCbF7htK3XNCzQKr0PPaTB3kDBoQK+LfmWFQdWsKNqBwDRQdFMSJ7ApORJxAXHmVxl99Fa7zHxey8Dlp3AefOB+QB5eXnHGTgvTsnuL6BiF0x/+IRfUtXQQpFXcPWE2IoGymqb250bE2ojLS6MvIExpMf2Jz0ujPS4UNLjwogPD5Tfd0L0MYcbD7PiwAq+2v8VKw6soLK5EoViZNxIbht9G1P6T2FU/CgCLObGSgm1or21Lxt73skCUcIkTpemsqHzQNrae9p6v6L+BEJqeCBxYYFkxIeRmx5LfLgRUmPDAokLC/I8HxMWKMN4fYzT5WRj+UZWHFjBtyXfsrZsLQ6Xg0BLIDn9crho0EVMSpnE8NjhWJT8tztJ+wHvpe1T3cdOm6x30c0KFkBoHGRd7Dmktaa0tpmiw0ZQ3eM9VPhw/RGLMiVFBpMWF8pZwxJIjwsjLTaUgXFhpMWFEhUiCzIJ0ZfZXXbWla7jmwPf8NX+r9hcsRmA2OBYpvSfwuT+k5mUMonY4FiTK21PQq1o43LBqudgwARIPuYULSFOmMPporLB7gmkbWHV6FGtqG/xDANuDalHWxgpOtTmDqOBDE4IZ3yGcT8uLJDY8CBPr2pceCAxoRJS/Y3Wmn21+1hxYAUrSlaw8uBKaltqARgRO4Ibsm5gUvIkxiWOIzgg2ORq/d4qYIhSKgMjzF4LzO2KC8t6F93D4XRxqLiIlC0fsGngDbzz8U5jqHB5A3sq6mmyt81vtVoU/aNDSI8L5eKxKUZgjQ1lYHwYA2Jk+xshRHsldSV8dcAYUvxdyXfU2euwKivZCdncM+4eJvef7PMfIPt1qJVPg7vYrs+hYidMf8jsSoSP01pT0+SgrLaJQzXNlLbe1jRzqLaJstpmyuuMwFrVaD9qSI3xhNQgMhPCiXOHVCOYukOqu2dVQmrvVNlUyXcl33mGFB+oPwBAclgyM9JnMCl5EvnJ+T73ibA/UUq9AkwH4pVSxcBjWuvnlFJ3Ax9jrHi8QGu90cQyBdBkd1Jc2UDR4YZ2Pa57yusprmzkDvUWP7U5uXPLGEqse0iPDSU9LpQpQ+I9Q4TTY0PpHxMivy+FEEfV7Gym8GChJ8juqt4FQFJYEucNPI8p/acwIXkCEYERJld64pQ+3l4RfiAvL08XFBSYXYb/e2Uu7PsO7t8EAUFmVyNMoLWmptHBodomI6DWNFFaa9yW1bY9Lq1tatcr0Co00Eq/yGASwo1hve3CqbtHNS48yB1SbbK/aR/U7Gxm9aHVrChZwbcHvmVLxRY0mghbBPnJ+UxKnsTElImkRaSZOndPKVWotc4zrYBeQNrmzjXZneworfOsIty6svDe8gZKaprafQgYERRAWpwxNDg9Noi71l2OPSaTxmsX0i8iGIsszCSEOAFaa/bU7OHrA1/z1f6vKDhYQJOziUBLILn9cpncfzJT+k9hUNQgn543f6y22a97akUXqtoH2z6CKfdJoO2FtDa2pWnXq+oOrh0fNzuODKvhQQEkRgSRGBnE2AHR9IsMIjEimET3bb/IIBIjgwkPkl8poj2XdrG1YisrSlaw4sAK1pSuodnZTIAlgOyEbO4aexeTUiaRFZdl+iITomvIKKo2WmuKKxtZvbeSNXurWLO3kk0lNdidbck1PjyQtNhQJg6KIy0utF2Pa2yY18JMWxfDdwdh6u+Jjgox6R0JIfxFvb2elSUrPUF2f52xbMLAyIFcMfQKJqdMJi8pj5CA3vH7RP6CEIaCBcZtbvfueSe6ltaaqgb7ET2rpV49rEbPajMtnYTViKAATzDNTYshMTLYHV6N237u2zAJq+IkHKg74BlO/F3Jd1Q2VwLGvnVXDb2KSSmTyOuXR6gt1ORKRXfoy3NqG1ocrNtXzZp9rSG2isN1xurCwTYLY1KjuWXKIEb3j2JgfChpsaFEBJ/gwkyFz0N4Egw7vxvfgRDCX2mt2Va5ja8PGHvGri5djcPlICQghAnJE7h55M2c0f8MBkQMOP7F/JD8pSrA0QyrX4Sh50N07/xB9zdaayob7J0P/3XPWy2taaastpkWZydhNTjAE0jHD4ztNKgmRgYRGii/AsTpq2mpYVXJKmNIccm37Kkxdp5JCElgaupUJiZPZGLyRBJCE0yuVIiuo7Vm9+F61uyt8vTEbj1U69nXOiM+jGlD4hmXHsO4AdEMT4o49SkXVXth28cw7QGwyurEQghDdXM1K0pW8PV+I8iWNZYBMDRmKDdk3cCUlCmMSxyHrQ/83pC/aAVsehcaDkO+bOPTE5wuTWltE/srGymubGR/VSMHq5s8w4DL3HNWvYentYoKsXmC6YSMMK+eVa+wGhEsK1uKbmV32llXts4zL3ZD+QZc2kVIQAjjk8Zz7bBrmZQyyefn5oju0VuHH9c02Vm/r9odYCtZs6+KqgY7YEzRGDsgmjunD2ZcWjRjB8QQGxbYdd989YugFOT8oOuuKYTwO06Xk03lmzwLPH1/+Htc2kVkYCSTUiYxOWUyk/tPJjE00exSe5yEWgGrnoXYwZAx3exKegWnS3OoponiykaKKxuM4FrZSHGVcf9AVeMRgTU6tC2sDkoI8wRU79uEiCCCbRJWRc/TWrOjaodnv9iCQwU0OhqxKiuj4kdx2+jbmJQyiTHxY/rEp8Hi2HrD8GOXS7OjrM4Ir+6e2O2ldZ5FnIYkhjMzqx85aTGMS4shMzEca3ct2uS0G6F2yEwZTSVEH3S48bCnJ/abkm+obq5GoRgVP4p5Y+YxOWUyo+JH9fl1Kfr2uxdQst5Y8fi834JFVqI9EQ6ni4Oe0OoOrO7wWlzVQElVEw5X+9CaGBFEakwIY1KjuWB0MqkxIaTGhJIaE0L/6BAJq8LnlDaUeubFflvyLYcbDwPGAhOXDL6ESSmTGJ803q+W+xfiaKoaWlizr4o1e4we2LV7q6htdgDGCJlxadFcODqFcWnRZA+IJiqkBz+82foh1B2CvB/23PcUwg9UN1dT1lCGRVlQSmFRlrYvjGNWZW17DgsWi/tWdXje65jZ7C47a0vXGkH2wNdsqdgCQFxwHGemnsnklMlMSplETHCMyZX6Fr8Otb11iFOPWvUsBITA2DlmV+Iz7E4XB6vb97S23t9f1UhJdZNnzhQYI8KM0BpKTloMqdlGYO0fHUJqTAgpElqFH6i311N4qJAVB4xVindW7wQgNjiWCUkTmJQyiYnJE0kOTza5UiFOj8PpYuuhWs9CTmv2VrLrcD0AFgXDkiKZPTbF3QsbTUZcmLlb5xQsgKgBkHmueTUI4SPqWur4fN/nfLT7I1YcWIFDO7r0+grVLvBalMVzzPtL0SEwd3iu4zFPqO4Ypr2OW5QFh3aw4fAG6u31BKgAxiaO5d6ce5nSfwpDY4ZiUdIBdTR+HWp7wxAnUzVWwfdvwJirIKTvfNrT4mgNrQ3tg2uV0etaUt2Id0erUpAUGUxqTAjjB8Z6wmprT2tydDBBARJahX9xuIyGs7U3dn3ZehzaQZA1iNx+uVyaeSkTUyZKIyr83uG6Zq/FnCpZX1xNQ4sTgLiwQMalxXBFbqrRC5sa7VurvZfvhF3L4OxHwSLtjOibGh2NfFH8BYt3L+bL4i9pcbWQHJbMDSNvICsuC7SxfZwLl3GrXWitcWkXTu007ns953m+w7F2z6PbXtvJc0c7v/Xa7V7bSV3ex53a6Tlfo7kg4wIm95/MhKQJhAeGm/3P7zd86De36HFrXwZ7A4zvXZ8JNDuclFQ1tetd9Q6vBztsbm/xhNZQJmTEth8aHBNCclQIgQHyR31fo7Vmc8Vmlu5Zyid7PqG4ttjskrpUa2OqUIyIG8GNI29kUsokxiaOJcgqe1WLU2fmKKoWh4vNJTWehZxW761kX0UjAAEWRVZKJFfnDWBcWjTjBsQwIDbEJ4YbHlXh82AJgHE3mF2JED2qxdnC1/u/5qOij1i2bxmNjkbiQ+K5athVzBo4i+yEbN/+f1f0OAm1fZXLZQw9Ts2H5DFmV3NSmuxODlS1rRzsPUR4f2Ujh2qPDK3JUUbv6hmD4+kf09rTGsKAmFCSooKxneo2C6JX0VqzqXwTS/YsYUnREorrirEqK/lJ+Zybfi6K3tWADosdxoSkCUQHR5tdiuhFenIU1cHqJtbsrfRsqfP9/mqa3Xty94sMIicthhsmppOTFsOo/lH+NRXE3gRrXoLhF0JEktnVCNHt7C47K0tW8tHuj/hs72fU2muJDormokEXcX7G+eQk5mCVEQviKCTU9lW7l0HFTpj+oNmVdKq6wc6einqKyhvYW95628CeinoO1TS3O9dqUaREB9M/OoQpQ+KPWIRJQqs4Fq01Gw5vYMmeJSzds5T9dfsJUAFMSJ7AbWNu4+wBZ0voE8IHNNmdbDxQ7ZkLu3pvJSXVTQAEWi2M6h/JDRPTGeeeC5sSHWJyxadp8yJorJAFokSv5nQ5WV26msW7F7N0z1IqmysJt4VzdtrZnJ9xPhOSJ2CzyKr64vgk1PZVq56D0HjIusSUb6+1prS2mT3lDewprzduK9oCbHWjvd35iRFBpMeFMiUzgbTYUE9Pa2psKP0igk59Q3vRJ7m0i+8Pf8+SIiPIltSXEGAJYFLyJG4fcztnp51NVFCU2WUK0eftKK3lP9/uZc2+KjYdqPZsh5YaE0LewFjGDYgmJz2GEckRvW9tg4IFxnZ7A6eZXYkQXUprzfrD61m8ezEfF31MWWMZIQEhTE+dzqyMWUzuP1mmwfghrTW6sRFnTQ3O6hqc1VW43PcjZ52HJSysW7+/hNq+qGqfsUXA5J9AQPf90nA4XRyoaqKovL5dYN1b3sDeigYa7U7PuVaLon90COlxoVw0JpmBcWGkxYWSHhdKWmwooYHyoypOj0u7WF+2no+LPmbpnqUcajiEzWLjjJQzuGvsXUwfMF2CrBA+prLBzqur9jImNZpbpgxyz4WNJjEy2OzSulfpZti7Amb+WrbbE72C1potFVtYXLSYxbsXc6D+AIGWQKb0n8L5GeczLXUaobZQs8sUgKulBVd1tTucVuOsrnaH02ojrNa4A6vnvnGuq7oabbd3es2QsdkEDR7crXVLUuiLCp83bvNuPu1LNdmd7K1oOKLHdU95PfsrG9vt1xoUYCEt1giqU4bEewLrwLgw+seEyBBh0eVc2sWa0jUs3bOUpXuWUtpQis1iY3L/ydybcy/TB0yXfVaF6GJduVBUTloMGx4/r++Nxil4HqxBkD3X7EqEOC27qnbxUdFHLN69mKKaIgJUABNTJnLXuLs4a8BZ0gZ3E+1w4KytNcKpJ6C27z31DqOec2pq0I2Nx7y2JSICa2Qk1qgoLFGRBCUOwRoVhTUqEov7uDXSeGyNisISGYWtX2K3v2cJtX2NoxlWvwhDZ0F02gm9pLrR7gmseysaKDrc2vPawMGapnbnRgQHMDAujFH9o7hwdPse134Rwebu9Sf6hNb5Oa2rFpc1lnk+DZ6ZO5MzU8+UJfKF6EZduVCU1aKgly3Qdlwt9bDuVRh5KYTFmV2NECdtX80+o0e2aDHbKrehUIxPGs8PRv6Ac9POJSa472wjeTq0y4Wrrs4Im1XVuGqqPfedNTXG43a9p9WeHlZXXd0xr61CQz3B1BoZiS09jeAjwmik8TjaOMd4HImy+uZUDwm1fc2mRVBfBuNv9RzSWlNW2+zuYT2yx7Wqof1QgoSIINJjQzkjM46BcWHtelyjQ22yxLrocQ6Xg9WHVrNkzxI+2fMJ5U3lBFmDmNp/KjMHzmRa6jTCbN07l0MIIbrEhreguVoWiBJ+5WD9QT4u+pjFuxezoXwDAGMTxvJg/oPMTJ9JQmiCyRWaz9XUhLOiAkd5Bc5K921FBY6KcpwVle77FW3htLbW2K3kKJTNhiW6NYhGYevXD+vQIVhaw2lkJNboI8OpNTISFRjYg++8Z0io7SMcThcl1U1EfvkPVGgaf9uaRNE3BZ6hw97zWy0KUqJDGBgXxgWjk0mPDSXdK7z61Mb0os9yuBwUHCpgSdESPt37KRVNFQRbg5ma6g6y/WV+jhDCDxUsgIQRMGCC2ZUIcUyHGw+zdM9SFu9ezOrS1QBkxWVxf+79nDfwPFLCU0yusHu5mpvbhdTOA6tx66yowNXQ0Ol1VGAg1thYrLExBMTEEpiW5hna2xpYrVFGGLVERXl6V1VwsHQkefHrdGLmBu++rLrBTuHeCgqKKtl4oIY95fUUVzYyRBfxUVAh/22/nn+v2GvMb40NZdLguLZhwrGhpMaEEhjQx+YvCb9gd9lZdXAVS4qW8Nnez6hsriQkIIQzU89kRvoMpvSfIkFWCOG/DqyBA6vh/N+D/LEqfFB1czWf7v2Uj3Z/xMqDK3FpF4OjBnP32LuZlTGL9Mh0s0s8Za6WFk8APSKcVlbg7BBUXfX1nV/IZiMgNhZrbCwBsbEEpqcTEBuDNTbOCK5xcVhj3LexsVjCwiScdgG/DrU9ucG7r9Jas7eigYKiSgr2VFK4p4Jth4xx9AEWxZB+EWSlRHL+6GSu2P8mzgPB3HrHozySmCTzW4VfsLvsfFfyHUv3LOXTvZ9S3VxNaEAoZw44k5npM5ncfzIhAX6+H6UQQoCxQJQtFLKvMbsSITzq7fV8tvczFhct5psD3+BwORgQMYBbRt3C+RnnMyRmiNkldkq3tOCorDyix/RovalHnYcaEEBATAzWuDgCYmMIGTDACKetITU2Fmus8Zw1Lg5LeLiEVBP4dajti+xOFxsP1FBQVEHhHiPIltU2A8YiTTlpMVycnUJueixjB0QTEuiezN1YBX/4EMZcRXJSsnlvQIgTYHfaWVGygiVFS/h83+fUtNQQZgtj+oDpzEyfyRkpZxAc0Mu39BDCT8koqlPUVA3fL4RRV0CwbC8mzNXoaOTL4i9ZXLSY5cXLaXY2kxSWxPUjrmdWxiyyYrN6PLi1Dvd1Vla6w2qlEU4rjF7Ujr2prtrazi9ktXqG+lrjYgkZNcoTWI/oTY2NxRIZKSHVD0io9XHVjXZW762ksKiSgj0VrN1XRZPdmDSeGhPClMx4ctNjyBsYw9DEiKP3vq57BewN7RaIEsKXtDhbWHFgBUv2LOHzvZ9Ta68l3BbOWQPOYubAmUxKmSSbsQvhB2QU1Sla/zrY62WBKGGaFmcL3xz4ho92f8SyfctocDQQFxzH5UMu5/yM88lOyMaiumZ6mna5jJV7K6uM+ajtgqo7rHo9dlRWoo8yJxWLxRjq6+5NDR6ZhdUdWL2HAbf2ploiI1Gy/3OvI6HWh2itKa5spGCPMR+2oKiSbaW1aG1sazAyJZI5+WnkpceSNzCGfie6+bzWsOpZSB0PKWO79T0IcTKanc18s/8bluxZwrJ9y6iz1xERGMFZaWdx3sDzmJg8kUBr71uhTwgh2tHaGHqcPBb655hdjehDHC4HK0tWsrhoMZ/s/YTallqigqI4P+N8zs84n7x+eVgtx9/CxdXUZITPigp3UD0ymDorKnBUuR9XVR11ZV8VEmIE1JgYrLGxBA7KMHpVY2KwxhrHA9zPWWNisEZFSUgVEmrNZHe62FxS454PawTZ0tahxEEBjEuP4cIxyeSlx5A9IPrUVx3etQzKd8Bl87uueCFOUZOjia/3f82SPUv4ovgL6u31RAZGMiN9BjPSZzAxeSI2q83sMoUQoufsWwmlG2H2X8yuRPQBLu1i9aHVLC5azNI9S6loqiDMFsbZA85mVsYsJvabgKWuAWdlJc1r1rqDaqURVisqcFZV4mgXVKuO3YsaHW0M6Y2OIShjENbcWKwx0UbvaUyM0ava+jg6GkuIrJMhTp6E2h5U02Rnzd4qCoqMALt2X5VnK53+0SFMGhxHXnoMuemxDEuKcG863wVWPQuhccZG7kKYoNHRyFf7v2Jp0VK+KP6CBkcD0UHRzBo4ixnpM8hPzsdmkSArhOijChZAUKQxn1b4La21cYv2PPbcR+O+2+5YZ6/xXO9o1+l43lGed2kXNDWjGxrQDY0cKt1N4dbP2bJ7JbqqhthGK/dakxnoHElMcwCuyvU4K79g57F6UUNDvXpRYwgaPAhrdGuv6ZFBVXpRRU+RUNtNWocSG4s5GSF26yFjKLFFQVZKJNeMH0DewBjy0mNJiuqmRW+qi2HrhzD5XgiQ+Yii5zTYG/hy/5csKVrCl/u/pNHRSExQDBcMuoCZ6TPJS8qTICuEEA0VsPFtyPkBBIWbXU2f43Q52V+3n51VO9lZvdO4rdrJnpo92F32Ew6jXUFpTVALBNshpBmCPfc1Ie77wS3u5+zuY83uc1oguFl77oe0GOdavMqzARPdXwBYwBpTjzXGhiUmFtvgwV7De1vDaYyxgJL7viVYFmkUvklCbRdxOF1sLqk1AuweY2GngzVNAIQHBTAuLZrzRyWTNzCGsaczlPhkFb5gzNWRhSd8Sl1LHTurd7Kjcoen4exNDtYf5Kv9X9HkbCI2OJbZg2Yzc+BMcvvlEmCRXzvi6LTTiW5qwtXU5Ln1vt/utrEJ3dz+Nv6O2wmIjzf7bQhx4ta9As5myLvZ7Ep6tdbwuqNqR7sAu7t6N83OZs95/UL7MTh6MDn9cgi2GgFOKYWi/eg5pRTKpbE22bE1ObA22QloshPQaMfaZMfa1EJAk3HfONZiHG+0E9DUgrXdsRasTS1Ymh2oE8zIziAbruBAnCGBuEICcYUG4owLxBUShCs4kPqQIGpD2h67QoMJDI9k5JAziExMJSBGFkwSvYv8dXmKaluHErv3hl2zt4qGFmMocUpUMPkZxmJOuekxDE+K7LqhxCfD0QKF/4KhsyA6ree/v6DB3sDu6t1sr9rOzqqd7KjawY6qHRysP+g5J9AS2OtW9Q0LDOOSzEs4b+B55CTmnNAiE8J3aZfLCJHNzejGRlxNzeimzm9dTY3oxiZczU3HuXVfy3NNI6hq+6l9wKNsNlRICDHXXiOhVvgPrY2hxwMmQr+RZlfTKzhdTorritlRtYNdVbuM2+pdR4TXpLAkBkcNJj8pn8HRgxkcPZiMsDQCi8to2rSJ5tXbcNaU4apvwNVQb9zW1+NqcN/W16Obmk64LktoKJawMOMrNBRLeDyWfmFHHm+9H+a+H+p1v/Wc0FCUVdpVIbxJqD1B+6saPXvDriqqZOvBGlzuocQjkiO5KjeV3IGx5KXHkBLtIxPcN70D9aWQL9v4dLdmZ7MRXiu3e4Yuba/azoG6A56hSYGWQDKiMsjtl0tmdCaDowaTGZNJ//D+XbZEvhDetNNpbJlQXm4s6lFR7tnPz1lVafxh1hpEO71tQjc2oltaTq2AgAAswcGokGAsQcFYQoJRQcFYgoOxRkZi6ZdoPG497n0bbJx3tNuOx+QPPN8h+9SehKIvjYUcp/0/syvxO06Xk321+zw9rq0hdnf1blpcbb+zksKSGBw9mAlJEzzhdVDUIEIJpHn7dpo2baJp0/c0bXqNA1u2opuN4KuCgozhtl5B05aS0hY4vcNmJ8HUGhaGCg01bkNCpEdUiG4mobYTDqeLLQdr3QHWCLIl1cancWGBVsalxfDjs4eQNzCGcWkxhPfUUOKTUb0fFj8EiSNh0NlmV9Nr2J12imqKPD2urQ3pvtp9xoIMQIAKYGDUQEbHj+bSzEvJjM4kMzqT1IhUGXorTkvbvn6VnQZVz8bzlRU4yiuOuWWCNSoKi/uPLUtQECokBGt4BCo+wR0Ug7AEhxz19rhBNCgIZZM5032R7FN7EgoWQEgMZF1idiU+y+FysK92n6fXtTXEFlUXtQuvyWHJDI4ezMTkie3Ca3hgOK7GRpq3bqVx1SaaNr1O6aZNNG/fAe6RIZaICIKzsoiZO5fgrCyCR2YRmJ4uH5YJ4UfkL2y3fRUNLCwspnBPJWv2VlLvHkqcHBVMbnoMeekx5A2MZXhSBAFWH/+0zdECb9wIjia46gWQTwdPmsPlYG/tXnZU7mg3bHhvzV4c2gGAVVkZEDGAoTFDOT/jfAZHD2ZI9BDSItNkASRxQrTWuGprcZSXG/v7lZcb+/dVlBuhtKLCCKytt5WV4HR2ei1LVJRn4/nAgQMJyck1tlCIjTNu4+KwxsQSEGdsmaAC5Ne/EKaqK4XN78GEO8Ami++0htfW0U6t8153V+9ut+5ESlgKg6MHc0bKGUZ4jRrMoOhBhNnCAHDW1dG8eTNN36yhZtNLRoDducvzAZ81JobgkSMJv3mqJ8DaUlNRyoRpYkKILiN/1bgdqmniL59tZ3hSJJfnpBqrEg+Mpb+vDCU+GUsegeJVcNW/IGGo2dX4tNaFIzrOeS2qLvI0ogpFakQqmdGZnJN2DoOjB5MZnUlGVAaB1kCT34HwJVprXPX17l5UrzDq3XvqHVQrKz09BR1ZIiI8odQ2YAAh2dlYY92htDWcxhpfATEx0isqhL9Z829wOSD3JrMr6VGtHxp7el7d4dW73QXoH96fQVGDmJwymUHRgzztbmt4BXBUVhoBdtMrVG3aRNPGTbTs2eN5PiAxkeCsLCJmnkfwyCyCs7II6NdPAqwQvZDPhVql1CDgESBKa31lT33fsQOiWf/YTCKC/fwPw/VvwMr5MOlu2ZfWi0u7KKkvYUdl+2HDu6p3tVs4IiUshcyYTKb0n+IZNpwRlUFIgB9+uCG6hKulBWdZGY7Dh48aVB0V7h7W8vKjLnRkCQszAmhsLLbkZIJHjTS2TIgzjllj44xtE+LijHlcgfKBiRC9lstp7E6QMQ3ih5hdTbewu+zsqzHmvHov2lRUU4TD5fCc1z+8P4OjBzOl/xRjrQl3uxtqC213PUdZGU3fFHB40yZjHuzGTdgPHPA8b+vfn+CsLKIuu9TogR0xgoCEhB57v0IIc/VIqFVKLQAuAkq11qO8js8C/gxYgWe11k9orXcBtyilFvZEba0CrBYifH1Y8fGUbob37oG0SXDu42ZXYwqtNYcaDrULrjsqjTk4jY5Gz3mJoYkMiR7C+KTxnvDqPXxJ9H7abjeCamkpjrIy7KWlxv3SMvetcdxZWdnp61VIiDuMxmJLSCR42HBPL2q7oOruUbUE9a4VroUQp2HnZ1C1F2b8yuxKTpvWmuLaYjZXbG63z2tn4TUzOpNpqdPaVhuOPDK8aq1xlJRQu+kbT3ht2rQJR1mZ55zAgQMJGTuWmOvmegKsNTq6p96yEMIH9VRP7QvA34AXWw8opazA34EZQDGwSim1SGu9qYdq6l2aauC16yEw3JhHa/XzHufj0FpT3lTuWW3Ye+GmOnud57y44DgyYzK5fMjlnjmvg6IHERkYaWL1ojtph8PoPS0txVFW6gmodndIbQ2tzooKYzsNb1YrAfHxBCQmGkN+c3OwJSYSkJCANT6+bV5qrLEiphBCnJKCBRCWCMMuNLuSk9YaYlcdWsWqg6tYeXAlpQ2lgDFdxzu8tn5g3Fl4BWPxu5Y9e9wrELcFWGdVlXGCxULQ4MGEnXGGZ/hw0PDhWMPDe/AdCyH8QY+EWq31cqXUwA6H84Ed7p5ZlFKvApcAEmpPltbw7l1QsRtufA8iksyuqMvVtNRQcLCAVQdXsal8Ezurd1LdXO15PjoomszoTC4cdKGn5zUzOpPo4GjzihZdSjudOCsq2npUvQKqJ7iWleIsrzhyxV+LhYC4OCOs9utHyOjRBCQmEpCYYBxrDa6xsbLapRCie1UXw7bFMOU+CPCPaQbFtcWsOmiE2FWHVnn2Wo8NjiU/KZ/xSeMZFT/qmNN1tNNJy+7d7cJr0+bNuOrcH0TbbAQPGULEjHON3tesLIKGDsUSItN/hBDHZ+ac2v7APq/HxcAEpVQc8BtgnFLqIa31bzt7sVJqHjAPIC0trbtr9W0r/g6bF8GM/4aBk82upks02BtYXbqalSUrWXlwJZsrNuPSLoKtwQyPHc6M9Bme4Do4ejBxwXGy8IOf0i4Xzqqq9uHUE1wPtx0/fLjTlX+t7rAakJhA8MgsAhLawqpxP5GAuFhZ7VcI4Ru++qPxYXTOjWZXclQH6g6w8uBKVh1cRcHBAg7UG3NXY4NjyeuXxy2jbiE/KZ+MqIxO215tt9O8c2dbeN20iaYtW9CNxjQgFRxM8LBhRF08uy3AZmaiZC0BIcQp8rm/8rTW5cAdJ3DefGA+QF5enj7O6b1X0dew9BcwYjac8WOzqzllzc5m1pWuY+VBI8R+X/Y9Du0gwBJAdkI2t4+5nfykfMYkjJEVh/2E1todVsu8elZLPcOC7a3zV8vKwOE44vXWmBgCEoxwGjRkyBG9qgGJiQTEx8uqv0II/7HyGVj1rLGNT0y62dWgXS5wOimp3k9hySrWHCxgXclqDtWWYNEQExDJxLjRZKfOZnRMFqmhKSiXC213ondX0+hajXY40Q479n3FngDbvHWrZ9E8S1gYwSNGEHP1VZ4AG5iRIR80CiG6lJm/UfYDA7wep7qPiRNVexAW3gwxA+GS/wM/6qm0u+xsPLzRCLElK1lTuoYWVwsWZWFk3EhuHHkj+cn5jEscJysPm8jV3IyrpgZnbR2uulqctbW4Wu/X1LqP1eGqrcVZZzznrK3BVVWNo6ys05WALVFR2BITCEhIJCg/w93LmujpbbUlJmJNSJDVf4UQNG3bRsW//mU8cLdxbT2Dqt3xdvcVJ/aaI245vfM7vs77WMUu2PQOxI6FyH6w8Sm004V2OsDhRLuc4HSiHU5wGbfa6T7mdBrntTvfBQ6H8ZzLaRzr9HyncZ73+e4v5bWuQKb76yrv/wBUAsuB5TiAIo7NEhVFcNYIYm/8gSfA2tLSUBY/X4hTCOHzzAy1q4AhSqkMjDB7LTD3ZC6glJoNzM7MzOyG8nyc0w5v3AzNtXDDOxDs2wsfOV1OtlZuZWXJSr47+B2rD62mwdEAwLCYYVwz/BomJE0gp18OEYERJlfbO7haWnDV1hqBs10o7RBOW0Np63PucOqqrT3q9jQeSmEJCzP2VI2IwBIRgS0hEUtmptGj6t2r6r5vCQ7umX8AIYQpurJtdlZVUf/1N22LunW41XgN1Gq9e5Rzj3xt957f7jUul7EnLeFgKYcVC4z5+wEBKIul3X0CAozHVgvK2nrf6j7HirIY95XNhgoKMo5ZA455fqNu4VDzYQ41lXKg8RBVjhpcCgJsISRH9iclagADotOJD++HxRqACvC6htXqvq4FWp+zWFCtdVqsKKuFgORkbP37y1QgIYQplO64+md3fBOlXgGmA/HAIeAxrfVzSqkLgD9hbOmzQGv9m1O5fl5eni4oKOiiav3Ex4/Air/B5c/AmKvNruYIWmt2Vu3ku4PfsbJkJQWHCqhpqQEgIyqD/KR8JiRPIK9fHjHBMSZX63u03Y6zrkPYbBdEvXpFO4bTujpcNTXolpbjfh9LaKgRSCMjsIRHYIkIxxoeceSxiAgs4e7biEisEeFYIiKwhIXJJ/CiV1JKFWqt88yuw5/1yba5M5VF8Oy5YAuBWz+F8MRu/5ZlDWWeRZ1WHVzFnpo9AEQERpDbL9ezuNPQmKFYlPwOF0L4h2O1zT21+vGcoxz/EPiwJ2roVTa+YwTa8bf5TKDVWrOvdh/fHfyOVSWr+O7gd1Q0VQDG3nTnpp/L+KTx5Cflkxja/Q26r3M1NtK4/nsaCgtoWv89jsqKdiFVNzUd9xoqJMTTO2oND8caFYUttb8RSiPdPafhEW0BNDwca2Rk27HwcFnpVwghulNDBbx0lTG66qYPuy3QHm48TMHBAs/iTkU1RQBE2IwQe9XQq8hPymdozFCsFvm9L4Toffx6ln6fHH58eDu8ezf0z4Pz/sfUUg7WH2TlwZV8V/IdKw+u9CzxnxCSwKSUSUxImsD4pPGkRqSaWqcvcFRW0rhmDQ2FhTQWFNK4aRPY7aAUgYMHYUvshy0p2d0rGtm+xzQi3AiikRFt4TQ8XBZIEkIIX+ZoNvaPrywypgklDO2yS5c3lrPq0CrPVne7qncBEGYLI7dfLlcMuYLxyeMZHjNcQqwQok/w61CrtX4PeC8vL+82s2vpEc11RgMZEAhX/6vH97dr/SS4dUjx3tq9AMQExZCXlMeto24lPzmfgZED+/ycGvv+/TQUFtJQuJqGwgJaduw0nrDZCBk1iribbiQkN5fQceOwRkWZW6wQQoiu5XLBOz+CPV/DFc+d9nZ7FU0VngC76uAqdlYbbUpoQCg5/XK4NPNSxieNZ3jscAIsfv2nnRBCnBL5zecvtIb37oXD2+D6tyCq+3s/q5urKThU4NkrdkfVDgDCbeHk9cvj2uHXkp+Uz5CYIX16To52uWjevoPG1YU0FBTSsHo1jpISACzh4YSMG0fURbMJzc0hePRoWShJCCF6u89+BRvehHMfh9FXnvTLK5sqKTxU6BlO3Nr+hgSEkNMvh9mDZ5OflM+IuBESYoUQAgm1/mPlM7BhIZz9cxh8Vrd8iwZ7g6cRXXlwJZvLN6PRBFuDyemXw0WDLmJC8oQ+/0mwq6WFpg0bPEOJG9aswVVjLIIVkJBASF4uobfcQmhuDkFDh8q8VSGE6EsKFsBXf4Tcm2HyT07oJdXN1UZPrHthp22V2wAjxI5LHMeFgy5kfNJ4suKysFlk6okQQnTk18mkz8yp3bcSPn4Yhs6CKfd32WWbHE2sK1vn2St2w+ENOLQDm8VGdkI2Pxr7I/KT8hkTPwabte82os7aWhrXrjV6Yd0LO7WuLByYkUHkeTMJycklNC8XW2pqnx96LYQQfda2JfDBf8GQmXDBU0fdP766uZrCQ4We4cTbKrd5PkQelziOe8bdw/ik8YyMHykhVgghToBfh9o+Mae2/jC8fiNE9YfLnobT2D7F7rKz8fBGz8JOa0vX0uJqwaqsjIwfyU2jbiI/KZ+xiWMJCQjpwjfhX+yHStsNJW7eutWYH2W1EpyVRczcuYTk5hCam0tAbKzZ5QohhPAFB9bCGzdB0mi48nmwHvknVoO9gd+u/C3v7ngXjSbIGsTYxLHcNfYu8pPzGRU3qk9/iCyEEKfKr0Ntr+dywsIfQmMF3LIUQk5uP1eXdrG5YjMrS1by3cHvWH1oNY2ORhSKYbHDuHb4tUxInkBOYg7hgeHd9CZ8m9aalt2724YSr16Nfd8+wNgyJ2RsNvF33klobg4h2dlYQkNNrlgIIYTPqdoLL18NobEw93UIOrJN3V29m/uX3c/Oqp1cN+I6zk0/l9Hxowm09uyij0II0RtJqPVln/8Gdn8Bl/wdksec1Esb7A389Iuf8uX+LwEYFDWISwZfwoTkCeT1yyM6OLobCvZ92m6nacsWGgoKjd7YwtU4K4z9dK2xsYTm5hAzdy6hebkEDx8u2+YIIYQ4tsYqYy9aexP8YBFEJB1xyuKixTz29WMEWYN4esbTnJFyRs/XKYQQvZhfh9pePad260fw5f9Czg9g3PUn9dKKpgru+uQuNlVs4r9y/4sLB11IQmhCNxXq21z19TSuX+8eSlxI47r16IYGAGwDBhA+bZp7KHEegRmyFZEQQoiT0LoXbflOuOEtSBze7mm7085TBU/x8paXyU7I5qkznyIp7MjQK4QQ4vT4dajttXNqK3bBW7dDcjac//uTeun+uv3cvvR2DtYf5E/T/8RZad2zUrKvclRUtBtK3LRpEzidoBRBw4cTffnlxlDinFxs/RLNLlcIIYS/0hoW/RiKvoTL/397dx9nY53/cfz1MYybcRu7FSIxJSmVIW2xtihqtZFEpRWlSbcidMO23Wlt2cqymkpy1xJWRAqFVq3FEm1FbtZNsUPDlJsxZub7++PMr50wnJm5zrnOdeb9fDzO4zHnOte5zvvjOud8fc519yo0aPOTh3fu38nAJQNZu2ctPZv0pH/z/jrpk4hIhAS6qY1LRw7B1NtCZ0zsNgHKhX9N0/UZ60ldmEp2bjavXvUqF/38oggG9Z9zjiM7dvxvV+KVq8jesgUAS0yk4gUXUPOOO6iU0pyKF15IQpUqPicWEZG48dEzsHYqXPE4XNDtJw/9/Zu/M+TjIeTk5TCy7Uja12/vU0gRkdJBTW2+3H37yPrqK39DOAfLXoZNG+CKYbB+F7ArrKeuz1jPqNWjOLdsRfo3H0DtzYc5sPkfkc3rh7w8Dm/ewsFVKzm06l/kpKcDUKZqVSpdfDHVunSmUvMUKjQ9jzKJOvmGiIhEwKo3YekfQ4cItR744+TcvFz+8tlfSFubRnKNZEa2HUn9qvV9DCoiUjoEuqn18pjarK++Yluv20seyhM14aPRRXpGRWAQAAfIeeMxtkUgVSwpe9ppVGrRIrQV9uLmlE9uhJXgckciInJiZpYELAGecM6963ce32xcCO/2h4ZXwrUjf7wW7XeHvmPwx4NZvnM51ze6nkcvebRUXx5PRCSaAt3UenlMbYUmTag/cYIHqYppzwaYOzB0luP2T4KF16At3LaIyV9MIrlGMvdf/ACVyyVFOKj/ytWuTbk6dfyOISISCGY2Dvg1kO6ca1pgegfgJSABeM0599xJFjUYmBaxoEGwc23o2vGnNoFub0L+NWVXp69m4OKBZGZn8uQvnqRzcmefg4qIlC6Bbmq9lFC1KpVatPDnxQ9mwCu9oWFN6DsFkmqe9CnOOUatHsWr+yfzq8uu4PE2I6hQNvzjb0VEpNQYD/wZ+PGXWzNLAEYD7YEdwAozm02owR1+1PN7A82AL4DSO9Bk7ghdi7ZCNbj5bShfBeccE76YwJ9W/YnalWszqd0kGp/S+OTLEhERT6mp9VteHsy8E/bvgt7zw2poc/JyePLTJ/nbxr/R9eyuPHbJY5Qto1UpIiLHcs4tNbMzj5rcEtjonNsMYGZ/BX7jnBtOaKvuT5hZWyAJaAIcMrN5zrm8SOaOKVmZoWvRZh8IjdVVT+eH7B8Yumwoi7Ytol29djx52ZNUSdQJCUVE/KBOyG9L/xg6PufakVCn+UlnP5RziIeXPMySHUtIbZZKv2b9dG1VEREpqjrA9gL3dwCXFDazc+4xADPrBewprKE1s75AX4B69ep5ldVfOdkwtWfoMKFbZ8Cp5/FVxlc8tPghdu7fycMpD9OzSU+NxSIiPlJT66evF8Li4dCsB6T0Puns+7L2ce+H97J291qGthpKt3O6nfQ5IiIiXnHOjT/J42lAGkBKSoqLRqaIcg7mPABblsD1f4Gz2jLz65k8849nqF6hOuM6jIv7y+eJiARBoJtaL89+HHX7tsHMO+DnTX5y9sTC7Ny/k7sW3sU3P3zDyLYjaVe/XZSCiohIHPoGOKPA/br506Sgxc/BZ1Og7SMcatqZZ/7+OO9seodWp7fiudbPUbPiyQ8ZEhGRyAv0NVCcc3Occ32rVavmd5SiyTkM026DvFy4aSIkVjrh7Bv2buDWebey5+AeXmn/ihpaEREpqRVAspk1MLNEoDsw24sFm1knM0vLzMz0YnH+WT0ZljwHF97Cfy66iVvm3cLsTbNJbZbK2HZj1dCKiMSQQDe1gfXeYPh2NXQeCzUbnnDWlbtW0uu9XgCM7zielNNSohBQRETihZm9BXwKnGNmO8ysj3MuB7gXeB/4EpjmnPu3F68X2B+cC9r0Ecy5H85qywcXdKL73B7sPribMe3GcM+F95BQJsHvhCIiUkCgdz8OpDVTYNUbcNmD0PjaE866aOsiBi0dRJ0qdRjbbiy1K9eOTkYREYkbzrkehUyfB8yLcpzYt+tzmNqTI7XOZmRyCpM+HsQFP7uAF375AqclneZ3OhEROQ41tdG0ax282x/ObA1XDD3hrNPWT+OZ5c/QtFZTRl8xmuoVqkcno4iISAkE+nwX338Lk29kV8UqDDzjDD7bMI1bz72Vh5o/RLmEcn6nExGRQmj342g5tC90SYCKNaDrOEg4/u8JzjnGrBnDU/94itZ1WvPaVa+poRURkcAI7O7HWd/D5G4s4xA3nlqdr7/fyvO/fJ7BLQeroRURiXHaUhsNeXkw627I3A695kLlnx93tty8XJ5e/jTTN0ync6PODLt0GGXLaBWJiIhEVO4RcqfdxivZ2xlbqxoNk05jZNuRNKjWwO9kIiISBnVM0bDsRVg/Dzo8B/VaHXeWrJwsBi8dzIfbP+TO8+/kvovu04XcRUREIs05Mmb3Y8jBz/m0elWua9iJx1s9TsWyFf1OJiIiYQp0UxuI43Y2L4EPn4LzusAlqcedJfNwJvd/eD+r01fzSMtHuPncm6McUkRExBuBGJsLWPPBIAZk/J19FZN44tJhdEnuoh+VRUQCJtDH1Mb8cTvffwvTe0PNZLhuFBxnkNx1YBe95vdi3Z51jPjlCDW0IiISaDE/NudzzjFxwUPcvvM9EsslMenat7jh7BvU0IqIBFCgt9TGtJxsmPZbyMmCmyZB+crHzLJp3yZSF6byQ/YPjG03lpant/QhqIiISOmyP3s/wz64mwXfreFXVOLpru9RtVJNv2OJiEgxqamNlAVDYcc/oesb8LOzj3l4Tfoa7ll0D4kJiYzvMJ7GpzT2IaSIiEjpsj5jPQMW3ceOA98yILs8v731faxSDb9jiYhICQR69+OYtW46LB8LrfpB0y7HPPzRto+444M7qFGhBhM7TlRDKyIiccPMOplZWmZmpt9RjjFr4yxumXszB/d/y+uZOfTqNksNrYhIHFBT67X0L2H2fXBGK2j/5DEPz9gwgwcXP0hy9WQmdJxA3Sp1fQgpIiISGbF4TG1WThbDlg1j6LKhXHgkj2m79tK821Sofobf0URExAPa/dhLh3+AqT0hMQluHA8FLtbunCNtbRp/XvNnLqtzGSN/OZJK5Sr5l1VERKQU2Pb9Nh5a/BDr966nr51Cv+3rSLh5KpzezO9oIiLiETW1XnEO3rkHMjbDbe9A1dN/fCg3L5fh/xzO1PVTua7hdTzxiycoV6bcCRYmIiIiJbVw60KGLhtKgiUwpvIFtF73Lvz6RUhu73c0ERHxkJpar/xjDHzxTmiX4watf5x8OPcwj3z8CAu2LuD2prfT/+L+ulyAiIhIBB3JO8KLq15kwhcTOL/W+Txf8RxqLx0Jlz8EKbf7HU9ERDymptYLWz+BD4ZC41/DL+7/cfL32d/zwIcPsPK/KxnUYhA9m/T0MaSIiEjkmVknoFOjRo18ef3/HvgvDy99mNXpq+nRuAcDyzcgcVYqNO0KVwz1JZOIiESWmtqS+uG/8PbtUONMuH4M5G+FTT+YTurCVLZkbmFEmxF0bNDR35wiIiJR4JybA8xJSUm5M9qv/em3nzLk4yFk5WTxxzZ/pINVhomdof7loTG6jM6PKSISjwLd1Pr9azC5OTD9dsjKhJ4zoULoTI9bMreQuiCVfYf3MebKMVxa+1J/8omIiJQCeS6PtLVpjFkzhobVG/JC2xc4K/sIvN4eajSA7pOgbHm/Y4qISIQE+idL3y8bsOj3sHUZdHoJTj0PgLW713Lbe7eRlZvFGx3eUEMrIiISQXuz9tJvYT9GrxnNtWddy+RrJnNWmSSY1BUSysMtb0NFXYtWRCSeBXpLra++mA2fvAwpfaDZTQAs3bGUgUsGUqtiLV5p9wpnVNX170RERCLls92fMXDJQL479B3DLh1G1+Su2JGDMKUbHNwDveZCjfp+xxQRkQhTU1scezbCrH5Qpzl0GA7ArI2zeOKTJzi7xtmMaTeGWhVr+RxSREQkPjnnmPLVFJ5f+TynVjqViddM5Lya5+UfFtQbdq2F7m9BnYv9jioiIlGgpraosg/AtJ6QUA5ufBOXkMjr617jpX+9RKvTW/Hir14kqVyS3ylFRER8EenzXezP3s/vPvkdH2z9gLZ12/L05U9TrXy10PXi3xsEG+bDNc/DOR0i8voiIhJ7An1MbdQ5B3MegPQvoevr5FWrwx9W/IGX/vUSHRt0ZMyVY9TQiohIqRbJ811s2LuBHnN7sGjbIvo3789LV7wUamghdEjQytdDl9ZrGfUTL4uIiI+0pbYoVrwG696GXz1O9pmX8+jSQbz/n/fp2aQnA1MGUsb0G4GIiEgkzN40m6c+fYrKiZV59apXaXFai/89+PkMWDAMzusM7X7vX0gREfGFmtpwbV8B8x+B5KvZf8ldPLiwH8t3LWdA8wH0atrL73QiIiJx6XDuYYYvH86Mr2fQ4rQWjGgz4qfnrdj6CfwtFepdCteP1bVoRURKITW14TiwB97+LVQ9nT0dn+XuBX3YuHcjz17+LJ0advI7nYiISFza/v12BiwZwJcZX3Ln+XfS78J+lC1T4L8ue76Gt3pA9XrQfQqUq+BfWBER8Y2a2pPJy4UZfeDAHrZ2n8BdH91LRlYGo64cxeV1Lvc7nYiISNzKOJxB+sF0Rl85mjZ12/z0wf3pMOkGKFMWbpkOlU7xJ6SIiPhOTe3JfPQsbF7M5+0epd+qZwEYd/U4mtZq6nMwERGR+NbsZ82Yf8N8KpQ9agts9gGYclOose31LpzSwJ+AIiISE3TgyYmsnw8fP8+yptfSe+sMKpWrxMRrJqqhFRERiZJjGtq8XJhxB3y7Grq+DnVT/AkmIiIxQ01tYTK2wN/6MqdOY+49+AX1qtRj0jWTqF+1vt/JREREYpaZdTKztMzMTO8X7hzMHwLr50HHP0Dja71/DRERCRw1tcdz5BBM68mbSYk8mniQi09tzhsd3vjp2RZFRETkGJG8Ti2fjoZ/pkGre+CSu7xfvoiIBFLMHVNrZknAGCAbWOycmxztDHlzB/BC9nYmVKvK1WdezbOXP0tiQmK0Y4iIiMj/+/cs+OBxOPc6uOppv9OIiEgMicqWWjMbZ2bpZvb5UdM7mNl6M9toZkPyJ3cBpjvn7gSui0a+go6sGMej38xnQrWq3Nz4Zka0GaGGVkRExE/blsPMvlC3BXRJ07VoRUTkJ6I1KowHOhScYGYJwGigI9AE6GFmTYC6wPb82XKjlA+AA9s+4d5VzzG3chIPXHgfQ1oOoYxp4BQREfHNd5vgre5QrQ70+CuUq+h3IhERiTFR6dicc0uBjKMmtwQ2Ouc2O+eygb8CvwF2EGpso5YP4LuMTfRecBfLKyTyVMpg7mjWFzOL1suLiIjI0Q7sCV2L1ix0Ldqkmn4nEhGRGOTnZsg6/G+LLISa2TrATOAGM/sLMKewJ5tZXzNbaWYrd+/eXeIwW7Z+xPYyjpcv7M/1591a4uWJiIhICe1aC4cyQltoazb0O42IiMSomDtRlHPuAHB7GPOlAWkAKSkprqSvm3LRHbzf8BqqVK1d0kWJiIiIFxpeAQ+ugwoROJOyiIjEDT+31H4DnFHgft38ab5RQysiIhJj1NCKiMhJ+NnUrgCSzayBmSUC3YHZRVlARC/wLiIiIiIiIjEvWpf0eQv4FDjHzHaYWR/nXA5wL/A+8CUwzTn376IsN6IXeBcREREREZGYF5Vjap1zPQqZPg+YF40MIiIiEnlm1gno1KhRI7+jiIhIKRHoi7Bq92MREZHYor2oREQk2gLd1GrgFBERERERKd0C3dSKiIiIiIhI6aamVkRERERERAIr0E2tjqkVEREREREp3QLd1OqYWhERERERkdIt0E2tiIiIiIiIlG5qakVERERERCSwAt3U6phaERERERGR0s2cc35nKDEz2w1szb9bDShul1sL2HOc5YTzd2GvfaI84WQtST1w/JqKUkMs11Tc9VKUmsLNGs2awqmnONOLOs+J6H0XfzWFm9/Pz5KX3+H1nXM/K2YOQWPzSeg7UmNzcbLqfXcsjc0am0Occ3F1A9JK8NyVx1tOOH8X9tonyhNO1pLUU1hNRakhlmsq7nopSk3hZo1mTeHUU9Ka9L7z930XizWFm9/Pz1IkvsN18+YWa+s1Vj57BZdVmr5PSlpTNL5PilpTuN8nJalJ7zt/33exWFO4+f38LEXiO/x4t0DvflyIORFYTjh/F/baJ8oTTlav6im4rKLUEMs1FXe9FKWmcLPG2vuuONOLOk+49L4r/PWDVFO4+ePxsyQlF2vrNVY+ewWXVZq+T8K5X9i0E00v7nxFWY7G5pNPC+exoswTDo3NRZ92ounFna8oy/F0bI6L3Y+9YmYrnXMpfufwkmoKhnirKd7qAdUUBPFWj4TE43pVTcEQbzXFWz2gmoIgWvXE45bakkjzO0AEqKZgiLea4q0eUE1BEG/1SEg8rlfVFAzxVlO81QOqKQiiUo+21IqIiIiIiEhgaUutiIiIiIiIBJaaWhEREREREQksNbUiIiIiIiISWGpqw2RmZ5nZ62Y23e8sXjGz683sVTObamZX+Z3HC2Z2rpmNNbPpZna333m8YGZJZrbSzH7tdxYvmFlbM/s4fz219TuPF8ysjJk9Y2ajzOy3fucpKTNrnb9+XjOzT/zO4wUzq2dms8xsnJkN8TuPeENjczBobI59Gptjn8bmEysVTW3+P1S6mX1+1PQOZrbezDae7B/SObfZOdcnsknD51FNs5xzdwKpwE2RzBsOj2r60jmXCnQDLotk3pPxop58g4FpkUlZNB7V5ID9QAVgR6Syhsujmn4D1AWO4HNNHn2OPs7/HL0LvBnJvOHwaB2dD0x3zvUGLopYWAmbxubj09gcWRqbC6WxOYI0NhfKs7G5VJz92MzaEPqgTnDONc2flgBsANoTeqOvAHoACcDwoxbR2zmXnv+86c65rtHKXhiPa3oBmOyc+1eU4h+XVzWZ2XXA3cBE59yUaOU/mhf1AM2AmoQGmT3OuXejk/74PKppj3Muz8xOBUY6526JVv7j8aim3sBe59wrfn9HePzdMA3o45z7IUrxj8ujdZQLTCf0H7eJzrk3opNeCqOxWWOzHzQ2a2z2g8bmyI/NZYv7xCBxzi01szOPmtwS2Oic2wxgZn8FfuOcGw7E/K4kXtRkZgY8B7zn96AJ3q0n59xsYLaZzQV8Gzg9WkdtgSSgCXDIzOY55/IimftEPP4s7QXKRyRoEXi0nnYA2fl3cyMY96S8WkdmVg/I9HvQBM/W0UDgd/nLmg6oqfWZxmaNzX7Q2Kyx2Q8amyM/NpeKprYQdYDtBe7vAC4pbGYzqwk8A1xkZo/kr5xYU6SagPuAdkA1M2vknBsbyXDFVNT11BboQugLeV4kgxVTkepxzj0GYGa9yP8VNaLpiqeo66gLcDVQHfhzRJMVX1E/SzOBUWbWGlgayWDFVNR6APoQ241fUWuaDzxhZjcD/4lgLikZjc0am/2gsVljsx80Nns4NpfmprZInHPfETq+JW44514GXvY7h5ecc4uBxT7H8JxzbrzfGbzinJtJaKCJG865g4QGmrjhnPud3xm85Jz7HPB991TxlsbmYNDYHPs0NgeDxubClYoTRRXiG+CMAvfr5k8LMtUU++KtHlBNQRBv9UB81iTxuV5VU+yLt3pANQVBvNUDPtZUmpvaFUCymTUws0SgOzDb50wlpZpiX7zVA6opCOKtHojPmiQ+16tqin3xVg+opiCIt3rAz5qcc3F/A94CdvK/U3r3yZ9+DaEzdG0CHvM7p2qKr5rirR7V5H/W0lhPvNakW3yuV9UU+7d4q0c1+Z+1NNYTizWVikv6iIiIiIiISHwqzbsfi4iIiIiISMCpqRUREREREZHAUlMrIiIiIiIigaWmVkRERERERAJLTa2IiIiIiIgElppaERERERERCSw1tSIiIiIiIhJYampFREREREQksNTUipQCZtbBzNbk35abmT77IiIiIhIXzDnndwYRiTAz+xpo45zb6XcWEREREREvaWuNSOkwD1hrZi/6HURERERExEtl/Q4gIpFlZr8ADDjdOZfjdx4RERERES9pS61I/LsR2OCcy7GQqn4HEhERERHxio6pFYlzZtYSeB1wwCGgn3Nulb+pRERERES8oaZWREREREREAku7H4uIiIiIiEhgqakVERERERGRwFJTKyIiIiIiIoGlplZEREREREQCS02tiIiIiIiIBJaaWhEREREREQksNbUiIiIiIiISWGpqRUREREREJLD+DyhXz9zxz7TbAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 1152x432 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def CG_cond(A, b, eps, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Conjugate Gradient algorithm estimate solution x to the problem\n",
|
|
" Ax = b, after diagonal conditioning A, and the number of iterations k it\n",
|
|
" took to decrease maximum norm error below eps or to exceed iteration maximum\n",
|
|
" k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" D = np.diag(np.diag(A))\n",
|
|
" C = np.sqrt(linalg.inv(D))\n",
|
|
" A_tilde = C @ A @ C\n",
|
|
" b_tilde = C @ b\n",
|
|
" \n",
|
|
" x_tilde, k = CG(A_tilde, b_tilde, eps, k_max)\n",
|
|
" x = C @ x_tilde\n",
|
|
" return x, k\n",
|
|
"\n",
|
|
"# Sorry for copying.\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots(1, 2, sharex = True, figsize = (16,6))\n",
|
|
"\n",
|
|
"for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\"), (CG_cond, \"Conjugate Gradient\")]:\n",
|
|
" k_list = []\n",
|
|
" t_list = []\n",
|
|
" for eps in eps_list:\n",
|
|
" t_start = time.time()\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" t_list.append(time.time() - t_start)\n",
|
|
" k_list.append(k)\n",
|
|
" ax[0].plot(eps_list, k_list, label=alg_name)\n",
|
|
" ax[1].plot(eps_list, t_list, label=alg_name)\n",
|
|
" \n",
|
|
"ax[0].set_xscale(\"log\")\n",
|
|
"ax[0].invert_xaxis()\n",
|
|
"ax[0].set_yscale(\"log\")\n",
|
|
"ax[0].set_xlabel(\"$\\epsilon$\")\n",
|
|
"ax[0].set_ylabel(\"$k$\")\n",
|
|
"ax[0].legend()\n",
|
|
"\n",
|
|
"ax[1].set_yscale(\"log\")\n",
|
|
"ax[0].set_xlabel(\"$\\epsilon$\")\n",
|
|
"ax[1].set_ylabel(\"runtime (seconds)\")\n",
|
|
"ax[1].legend()\n",
|
|
"\n",
|
|
"fig.show()\n",
|
|
"\n",
|
|
"print(\"The number of iterations is brought down by a lot due to the conditioning, bringing down the runtime.\")"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.8.10"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|