1099 lines
114 KiB
Plaintext
1099 lines
114 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",
|
|
" 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",
|
|
" k_list = []\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" k_list.append(k)\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(\"$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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACdnElEQVR4nOzdd3Rc1b328e/WqPdiyZYtS5Yt9y7JDRswxWCKA6E7hQ4hCSnkDSmX9HKTm+QmuSEkxIAhkEAgJqEXU0Mz2HIBXHDvTbZ615T9/nFGI8mWu6SjkZ7PWrNm5syZ0W8o2nrObsZai4iIiIiIiEg4inC7ABEREREREZGTpVArIiIiIiIiYUuhVkRERERERMKWQq2IiIiIiIiELYVaERERERERCVsKtSIiIiIiIhK2It0uoDP069fPDhkyxO0yRESkl1i+fPlBa22m23WEM7XNIiLSmY7WNveKUDtkyBBKSkrcLkNERHoJY8x2t2sIV8aYecC8goICtc0iItJpjtY297jhx8aY2caYt40x9xpjZrtdj4iIiBw/a+2z1tpbU1JS3C5FRET6iG4JtcaYhcaYUmPM6kOOzzXGrDfGbDLGfCd42AK1QCywqzvqExERERERkfDUXT21DwFz2x4wxniAe4ALgDHAfGPMGOBta+0FwLeBH3dTfSIiIiIiIhKGuiXUWmvfAsoPOTwV2GSt3WKtbQb+AVxirQ0EX68AYrqjPhEREREREQlPbi4UNQjY2eb5LmCaMeYy4HwgFfjjkd5sjLkVuBUgNze366oUERERERGRHqvHrX5srf0X8K/jOG8BsACguLjYdnVdIiIiIiIi0vO4ufrxbmBwm+c5wWMiIiIiIiIix8XNULsMGG6MyTfGRAPXAM+cyAcYY+YZYxZUVVV1SYEiIiIiIiLSs3XXlj6PAUuAkcaYXcaYm6y1PuB24GVgHfCEtXbNiXyu9sITERERERHp27plTq21dv4Rjr8AvNAdNYiIiIiIiEjv4+bw41Om4cciIiIiIiJ9W1iHWg0/FhERAAJ+qNkHfp/blYiIiAhAIADVe5w2uov1uC19RERE2vH7oHaf0zBW7w7eH/K4Zi8EfHB7CfQb7nbFIiIivZ+1UHcQKrc7t4rgfeUO53HVTvA3w1dXQXp+l5YS1qHWGDMPmFdQUOB2KSIicjJ8zU4g7TCwBp/X7gcbaP++qHhIHujchswKPh4EcenufA8REZHexlporDw8rFbuaH3urW//nvgMSM2FAeNh9MXO49iuH1Ub1qHWWvss8GxxcfEtbtciIiKH8Da0CalH6GWtKz38fTHJrYE1a7QTVltCa8vx2FQwptu/Ul9ljEkA/gQ0A29aa//uckkiItIZmmrb9LLuOPxxU3X782NSIC0XMgpg2DlOaE3Lg9Q8SB0MMUmufI2wDrUiIuKS5rpDgmqbwFoV7GVtKD/8fbGpreE0e2L7oNryODa5279OX2SMWQhcDJRaa8e1OT4X+D/AA9xvrf0lcBmwyFr7rDHmcUChVkQkHHgboHLnkYcIH9pWR8UHA2ou5M1ofZwWvI9Lc+d7HINCrYiItNdYdfT5q9W7nXMOFZ/hhNKUHBg89ZDe1UGQnA3RCd3/feRIHgL+CDzccsAY4wHuAeYAu4BlxphngBzg4+BpXb/ih4iIHB9fM1TvOvIQ4dr97c/3RDvhNDUPsie1htXUIc7j+IywHAkV1qFWc2pFRE6Qr9lZuKF8C1Tt6ji0Ntcc8iYDiVlOOE0f2n4Oa0sva9JAiIp15SvJybHWvmWMGXLI4anAJmvtFgBjzD+AS3ACbg6wijDfOUFEJKwE/E7bfKQhwjV72q87YTzOxeW0PBg+pzWstgTZxP4Q0ft+jYd1qNWcWhGRDngboWKbE1wPvVXtPKTxi4DEAU4wzRwJw84O9rYOag2tiQMgMtq1ryPdahCws83zXcA04A/AH40xFwHPHunNxphbgVsBcnNzu7BMEZFepKESyjdD2Rao3Na+17Vql7O6f4hx2ufUXMg/vTWstgTXpIHgCeuId1L63jcWEekNmmqhYmubwLq19b56N2Bbz41NgfRhkDMFJlzt9Lam50PKYOeKbR9s/OTEWGvrgBuO47wFwAKA4uJie4zTRUT6juY6KNscDK+b2z+uP9j+3MT+TkAdVAzjLm8Nrqm5TtutC82H0V8yIiI9VUOlE1Qrth4SXLccPkcmIdMJq/mnQ1p+MLgGw2u8trmR47YbGNzmeU7wmIiIHIu30WmzOwqvNXvbn5uU7awgPOoiyBjmPE4f5vS4RsW5U38YU6gVEXGLtVBf3vEw4fIth69ImDTQCanD57QJrUOdEKsVg6VzLAOGG2PyccLsNcBnTuQDtN6FiPRqfq8zLLhsU5vwuskZOly1k3YjpeL7OYF12NlOex0Kr0O1cGInC+tQq4ZTRHo8a51e1Q6D6zZoaruKsHGGFaXnw5hLDgmuQyA63qUvIb2RMeYxYDbQzxizC/ihtfYBY8ztwMs4W/ostNauOZHP1XoXIhL2An5nLuthQ4U3OfNdbZtF4GNSnLCaOw0yPuv0tmYMde7jUl37Cn1NWIdaNZwi0iMEAs481raBtWJr63Bhb33rucbjDC1KHwqDp7X2tKYPdY5Hxrj3PaRPsdbOP8LxF4AXurkcEZHuZa0zJPjQ3tbyzU7b7W9uPTcq3gmuAybA2E+3DhXOGBa2W+D0NmEdakVEuo3f17oVzqHzWyu2gb+p9VxPdDCo5kP+Ga1zW9OHOj2xnijXvoZIV9MoKhHpMayFuoNtelw3ta4yXL4FvHWt53pinLY6owCGn9d+nmvSAAXXHk6hVkTkUI1VsO9j2Pth661sU/sl9SPjnJDabziMOL/9UOHkgRDhca9+ERdpFJWIdLuGykMWZ2oTXttO8zEeZzpPxjBnz/WMYa3hNXmQ2u4wplArIn1bXRns+7B9gC3f0vp60kDInggjL3QavpahwrpqKyIi0r18TbB/NexZBXtXQeknTnitL2tzUnB9ioxhMOHK4DDhAud5aq5GS/VSCrUi0nfU7GsfXvd+GFypMCg1zwmwkz4L2ZMgewIkZrlWroiISJ91aIDdswpK17aOmopNhf7jYNTFTmBtCa9pQyAq1rWyxR1hHWo1b0dEOmStE1YPDbChvV2N0/ANngZTb3WCbPYEiEtztWyR3kBts4icsOMJsAMnw2lfcS46D5zkXIjWiCkJMtbaY5/VwxUXF9uSkhK3yxARNwQCzkrDe1e1D7ANFc7rxgOZo4LBNXgbMA5iklwtW3o2Y8xya22x23WEM7XNItKh4w2wAycpwEo7R2ubw7qnVkT6GL8PyjYe0gP7ETTXOK97oiFrDIz+VDDAToL+YyAqztWyRURE+qQOA+w6CHid19UDK51EoVZEeiZfMxxY1z7A7lsNvgbn9ah4Zy7NxGtae2AzR0FktLt1i4iI9EXHHWBvV4CVTqdQKyLu8zbA/jXthxDvX9vaEMYkO6F1yk2tATajQEvvi/RAmlMr0gcowEoPo1ArIt2rqebwPWAPrAfrd16PS3cav9Nubw2wqUMgIsLNqkXkOGmfWpFeRgFWwoBCrYh0nfpy2PdR+wBbtqn19aRsJ7SOntcaYJMHqSEUERFxg6/JGTm1Z6UCrIQVhVoR6Tw1++GT52DLG06ArdzR+lpqrhNaJ17jNIQDJkBSf9dKFRER6dMUYKUXCetQq3k7Ij1A5U5Y9yysewZ2vA9Yp9HLmQJTbg5uoTMB4tPdrlRERKRv8vuckVMKsNJLhXWo1bwdEZeUbYa1TztBds9K51j/cTD7O852Olmj1RCK9FG64CzSQwT8sP09WPNvp72uO+AcV4CVXiisQ62IdBNrnY3R1z7j9MqWrnGODyyEc3/kBNmMYa6WKCI9gy44i7go4IcdS5wgu/YZqCt1tsAbcT6MuhhyihVgpVdSqBWRjlkLe1Y4IXbtM1C+GTCQOwPm/tJpHFMHu12liIhI3xYIwM73g0H2aajdD5FxMOI8GPtpGH4eRCe4XaVIl1KoFZFWAT/sXOoMU1r3LFTtBOOB/DNgxpedIKvFnURERNwVCMCupa1BtmYvRMY6AbYlyMYkul2lSLdRqBXp6/xe2PaOE2Q/ed65wuuJhmFnw+zvwsgLtMiTiIiI2wIB2F3iBNk1T0HNHvDEwPA5TpAdMVdBVvoshVqRvsjXBJvfcILs+hegocKZczN8jjM/dvh5EJvsdpUiIiJ9m7Wwe3lrkK3e5Vx4LpgDY3/izJVVey2iUCvSZzTXwaZXnfmxG16G5hqISYGRc2H0PBh2DkTHu12liIQ5rX4scoqshd0rYG0wyFbtDI6gOgfO+YHTbsemuF2lSI+iUCvSmzVWOQF27dOw6TXwNUBcOoy9FMZcAvlnQmS021WKSC+i1Y9FToK1zhZ5a/4Na5+Cyh0QEeVMBTrrLmcqUFyq21WK9FgKtSK9TV0ZrH/e6ZHd8qazsXriAJj8ORjzKcg9DTz6X19ERMRV1sLeD4NDi/8NldshIhKGngVnfgdGXQhxaW5XKRIW9JetSG9QvRc+ec6ZI7vtXbB+SM2FaV9wemQHFUNEhNtVioiI9G3Wwr6PW4NsxVZnl4Ghs+GMO2HURVqcUeQkhHWo1bwd6dMqtjvb7qx7xtmGBwsZw2HW153FnrInanN1ERERt1kL+9e0Btnyza3b5c26w1nXQkFW5JSEdajVvB3pcw5udObHrnvGGbIE0H88nPVfTpDNGuVufSIiIuIE2dJ1rUG2bCOYCCfInvYVJ8gm9HO7SpFeI6xDrUivZy3sX+3Mj133LBxY5xwfVAxzfuI0iulD3a1RREREHKXrnBWL1/wbDq53guyQWTDjSzBqHiRmul2hSK+kUCvS07TsSbfuGSfMVmx1GsXc0+CCXznzbVJy3K5SREREAA6sbw2yB9YBxgmy0251RlElZrldoUivp1Ar0hME/LDjfSfIrnsWqnc7KyDmnwEzvwajLtbVXREJC1rvQvqEgxtbhxaXrgUM5J0GF/7GGUWVNMDtCkX6FIVaETd5G+CVH8Kaf0HdAfDEQME5cPb3nc3VtZS/iIQZrXchvVbZZqe9XvOUMzUIIHeGM4pq9KcgOdvV8kT6MoVaEbcEAvCvW52e2TGXOHvIDj8PYpLcrkxERETACbJrn3J6ZPd97BwbPA3m/tIJsimDXC1PRBwKtSJuee3HznDj834Op93udjUiIiIC4PfBB/fCx0+07jSQMwXO/2/nIrTWtRDpcRRqRdyw/CF49/dQfBPM+LLb1YiIiEiLV34A798DAwvhvJ85QTY11+2qROQoFGpFutvm1+G5b0DBuc48HGPcrkhEREQAVjzsBNppt8EF/+N2NSJynCLcLkCkT9m/Fp64DrJGwxUPgkfXlURERHqE7e85F52HnuVMDRKRsKFQK9JdavbDo1dDVDx85nGITXa7IhEREQGo2A6Pfw7S8uBKXXQWCTf6P1akOzTXw2PXQP1BuOEFLTIhIiLSUzTVwGPzIeCD+Y9rOz2RMKRQK9LVAgH41y2wZyVc8ygMnOx2RSIiIgKt2+sd+AQ+twj6FbhdkYichB45/NgYk2CMKTHGXOx2LSKn7NUfwCfPwdxfwKgL3a5GREREWrz+U1j/gtNGDzvb7WpE5CR1S6g1xiw0xpQaY1YfcnyuMWa9MWaTMeY7bV76NvBEd9Qm0qVKFsJ7d8OUW5yVFEVEejljzDxjzIKqqiq3SxE5uo+egHd+C0XXw9Rb3a5GRE5Bd/XUPgTMbXvAGOMB7gEuAMYA840xY4wxc4C1QGk31SbSNTa9Cs9/E4afB3N/qa17RKRPsNY+a629NSUlxe1SRI5s13J4+nbImwUX/FpttEiY65Y5tdbat4wxQw45PBXYZK3dAmCM+QdwCZAIJOAE3QZjzAvW2kB31CnSafavgSeuh6wxcMVCraIoIiLSU1Tthn/Mh6QBcNXDEBntdkUicorc/Et7ELCzzfNdwDRr7e0AxpjrgYNHCrTGmFuBWwFyc3O7tlKRE1GzD/5+FcQkOlv3xCS5XZGIiIiAsxvBPz4DzXVw7dOQkOF2RSLSCXps95G19qFjvL4AWABQXFxsu6MmkWNqrnO27mmoCG7dM8jtikRERATAWnj6S7D3Q5j/D8ga7XZFItJJ3Ay1u4HBbZ7nBI+JhKeAH568xWksr3kMBk5yuyIRERFp8davYc2/Yc5PYOTcY58vImHDzS19lgHDjTH5xpho4BrgmRP5AK2wKD3KKz+A9c87i0KpsRQREek51j4Nb/wcJs6H077qdjUi0sm6a0ufx4AlwEhjzC5jzE3WWh9wO/AysA54wlq75kQ+VyssSo+x7H5Y8keY+gWY9gW3qxEREZEWez+Ef98GOVPh4t9rpWORXqi7Vj+ef4TjLwAvdEcNIl1m4yvwwp0wYq6zebuIiIj0DDX74bHPQFw6XP03iIp1uyIR6QJuDj8+ZRp+LK7b9zH883roPw4ufwAiPG5XJCIiIgDeRnj8s9BQDvMfhaT+blckIl0krEOthh+Lq6r3wqNXQ0xycOueRLcrEhEREXBWOn7u67BrGXz6Xsie6HZFItKFeuyWPiI9WnMdPHY1NFTCjS9B8kC3KxIREZEW7/0BPnwMzroLxlzidjUi0sXCuqdWw4/FFQE/PHmzM/T4ygche4LbFYmIiEiL9S/BKz+EsZfBGXe6XY2IdIOwDrUafiyuWPw9WP8CXPArGHG+29WIiIhIi/1r4cmbnOHGl9yjlY5F+oiwDrUi3e6DBfD+n2DaF2HqLW5XIyIiIi3qyuCxayA6EeY/BtHxblckIt1EoVbkeG14GV76Noy4AM7/udvViIh0K2PMUGPMA8aYRW7XInIYXzM8cS3U7INrHtVaFyJ9TFiHWs2plW6z9yNYdCMMGA+X36+te0QkrBhjFhpjSo0xqw85PtcYs94Ys8kY852jfYa1dou19qaurVTkJFgLL3wTtr/jDDnOKXK7IhHpZmEdajWnVrpF9R5n657YFJivrXtEJCw9BMxte8AY4wHuAS4AxgDzjTFjjDHjjTHPHXLL6v6SRY7T0gWw4q9w+v+DCVe6XY2IuEBb+ogcTVMtPHoVNFXDjS9DcrbbFYmInDBr7VvGmCGHHJ4KbLLWbgEwxvwDuMRa+wvg4m4uUeTkbHoNXvoOjLwIzvqe29WIiEvCuqdWpEsF/M4KivvXwJUPwYBxblckItKZBgE72zzfFTzWIWNMhjHmXmCyMea7RzjnVmNMiTGm5MCBA51brcihDm6Ef94AmaPhsgUQoT9rRfoq9dSKHMnL/wUbXoKL/heGz3G7GhERV1lry4DbjnHOAmABQHFxse2OuqSPaqhwVjr2RDorHWtqkEifFtaXtLRQlHSZ9++FD+6FGbfDlJvdrkZEpCvsBga3eZ4TPCbSs/l9Tg9txXa4+m+Qlud2RSLisrAOtVooSrrE+hfh5e8683Pm/MTtakREusoyYLgxJt8YEw1cAzxzqh+qC87S5RbfBVvegIt/C3mnuV2NiPQAYR1qRTrd3g9h0U0wYAJcfp+27hGRXsEY8xiwBBhpjNlljLnJWusDbgdeBtYBT1hr15zqz9IFZ+lSJQ86I6mmfxkKr3W7GhHpITSnVqRF1W5n6564NPjM4xCd4HZFIiKdwlo7/wjHXwBe6OZyRE7Otnec/WgLztVIKhFpRz21IgBNNU6gbaqFzz4BSQPcrkhEJCxp+LF0ifKt8PjnIX0oXLHQWSBKRCRIoVbE74NFN0LpWrjqIeg/1u2KRETCloYfS6drrIbH5oMNwPx/QKz+2xKR9sI61OpqsJwya51N2zcuhot+4wxpEhERkZ4h4Id/3QIHN8BVf4WMYW5XJCI9UFiHWl0NllP2wb2w7D447StQfKPb1YiIiEhbr/3E2TP+gv+BobPdrkZEeqiwDrUip+ST5+Gl78Koi+FcLTghItIZNIpKOs2H/4B3fw/FN8HUW9yuRkR6MIVa6Zv2rIQnb4aBk+Gy+yBC/yuIiHQGjaKSTrFzGTzzFRhyutNLKyJyFPpLXvqeql3w6DUQn+EsOBEd73ZFIiIi0qJqF/zjM5A8CK56GDxRblckIj2c1kOXvqWxGv5+FXjr4fMvQ1J/tysSERGRFs11zkrHvka4/jmIT3e7IhEJA+qplb7D74NFN8CBT5wVFPuPcbsiEZFeR3Nq5aQFAvDv22D/amcv2syRblckImFCoVb6BmvhxTth06tw8W9h2NluVyQi0itpTq2ctP/8D6x7Bub8FIbPcbsaEQkjYR1qdTVYjtuSe6BkIcz8GhRd73Y1IiIi0tbqf8F/fgmTPgczvux2NSISZsI61OpqsByXdc/B4u/B6E/BOT9yuxoRERFpa89KeOpLMHi6M5rKGLcrEpEwE9ahVuSYdq9wtu4ZVAiXLdDWPSIiIj1JzT547DOQ0A+u/htExrhdkYiEIf2FL71X5Q547BpIyHS27omKc7siEZFeT1OD5Lh5G5ytexqrYP5jkJjpdkUiEqYUaqV3aqyCR68GbyN89p+QmOV2RSIifYKmBslxsRae+SrsXu6MpBow3u2KRCSMaZ9a6X38Xvjn9XBwA3x2EWSNcrsiERERaeud38HHT8DZ34fRF7tdjYiEOYVa6V2shRfuhM2vw7w/wLCz3K5IRERE2vrkeXjtJzDuCjj9/7ldjYj0Ahp+LL3Le3fD8gdh1h1QdJ3b1YiIiEhb+1bDk7fAwMlwyR+10rGIdAqFWuk91j4Dr/wAxlwKZ//A7WpERESkrbqD8Nh8iE2Gax7VAo4i0mkUaqV32LUc/nUr5BTDp+/V1j0iIi7R6sfSIV8zPP45qCt1Am1yttsViUgvor/8JfxVbIfHrna2ArjmMV35FRFxkVY/lsNYC8/fATuWwKV/cvaOFxHpRGEdanU1WGiohEevcq4Af3aR9rgTERHpad7/M6z8G5zxLRh3udvViEgvFNahVleD+zi/F/55HZRtgqsfgcyRblckIiIibW18FRbfBaPnwezvul2NiPRS2tJHwpO18Pw3YMubcMk9MPRMtysSERGRtg5sgEU3QNZY+PRftN6FiHQZ/XaR8PTu/8GKh5397SZ/zu1qREREpK36cme9i8gYmP8YRCe4XZGI9GLqqZXws+YpePWHMPYyOOt7blcjIiIibfm98M/roWoXXPccpA52uyIR6eUUaiW87FwG//4C5EyFS/+soUwiIiI9zUvfha3/cdrp3GluVyMifYASgYSPgB+evBES+ztDmaJi3a5IREQOoZ0J+rhl98Oy++C0r8Kkz7hdjYj0EQq1Ej42vAyVO+C8n0JCP7erERGRDmhngj5sy3/ghW/B8PPh3B+5XY2I9CEKtRI+lt0HSdkw8iK3KxEREZG2yjbDE9dCv+Fw+f0Q4XG7IhHpQzSnVsJD2WbY/DrM/i/w6D9bEel7jDERwERgINAArLbWlrpblQjQWAWPzQcTAfP/AbHJblckIn2M0oGEh2UPQEQkFF3ndiUiIt3KGDMM+DZwLrAROADEAiOMMfXAX4C/WmsD7lUpfdrTX4byzfD5pyA93+1qRKQPUqiVnq+5Hlb9DUZ/CpIGuF2NiEh3+xnwZ+AL1lrb9gVjTBbwGeDzwF9dqE36ugMbYN2zzkiq/NPdrkZE+iiFWun5Vi9yhjZNudntSkREup21dv5RXisFft991YgcYvlDEBEFxTe6XYmI9GFaKEp6Nmth6X2QNQbyTnO7GhER1xhjrjTGJAUff98Y8y9jTKHbdUkf5m2AVX+H0fMgMdPtakSkD+txodYYM9oYc68xZpEx5otu1yMu21UC+z6CKTeBMW5XIyLipu9ba2uMMbOAc4AHcIYli7hjzVPQWKleWhFxXbeEWmPMQmNMqTFm9SHH5xpj1htjNhljvgNgrV1nrb0NuAqY2R31SQ+27D6IToIJV7tdiYiI2/zB+4uABdba54FoF+uRvq5kIWQMhyGz3K5ERPq47uqpfQiY2/aAMcYD3ANcAIwB5htjxgRf+xTwPPBCN9UnPVHdQVjzb5g0H2KS3K5GRMRtu40xfwGuBl4wxsTQA0dcSR+xbzXsWgrFN2gklYi4rlsaQ2vtW0D5IYenApustVustc3AP4BLguc/Y629APhsd9QnPdSKh8HfDMU3uV2JiEhPcBXwMnC+tbYSSAfudLWiDhhj5hljFlRVVbldinSl5Q+CJwYmHnEdMxGRbuPmFd5BwM42z3cBg4wxs40xfwhejT5iT60x5lZjTIkxpuTAgQNdXat0t4AfSh6EIadD1ii3qxERcY0xJt0Yk46zN+2bQFnweRNQ4mZtHbHWPmutvTUlJcXtUqSrNNXCh4/DuMsgPt3takREet6WPtbaN3Ea7WOdtwBYAFBcXGyPcbqEm42LoWoHnPdTtysREXHbcsACBsgFKoKPU4EdQL5rlUnftHoRNNdogSgR6THc7KndDQxu8zwneEzE2cYnKRtGXeR2JSIirrLW5ltrhwKvAvOstf2stRnAxcBid6uTPqlkIfQfBzlT3K5ERARwN9QuA4YbY/KNMdHANcAzJ/IBmrfTS5Vths2vQdH14IlyuxoRkZ5iurU2NC3HWvsioA28pXvtXgF7P3TaaC0QJSI9RHdt6fMYsAQYaYzZZYy5yVrrA27HWfRiHfCEtXbNiXyu5u30UiULISLSaTBFRKTFHmPM94wxQ4K3u4A9bhclfUzJQohK0FZ7ItKjdMucWmtth0vjBa84a9seadVcDyv/BqPnQdIAt6sREelJ5gM/BP4dfP5W8JhI92iohNVPwvgrITbZ7WpEREJ63EJRJ8IYMw+YV1BQ4HYp0llWPwmNlTDlZrcrERHpUay15cDX3K5D+rCPngBvvbM3rYhIDxLWm7Zr+HEvYy0suw8yR0PeTLerERHpUYwxI4wxC4wxi40xr7fc3K5L+ghrnb1pB052biIiPUhY99RKL7N7ubP4xIW/0eITIiKH+ydwL3A/4He5Fulrdn4ApWvhU3e7XYmIyGHCOtRq+HEvs/Q+iE6Cide4XYmISE/ks9b+2e0ipI8qWQgxyTDucrcrERE5jIYfS89QVwZr/uUE2pgkt6sREemJnjXGfMkYk22MSW+5uV2U9AF1ZbDmKaeNjk5wuxoRkcOEdU+t9CIrHwZ/M0y5ye1KRER6quuC93e2OWaBoS7UIn3Jh4+CvwmKtECUiPRMCrXivoDfGdY05HTIGu12NSIiPZK1Nt/tGqQPshZKHoTB06H/GLerERHpUFgPP5ZeYuMrULlDvbQiIkdhjIkyxnzVGLMoeLvdGBPldl3Sy219C8o3Q/GNblciInJEYR1qjTHzjDELqqqq3C5FTsWy+yApG0Zd7HYlIiI92Z+BIuBPwVtR8JhI1ylZCHFpMOYStysRETmisA61WiiqFyjfAptehaLrwaMOBxGRo5hirb3OWvt68HYDMKW7frgx5lJjzH3GmMeNMed1188VF9Xsh0+eg0mfhahYt6sRETmisA610gssewAiIqHwumOfKyLSt/mNMcNanhhjhnKc+9UaYxYaY0qNMasPOT7XGLPeGLPJGPOdo32GtfYpa+0twG3A1SdRv4SbVX+DgE8LRIlIj6eFosQ93gZY+Tdn2HFyttvViIj0dHcCbxhjtgAGyAOON208BPwReLjlgDHGA9wDzAF2AcuMMc8AHuAXh7z/RmttafDx94Lvk94s4IflD0H+GdCvwO1qRESOKqxDrTFmHjCvoEC/bMPS6iehsRKm3Ox2JSIiPZ619jVjzHBgZPDQemtt03G+9y1jzJBDDk8FNllrtwAYY/4BXGKt/QVw2CIHxhgD/BJ40Vq7oqOfY4y5FbgVIDc393hKk55q8+vOIo5zfuJ2JSIixxTWw481pzaMWQtL74PM0TBkltvViIj0eMaYLwNx1tqPrLUfAfHGmC+dwkcOAna2eb4reOxIvgKcC1xhjLmtoxOstQustcXW2uLMzMxTKE1cV7IQErJg5EVuVyIickxhHWoljO1eAXtXOdv4GON2NSIi4eAWa21lyxNrbQVwS3f9cGvtH6y1Rdba26y193bXzxUXVO2CDS9B4echMtrtakREjkmhVtyx7D6IToQJWmtEROQ4eYJDgIHQnNhTSRy7gcFtnucEj50SbbfXC6x4xBlRpUUcRSRMKNRK96srg9X/gonXQGyy29WIiISLl4DHjTHnGGPOAR4LHjtZy4Dhxph8Y0w0cA3wzKkWqalBYc7vgxV/hYJzIS3P7WpERI6LQq10v5WPgL9JC0SJiJyYbwNvAF8M3l4DvnU8bzTGPAYsAUYaY3YZY26y1vqA24GXgXXAE9baNV1SuYSPDS9BzV4ovtHtSkREjptWP5buFfA7i0/kzYKs0W5XIyISNqy1AWPMQ8Dr1tr1J/je+Uc4/gLwQieUJ71FyUJIHgTDz3O7EhGR4xbWPbUa4hSGNr0KldudBaJEROS4GWM+BawiOOTYGDMpuK9sj6I5tWGsfCtsfg0KrwVPWPd7iEgfE9ahVsLQ0vsgcQCMnud2JSIi4eaHOHvLVgJYa1cB+S7W0yFdcA5jK/4KxuOEWhGRMKJQK92nfIvTU1t0PXii3K5GRCTceK21h3Z/Wlcqkd7H1+ysejzyAkge6HY1IiInRKFWuk/JQjARUKQtAkRETsIaY8xncLb2GW6MuRt4z+2iDqXhx2Hqk2eh/iAU3+B2JSIiJ0yhVrqHtwFW/g1GX6wrwCIiJ+crwFigCWc7n2rg624W1BENPw5TJQ9Cah4MPdvtSkRETphWAZDusfpf0FChbXxERE6StbYeuAu4yxjjARKstY0ulyW9wYENsO1tOOeHEKH+DhEJP/rNJd1j2f2QOQqGnO52JSIiYckY86gxJtkYkwB8DKw1xtzpdl3SCyx/CCKiYPLn3K5EROSkhHWo1bydMLF7OexZ4fTSGuN2NSIi4WqMtbYauBR4EWfl48+7WpGEP28DrPq7sytBYpbb1YiInJTjDrXGmJ91cMzTueWcGM3bCRNL74foRJhwtduViIiEsyhjTBROqH3GWuulB65+rAvOYWbNU9BYCcU3ul2JiMhJO5Ge2kHGmPktT4wxWcCrnV+S9Cr15bD6SSfQxia7XY2ISDj7C7ANSADeMsbk4SwW1aPognOYKVkIGcNhyCy3KxEROWknslDUF4CXjTGbca4MPwh8u0uqkt5j5SPgb9ICUSIip8ha+wfgDy3PjTE7gLPcq0jC3r7VsGspnP/fmh4kImHtmKHWGPMwsAJYCXwZeBTwAZdaazd1bXkS1gIBWPYA5M2E/mPcrkZEJCwZYz4HPGqtDbQ9bq21gM8YMwzItta+40qBEr6WPwieGJg4/9jnioj0YMfTU/sQMBG4AZgADAGWAZ8zxqy21i7qsuokvG16FSq3w7k/crsSEZFwlgGsNMYsB5YDB4BYoAA4EzgIfMe98iQsNdXCh4/DuMsgPt3takRETskxQ6219nXg9ZbnxphIYDRO0J0GKNRKx5bdB4n9YdTFblciIhK2rLX/Z4z5I3A2MBPnAnMDsA74vLV2h5v1HcoYMw+YV1BQ4HYpcjSrF0FzjRaIEpFe4UTm1AJgrfXh7I/3MfC3Tq9IeofyrbDxFTjzWxAZ7XY1IiJhzVrrB14J3no0a+2zwLPFxcW3uF2LHEXJg5A1FnKmuF2JiMgpC+t9aqUHK1kIJgKKrne7EhEREWlr9wrYuwqKb9ACUSLSKyjUSufzNjirHo+6CJIHul2NiIiItFWyEKIStH+8iPQaCrXS+db8GxoqtI2PiIhIT9NQ6ewfP/4K7R8vIr1GWIdaY8w8Y8yCqqoqt0uRtpbdD/1GQv4ZblciItJrGGP6G2MeMMa8GHw+xhhzk9t1SZj56Anw1jtDj0VEeomwDrXW2mettbempKS4XYq02L0Cdi93emk1T0dEpDM9BLwMtMzr2AB83a1ijkQXnHswa529aQdOdm4iIr1EWIda6YGW3e/M05moeToiIp2sn7X2CSAAod0I/O6WdDhdcO7Bdn4ApWu1jY+I9DoKtdJ56sudeToTr4ZY/TEjItLJ6owxGYAFMMZMB9QdKsevZCHEJMO4y92uRESkU53wPrUiR7Tyb+Br1AJRIiJd4xvAM8AwY8y7QCZwhbslSdioL4c1T0HRdRCd4HY1IiKdSqFWOkcgACUPQO5p0H+s29WIiPQ61toVxpgzgZGAAdZba70ulyXhYtWj4G+CIi0QJSK9j0KtdI7Nr0HFNjjnB25XIiLSKxljPMCFwBCc9vs8YwzW2t+6Wpj0fNY6Q48HT4f+Y9yuRkSk0ynUSudYeh8kZMGoeW5XIiLSWz0LNAIfE1wsSuS4bH0LyjfDmd92uxIRkS6hUCunrmIbbFwMZ9wJkdFuVyMi0lvlWGsnuF2EhKGShRCXBmMucbsSETkKf8ASYcBoW8wTplArp65kIZgIKLre7UpERHqzF40x51lrF7tdyNEYY+YB8woKCtwuRQBq9sMnz8G02yAq1u1qRCTIWsu2snpW7axg1Y5KVu6sZN3eagyG1Pgo0uKjSY2PIj0hmtT4aNLaHEuLjyYtISp4PJqUuCg8EX07CCvUyqnxNsKKR2DUhZAyyO1qRER6s/eBfxtjIgAvzmJR1lqb7G5Z7VlrnwWeLS4uvsXtWgRY9TcI+LRAlIjLquq9rNpVycodFazaWcmHOyupqHfW+ouP9jB+UAo3zMzHGKis81JR30xlvZeNpbVUBh/7ArbDzzYGkmOjSIuPOiQAO49TEw4PxekJ0cRGebrzH0GXUqiVU7Pm39BQDlP0t4uISBf7LTAD+Nha2/FfNiJtBfyw/CHIPwP6qedcpLt4/QE+2VvDqp0VrNxRyaqdlWw5WAc4AXR4ViJzxvRn0uA0JuemMjwrkUhPxFE/01pLTZMvFHhbQq/z2Etlm/sDtU1s2O+E4bpm/xE/MzYqol34bRt6O+oRTouPIjk2ioge2CusUCunZtn90G+E02CKiEhX2gmsVqCV47b5dajcAXN+4nYlIr2WtZbdlQ2s2lkZGka8encVTT5nPb9+idFMGpzG5UU5TBqcyoScFJJio0745xhjSI51QmVuRvxxv6/J56eq3ktFfUvvbzPlda2P2wbidfuqqQw+P0KnMBEGUuIODcDBUJwQHDLdcizBeb1fYkyXD49WqJWTt2cl7C6BC37lXHYSEZGutAV40xjzItDUclBb+sgRlSx0diYYeZHblYj0GrVNPj7aVRnqgV21s5IDNc6v5OjICMYNTOZz0/OYNDiVSYNTyUmLc3Xhp5hID1nJHrKSj39OfSBgqWn0HaVHuPXx3qpG1u2tpqLeS4O3417hV79xBgVZSZ31lTrUI0OtMeZS4CIgGXigpy+K0WctvR+iEmDiNW5XIiLSF2wN3qKDN5Ejq9oNG16CmV/XzgQiJ8kfsGwsrXF6YIMhdkNpDS3jZfL7JTCroB+TBqcyOTeVUQOSiY48+jDicBARYUiJjyIlPoohJBz3+xq9fifw1rX2/lbUNzMgJa4Lq3V0W6g1xiwELgZKrbXj2hyfC/wf4AHut9b+0lr7FPCUMSYN+A2gUNvT1JfD6kUwcT7EprhdjYhIr2et/bHbNUgYWfEwWAtF17ldiUjYKK1uZGWw93XVjko+2lUZmpOaEhfFpMGpzB03gEm5qUzKSSUtQReM2oqN8pCdEkd2N4TYQ3VnT+1DwB+Bh1sOGGM8wD3AHGAXsMwY84y1dm3wlO8FX5eeZtXfwdcIU252uxIRkV7NGPN7a+3XjTHPAofNcrLWfsqFsqQn8/tgxV+h4FxIG+J2NSI9UqPXz+rdVaza2doLu7uyAYDICMPo7OTQPNhJg1PJ75eg/WN7sG4Ltdbat4wxQw45PBXYZK3dAmCM+QdwiTFmHfBL4EVr7YruqlGOUyAAyx6A3BkwYNyxzxcRkVPxSPD+N65WIeFjw0tQsxcu0nRrEXDmiG4tq2NVMLyu3FnBJ3trQlvkDEqNY1JuKjfMHMLk3FTGDkzpVdvd9AVuz6kdhLOaY4tdwDTgK8C5QIoxpsBae++hbzTG3ArcCpCbm9sNpUrI5tehYiuc/T23KxER6fWstcuDDydZa/+v7WvGmK8B/+n+qqRHK1kIyYNg+HluVyLiioq65uCesJWhPWGrGpw9YROiPUwcnMqtZwx1emFzU8lKOv5FlKRncjvUdsha+wfgD8c4ZwGwAKC4uFjbG3SnZfc7qymO1og3EZFudB3OGhRtXd/BMenLyrc6F59nfwc8PfLPPJFO1ewLsG5vdXAYcQWrdlayrawecLafGdE/iQvHDwgOI06jICuxy7eXke7n9m+73cDgNs9zgsekp6rY7gxrOuObWk1RRKQbGGPmA58B8o0xz7R5KQkod6eqIzPGzAPmFRQUuF1K37Tir2AioPBatysROWXWWmqbfBysbeZATRMHa5va3W/YX8PqPdU0B/eEzUqKYdLgVK6aMpjJg9OYkJNCQozbcUe6g9v/lpcBw40x+Thh9hqchvu4qOF0QclCZ0/aouvdrkREpK94D9gL9AP+t83xGuAjVyo6Cmvts8CzxcXFt7hdS5/ja4YVj8DICyB5oNvViBxRXZPvsIB6oE1wbftaozdw2PsjDKQnxDAkI57rZuQxaXAak3NTyU6J1WJOfVR3bunzGDAb6GeM2QX80Fr7gDHmduBlnC19Flpr1xzvZ6rh7GbeRlj5CIy8EFJy3K5GRKRPsNZuB7YDM9yuRXq4T56F+oNQfIPblUgf1NDsd8LoIWG19b41tNYHt8lpyxhIj4+mX2IMmUkx5OXFk5kUE3re9j49IVpDiKWd7lz9eP4Rjr8AvNBddcgpWPsU1JfBVF1DEBHpbsaYy4D/AbIAE7xZa22yq4VJz1HyIKTmwdCz3a5Eeokmn781jNY4gTV0f0hYrW3ydfgZafFRoUA6aXDqIQE1msykGDKDQTXSE9HN31B6C7eHH58SDT/uZsvuh4zhkH+m25WIiPRFvwLmWWvXuV2I9EAHNsC2t+GcH0KEgoEcWbMvQFldEwdrmjlQ2xi8bxkC3D68Vjd2HFSTYyND4XTswORQSM08pFc1PSGa6Ej99yhdL6xDrYYfd6M9q2DXMpj7P874EBER6W77FWjliJY/BBFRMPlzblciLgkELAfrmthf1cTeqgb2Vzeyr7qRfVVN7K9uZH91Iwdqm6is93b4/qSY1qA6ekAy/QqiOxz+m5EYTUyk9nCVniWsQ610o2X3QVQ8TLzG7UpERPqqEmPM48BTQFPLQWvtv1yrSHoGbwOs+juMngeJWW5XI12g0et3QmpVS1B17kPHqhoprWnCF2i/y6UnwpCVFMOAlFiGZSYyfWjGYUN/W57HRimoSvhSqJVja6iAjxc5gTYu1e1qRET6qmSgHjivzTELKNT2dWuegsZKKL7R7UrkBFlrqWrwsrclpB4SWvdVOcG1ooPe1YRoD/1TYslOiWX6sAwGJMcyICWW/snOsQHJsWQkxmhBJekTwjrUak5tN1n5d/A1wpSb3a5ERKTPstZqSVvp2PIHnTUvhsxyuxJpw+sPcKCmib1VrT2qLUOC2x5r8rXfssYYyEiIYUBKDDlpcRQPSWNAshNWBwRDbP/kWJJio1z6ZiI9T1iHWs2p7QaBAJQ8AIOnw4DxblcjItJnGWMexOmZbcdaq+65vmzfatj5AZz/31rzohvVNvlCIbVtQG3by3qwtgl7yP+x0ZERTo9qciwTclI5f2wwrCbHMiAlhv7JsWQlxWpxJZETFNahVrrBltehfAucdZfblYiI9HXPtXkcC3wa2ONSLdJTLH8QPDEwscOdE+UEWWuprPeyu7IhNCR4X1VDaLGlltDa0fY1qfFRoR7VMdnJoaHBbXtZ0+KjMLr4INLpFGrl6JY9AAmZzuITIiLdrNkXoLrRS3WDl+pGH1UNLY+9wce+0OtVDV5+dcUEslPi3C67S1hrn2z73BjzGPCOS+VIT9BUCx8+DuMug/h0t6sJC9Zaqht87KyoZ1dFPbsqGoK31seHBta2iy0Nz0pkVkG/dsOAW0JrXLQWWhJxS1iHWs2p7WKVO2DDSzDrGxAZ43Y1IhKG/AFLbUsYbRM+nce+NuHUCa2Hvt7g9R/186M8hpS4KJJjo0iOi6LRGzjq+b3McEBL3fZlq5+E5hoo0nTrtqoavIcF1p3lzv3uigZqDgmtiTGR5KTFkZMWz4xhGeSkxTMoNY6BqVpsSSRchHWo1ZzaLlay0Lkvut7VMkTEPdZa6pr9rWHzCOGzo3Ba3eA97I/HQ0UYSA6F0khS4qIoSEps97zldedxZCjApsRFERMZ0WeG8hljamg/p3Yf8G2XypGeoGQhZI2FwVPdrqRb1TR6QyH10J7WnRX11DS2/72TEO1hcHo8OWlxTB+aEQywTogdnBZPclxkn/k9ItJbhXWolS7ka4IVD8PICyF1sNvViMhJstbS4PWHgmdNY/shu8ca0lvd6MMfOGxtonYSY5zwmRQbSXJcFINS4xiTnRwKoK3BNDIURlueJ8boj8njYZx/SGOttTvcrkV6iN0rYO8quPA3vW6BqNomnxNSyxuCw4TbB9iqhvbb28RHe0IhdcqQNHLS4kPPc9LiSNU8VpFeT6FWOrbmKagv0zY+Ii4LBCy1zT5q2vR+Vjf6guG0Ta9oo9c555DQWtPow3eMUBobFdEufPZLjGZoZkKbYb2HhtPWXtPEmEgiPVqls6tZa60x5nnAlWXojTGjga8B/YDXrLV/dqMOaaNkIUQlwISr3a7khNU1+Q4Jqq29rLsqGqg8ZE/WuChPqHe1MDeNnLS4UM9rTlq8Fl8S6cFqmmtIik7q8p8T1qFWc2q70LL7IaMA8s90uxKRsObzB6ht8rUGzQ5C55GPOcN3D90S4lDx0Z5Q+EyKbQ2lbY+1DadJsc6xlmAaE6nFTcLECmPMFGvtshN5kzFmIXAxUGqtHdfm+Fzg/wAPcL+19pdH+gxr7TrgNmNMBPAwoFDrpsYqZz7t+CshNtntag5T3+xjd7BXtaOe1vK65nbnx0RGhALqpMGph/W0ZiREK7SKhAl/wM/qstW8u/td3t3zLqsPrubly19mQMKALv25YR1qNae2i+z9EHYthbm/hAj1wEjfFghYqhq8VLYZolvdEOwpbRNGa9r0mLY91tG2D4dKinGG5bYM3x2YGseo2KR2Q3Zbwmjb+adJwWNR6inFWsvW6q3kJ+f35j9+pwGfNcZsB+oAg9OJO+EY73sI+CNOGAXAGOMB7gHmALuAZcaYZ3AC7i8Oef+N1tpSY8yngC8Cj3TCd5FT8dET4K2HYncWiPL5A+ytamR7WT07yuvZXl7HrjZzXMsOCa3RbULruEEpTk9rm+DaL1GhVSSc7a/bz3t73uPdPe+yZM8SqpurMRjG9xvPrRNuJTKi6yNnWIda6SLL7oeoeO15J71SIGCpbPBSXtfEwdpmyuuaKatrpry2mbK6ptBj53gTFfXeo84pbbvQUVKs0xM6pF98aDGjlmNtA2rbY4kxkVpV8yQdqD/A+3vfd2573qe0oZSnL32aoSlD3S6tq5x/Mm+y1r5ljBlyyOGpwCZr7RYAY8w/gEustb/A6dXt6HOeAZ4JDoN+9GRqkU5grTP0eOBk59ZFapt8bC+rY2d5fSi8ttx2VzS0m9YQ5TGhXtXzBqYcshBTHP0SY4jQ7zmRXqPZ38zy/ct5b897vLP7HTZVbgIgMy6TswafxaxBs5iePZ3U2NRuq0mhVtprqICP/gkTroK4VLerETkmf8BSWd8aTstqmylvCafB52V1TaHHFfXNHCmjpsRFkZEQTXpCNHkZ8RTmpYWepyVEtQmnrb2mCdEe9TB0k3pvPSX7S3h/7/ss2bMk1IimxqQyLXsaM7JnkBGb4XKVXcdau70TP24QsLPN8104PcEdMsbMBi4DYoAXjnLercCtALm5uZ1Qphxm5wdQuhY+dfcpfUwgYCmtaWJ7WR3by+sPC6+HDhFOjY8iLz2e8YNSuHhCNrnp8eSmJ5CbEc+A5FhdnBPpxay1bK/ezrt73uXd3e9Ssr+EBl8DURFRFGYV8qmiT3HawNMYkTbCtb+JFGqlvVWPgq9BC0SJa/wBS0V9x4G0pfe05XF53XGE1MRoMhKiye+XQFFeOv0SnZCanhBNRkJM6PW0hGgN4+1h/AE/a8rWsGTPEt7f+z6rDqzCF/ARHRFNYf9CLh56MTMGzmBU+igijP7ddSVr7ZvAm8dx3gJgAUBxcfExZoPLSSlZCDHJMO7yY57a6PWzMxhSD+1t3VleT5OvdV/nCAOD0uLITY/n/LEDyE2PJy8jntz0eAanx5MSF9WV30pEepg6bx0f7P0gNDd2d+1uAHKTcrm04FJmDpzJlAFTiI+Kd7lSh0KttAoEYNkDMHgaZB9ripbI8fH5A1TUe0OBtDWsOj2q5XXNoWHALSH1SAsjpcZHBcNoNMMyE5mS7zzOSIgmPTEm1KuakRhNWrxCarix1rKzZidL9ixhyd4lLN23lJrmGgBGp4/m82M+z4zsGUzOmkxsZKzL1Ya93UDb/dpygsekJ6svd3YnKLwWohOw1lJW1+wE1bLW8LozOM91f3VTu7cnRHvIzUhgWGYCZ4/KYnB6PHnB8DowNU6/M0X6sIANsL58fag3dlXpKnzWR1xkHNOyp3H92OuZOXAmg5N75lafYR1qtfpxJ9vyBpRvhtnfdbsS6eGstVQ3+jhQ08j+6iZKW+6rm9hf08iBmibKap3AWtngPWJITQuF1BgKMhPJCIZUJ5gGQ2qwZ1UhtXeqaKzgg70fhIYU76nbA0B2QjZz8uYwI3sGU7Onkh6b7nKlvc4yYLgxJh8nzF4DfKYzPlhtc+fy+gPsrmhgR3k98cvvpdjfxI/3TmXJ799iZ3k9dc3+ducPSI4lNz2e04dnhnpbW8JrulYRFpE2yhvLWbJnCe/ufpf39rxHWWMZAKPSR3Ht2GuZNWgWkzInEeXp+SM1wjrUavXjTrbsAYjvB2M+5XYl4hJrLdUNPvbXNDoBtbqR0hrn/kBN6/PSmkYavYHD3h8f7aF/ciyZiTGM6J/UPpwGe1QzEmOCITVK+5v2QU3+JlbsX8GSvUt4f8/7fFL+CRZLUlQSU7OncuO4G5k+cDq5Sbn647uTGGMeA2YD/Ywxu4AfWmsfMMbcDryMs+LxQmvtms74eWqbT1xVg/eQOa11oV7XPZUNwSkWltej/8ZyM5K3q/uTlx7H9KEZoSHCeRnx5KTFExulLbpEpGO+gI+PDnzEO7vf4b0977G2bC0WS2pMKjMGzmDWoFmcNvA0+sX1c7vUExbWoVY6UeVO2PAizLoDImPcrkY6mbXOtjTtelWDwfXQ523nWLVIjIkkKymGrOQYJg1OpX9yDFlJsWQF7/snx5CVHEtijH6lSHstw5mW7F3Ckj1LWFm6kiZ/E5ERkUzMnMiXJ32ZGQNnMCZjTLcs+d8XWWs7XMreWvsCR1n0STrfnsoGlm+vYN3e6nbzWyvrve3Oy0iIJjcjnqK8NC6bPIjB6fGMa/6QoYv3Ebj0x7w6SXvIi8jx2VO7h3f3vMt7u9/j/b3vU+utxWM8TMicwJcmfYlZg2YxOn00nojwviCmvyDEUbLQuS9yZ887OTnWWirrvYf1rJa26WF1elabaO4grCbFRIaCaVFuGlnJscHw6tz3D94nKKzKCdhTuyc0nPiDvR9Q0VQBQEFqAVeOuJIZA2dQ3L+4xywuIZ1Lw48dTT4/a/ZUs2J7BSt2VLBieyX7qhsBiIwwoUWZLhqfHeptbVlNuMMLhE98D+LSiBh7afd+EREJK42+RpbvX847u9/h3T3vsrVqKwADEgZw/pDzmTloJtOyp5EcnexypZ1Lf6kK+JpgxcMw4gJI7ZmTv/saay0V9d6Oh/8G562WVjdxoKaJZn8HYTU2MhRIpwxJ7zCoZiXHEB+tXwFy6qqbq1m2d5kzpHjv+2yvdnaeyYzL5PSc05mePZ3p2dPJjM90uVLpDn11+HFpdSPLWwLsjko+3l0VupiYkxbH1Px0CnNTKcpLZ1R20omtEVBbCp88B9NugygtkiYiray1bK3aGhpSXLK/hCZ/E9ER0RQPKOaK4Vcwc9BMhqYM7dXTevQXrcDap6H+IEzVNj7dwR+wlNY0sruigV0VDeyubGBfVWNoGPCB4JxVr//w1ZVS4qJCwXRafkKbntU2YTUplrjo8B5CIj2b1+/lwwMfhubFri5bTcAGiIuMY8qAKVwz8hpmDJzR6xtQ6bu8/gDr9rb0wlayfHsFuysbAIiOjGD8oBSuP20IhbmpFAZHwZySlY9AwAdF15968SIS9qqbq9ttt7Ovbh8A+Sn5XDniSmYOmklR/yLiIuNcrrT7KNQKLLsf0odB/my3K+kV/AHL/upGdlU0sKui3gmuFQ3sqnQe76lsOCywpsa3htWhmQmhgNr2PjMpRguAiCustWyq3BTaL7Zl03WP8TCu3zhuGX8LMwbOYEK/CWGxQqLIiSqrbQqF1xU7KvhoV2VosbwBybEU5aVxw8whFOWlMWZgMjGRnfi7OuCH5Q9B/hnQb3jnfa6IhI2ADbC2bG0oxH504CP81k9iVCLTs6dz64RbmTlwJgMTB7pdqmsUavu6vR/Bzg/g/F9AhFaiPR4+f4B9odAaDKzB8Lqrsp69lY34Au1Da1ZSDDlpcUzISeXC8dnkpMWRkxZPTlocg1LjFFalxymtLw3Ni31/7/scbDgIwJDkIVwy7BJmDJzBlAFTSIpOcrlS6WnCfU6tP2BZv6+G5TsqWBkMsdvK6gGI8hjGDExh/tRcivLSKMxNY2BqF/eEbH4dKnfAnJ907c8RkR7lYMNB3tvzHu/sfof397wfWp9ibMZYbhx3I7MGzWJ85niiInQxGcI81IZ7w9kjLLsfIuNgUoeLY/ZJXn+AfVXte1pbHu+ubGBvVSP+NqHVmJbQGk9hbho5E53AOig1jpy0OAYqtEoYqPPWsXz/cpbscVYp3ly1GYD02HSmDZjGjIEzmJ49nezEbJcrlZ4u3ObUVtY3s3JHJSt2VLB8ewUf7qwM7f3aLzGGwtxU5k/NpTAvjfGDUrr/93nJQkjIhJEXde/PFemhqpqqeHX7q7y07SX21O4hwkSEbh7jaXcfOh5xhOOHnO8xHowxhx+P8GAIHo84ys851ucdo46ADbCqdBXv7nmXT8o/AZx2eNagWcwcNJMZA2do3/YjCOtQG24NZ4/TUAkf/xMmXAlxaW5X022afS2htb59cK10el33VrXsCegwxhlelpMWx5Qh6aGw2tLTmp0a27lDzUS6gS/gY/XB1aHe2I8OfITP+ojxxFDUv4hLCy5l+sDpjEgbQYTRKA7pHQIBy6YDtazYXhEaSrz5QB0AngjD6OwkLi/KoTA3jaK8NHLS4tydF161Gza8BDO/DpHR7tUh4rI6bx2v73idl7a9xHu738NnfeQm5TI2YywBAgRsAH/A3/rY+gkEWh/7A3681tv6Wpv7lttRjwece4t1fo5tfd7ZIk0kk7Im8bXCrzFz4ExGpo9UO3wcwjrUyila9Sh462FK77om0OTzs7eysV3vatvwuq+6Edvmd1BEKLTGMy0/vf3Q4LQ4slPiiI7UL5O+xlrLuvJ1vLL9FV7d/iq7ana5XVKnamn4DYbRGaO5bux1zBg4g0lZk4jxaK9q6R2qG72sCvbCrthRycodFdQ0+gBIi4+iMDeNywqdEDtxcErPWxF+xcNgLRRd53YlIt2u0dfIW7ve4qVtL/HWrrdo8jcxIGEAnx/zeebmz2V0+mjXFyO01jrhNnh/vOE4QKD1cZvXLZZhKcNIjE509XuFox7221u6TSDgDD3OmQrZE9yu5oQ0ev3sqWxdObjtEOHdFQ3srzk8tGanOL2rpw3rx6C0lp7WOAanxTMgJfbEtlaQXstay9qytSzevpjF2xazq3YXHuNh6oCpnJt3LobetZLvyPSRTBswjdTYVLdLkV7EralB1lq2HKxrty/shtIarHVG3Izsn8S8iQMpzE2jMDeV/H4Jrv9BfFR+H6z4KxScC2lD3K5GpFt4/V7e2/MeL257kTd2vEG9r56M2AwuH345F+RfwITMCT2q19IYQ6Rx4lQUmtvqJoXavmrrm1C+GWZ/x+1KOlRV72V7eR3byurZUdZyX8/28jr2Vze1O9cTYRiYGsug1DhmDe932CJMCq1yNNZaVh9czeLti3ll+yvsrt1NpIlkWvY0bplwC2cPPluhT+QEdNfUoLomHx/uqgxtq7NiRwWV9V7A2au7MDeNC8dnU5iXyqTBqSTFhtkfnBtegpq9cNFv3a5EpEv5Aj6W7VvGS9te4tXtr1LdXE1ydDIX5F/ABfkXUNy/GE+EpnnJ0SnU9lXLHoD4fjDmEld+vLWW0pomtpfVs72szrkvbw2wVQ3edudnJcWQlxHPrIJMctPjQz2tOenx9E+KIVKhVU5AwAb4+ODHLN7mBNm9dXuJjIhkRvYMvjDhC5ydezYpMSlulykiQdZadpTXh3pgl2+v4JN91aH1DwqyEjlvTP/QXNhhmYlERPTgXtjjsfxBSBoIw89zuxKRTteyINKLW19k8fbFlDeWEx8Zzzm55zA3fy4zsmdoizg5IQq1fVHlTlj/QnDhia6bO+fzB9hT2ci2srp2gXVHWT07yutp8PpD53oiDINS48jLiOfiCdkMyUggNyOevIx4ctPje948Jwk7ARvgowMf8fK2l3ll+yvsr99PVEQUpw08jS9P+jKzB89WkBXpYT7eVcUfXt/Iyh0VHKxtBiAh2sOk3FRuP6uAyXlpTB6cSmp8L1tEqXwrbHrNGU3lUfsnvUPLFJ8Xt77IS9teYn/9fmI8MZyZcyYX5F/ArEGziI2MdbtMOQnW7ydQU4O/pgZ/VTWB6ir81dX4q6rxV1eRds01eJK6dgtA/absi5Y/6NwX33DKH9Xo9bOjvP6wHtftZXXsrmhot19rTGQEuelOUJ01vF8osA7JSGBQWpyGCEunC9gAK0tX8sr2V3hl+yuU1pcSFRHFzEEz+Vrh15g9eLb2WRXpwSyWjftrOGNEZqgXdkT/JDzh3gt7LCv+CiYCCq91uxKRU7axYmMoyO6s2UlkRCSzBs7ijqI7mD14NglRCW6XKID1+fBXVxOorm4XSAOhx9UEalof+6urCLQcr62l3YI2h0g66yyFWulkviZnNcURcyE197jeUtXgDQXWHeX1bDvY0vNaz77qxnbnJsVGMiQjgXGDUrhofPse1/5JseE/HEx6PH/Az4rSFaFViw80HCA6IppZg2ZxXtF5nJlzplYVFOlCnblQ1IScVN6886xTLyqc+JphxSMw8gJIHuh2NSInZUf1jlCQ3VS5iQgTwbQB07hl/C2a4tOFAs3NbUJp1eEBtara6U1tE0j91dUEqqoI1Ncf9bNNTAye5GQiUpLxJKcQldWfiIICPMkpeJKT8aQkE5Gcgicl2TkvORlPSgqepCRMXFyXf3eF2r5m7TNQdwCm3Bw6ZK3lQE1TsIf18B7XloU3WmQmxZCXHs9pBRkMyUho1+OaGh/Vs1eTlF7JF/CxYv8KFm9fzKvbX6WssYwYTwynDzqd84acxxk5Z+hKsEg30R7yp+iTZ6H+IBSd+mgqke60t3YvL297mRe3vcjasrUAFGYVcte0uzg371z6xfVzucKez1qLbWrqcAhvqMe0prp9IK2uCvWe2sbGo36+iY93AmjwFjVoELGjRwcDqRNW2z1OTgqF04iYnr3dn0JtH+HzB9hb1Ujy23/GxOfyx/UD2PZeSWjocNv5rREGBqbGMSQjgQvHZ5OXHk9em/CaEKP/bMR9voCPkv0lLN62mNd2vEZ5YzmxnlhOzwkG2UFnEB8V73aZIiInpuRBZyTVsLPdrkTkmA42HGTxtsW8tO0lVpauBGBcxji+WfxNzh9yPgMSBrhcobustdj6enzl5fjLy/GVleOvKHeetzwuC75WUYG/vBzb1HTUz4xITAz2mDq9oNFDhnQcSFOS8SQltfaeJiVhonvZ+gNthHU6cWsvvJ6uqt7L8h3llGyrYM2earaX1bGrooHhdhsvxiznp97P8ciSHc781vR4ZgzLaB0mnB5PTlo80ZGa3yo9jzfgZdm+ZSzetpjXd7xORVMFcZFxnJlzJnPy5jBr0CwFWREJXwc2wLa34ZwfQoTaYemZqpqqeG3Ha7y49UWW7ltKwAYYnjacr07+KnOHzGVw8mC3S+wyoZAaDKC+sjL85RX4yp17f3l5a4AN3h8ppJq4OCLT0vBkZODJ7EfMiBF40tPxpKYGe1KTDh/Om5SEiQzr+NZlwvqfioY4tW5zULKtgpLtFSzfXs6G/bUAREYYhvdPYszAZC4Yn83lu5/EvyeWm2/7HndlDdD8VgkL3oCXD/Z+wCvbX+G1Ha9R1VRFfGQ8Zw4+k/PyzmPmoJnERXb9XA0RkS63/CGIiILJn3O7EpF26rx1vL7jdV7a9hLv7XkPX8BHblIut4y/hblD5lKQFr4dTIH6enzlFfjLy4JBtOVxBf6yMnwV5e2C65GG+JrYWCLT051gmpFOzPDheNLTicxIx5OWjic9jciMDDxp6USmpxERr4vwnSmsQ21f5PUHWLOnmpJt5Szf7gTZAzXOFaCWzeY/NXEgRXnpTBqcSlx0cLPqhkr47Qsw4UqyB2S79wVEjoPX72XJ3iUs3raYN3a+QXVzNQlRCcwePJvz8s7jtIGnadl/EeldvA2w6u8weh4kZrldjQiNvkbe3v02L259kbd2vUWTv4kBCQP4/OjPMzd/LqPTR/fIdVQCDQ3te0wPGfLrqwgO/Q2ec8SQGhODJyOdyLRgSB02DE9GBpHpaU5IzUgPhtgMhdQeQKG2h6tq8LJiRwXLt1VQsr2cVTsrafQGAMhJi2NWQT+K8tIoHpLGiKykI/e+fvgYeOvbLRAl0pM0+5tZsmcJi7cv5o0db1DjrSExKpGzBp/FeUPOY8bAGcR4evYiBSKiqUEnbe3T0FjZKdvtiZyslovKL259kdd3vE69r56M2AwuH345F+RfwITMCUSY7h0ab71eZ7hvWVlw/mkZvoNlzn1Z+WFDfm1DQ4efY2JinJ7TYG9qzLBhTq9qehqR6RnB+3QnuKalYeLje2Rol44p1PYg1lp2VTRQst2ZD1uyrYINpTVYC54Iw9iBycyfmktxXjrFQ9Lon3ycPVXWwrL7IWcKDJzUpd9B5EQ0+Zt4b/d7LN6+mDd3vkmtt5ak6CTOyj2L84ecz/Ts6UR7eu+iBiK9kaYGnaSShZBRAENOd7sS6WP8AT/L9i/jpa0v8cr2V6huriY5OpkL8i/ggvwLKO5fjCfC02k/z1pLoK4+GErLDgurvvIy/AeDQ4HLyvBXVnb4OSYqKthz6gTRmKH57XtQ04JDf4NBViG1d1OodZHXH2Dd3urgfFgnyJa2DCWOiWRyXhoXTcimOC+NiYNTT37V4S1vQtkm+PSCzite5CQ1+hp5d/e7LN6+mP/s+g913jqSo5OZkzeHOXlzmJ49nShPlNtlioh0n32rYecHcP5/g/7olm4QsAE+PPAhL259kcXbFlPWWEZCVAJnDz6buflzmZE944TaYuv346+sbN+DWnYQX1k5vrKDzrDf8nL8Bw8edchvRHKyM+80I52YggIip011hvf2y3DCab9+TmDt14+IhASFVAlRqO1G1Y1eVu6opGSbE2BX7awMbaUzKDWOGcMyKM5LoygvnZEDkvB01kJOy+6H+AwYe2nnfJ7ICWrwNfDO7nd4Zdsr/GfXf6j31ZMak8rcIXOZkzeHqdlTiYpQkBWRPmr5g+CJgYnz3a5EejFrLWvL1/LS1pd4adtL7KvbR4wnhjNzzuSC/AuYNWhWu/UqAo2NrSG1TVgNhdSWXtbycvwVFRAIHP5DIyNbh/RmZBCTPwRPRr9gD+ohYTUtrVdvOSNdS6G2i7QMJXYWc3JC7Pr9zlDiCANjBiZz9ZTBFA9JozgvnQEpXbToTdUuWP8CzPwaRGo+onSfem89b+9+m8XbFvP27rdp8DWQFpPGhUMv5Ly88ygeUKwgKyLSVAsfPg5jPw3x6W5X02c1+hrZWrWVTZWb2FS5ie3V2/EGvO3OsdZ2+F7L4cc7OtbxoSN8Zgc/61R+DsDemt2UH9hJeqOHs2PHMSNuNqN8A4j4sBb/629zoOyp4FBgJ6wG6us7/JyIhAQ8/TKITM8gKi+XuMJCJ6QGg2vbHtWIlBT1pkq3UKjtJD5/gHV7a5wAu91Z2GlftTO0IjEmksm5qVwwLpviIWlMOpWhxCdq+UPOnNriG7vn58lxqW2uZXPVZjZVdNxwhrt9dft4Z/c7NPobSY9NZ97QeZw35DyK+hcRGaFfO3Jk1u/HNjYSaGwM3bd93O6+oRHb1P6+321fILJfP7e/hsjxW/0kNNeone4m3oCXHdU72Fi5kU0VToDdXLmZHTU7CFinpzEyIpLcpNwOV9k3dBzQOjreUZjr8P1HyHxtzzV+S0yjn9jGADGN/va3pgAxDX5imvzO/aGvN/iJbvIT2+AnImABP7ACWEEFQERE6wJKGenETZwYDKn92oXVlh7XiFjtPiA9j/66PEk1LUOJg3vDrtxRSX2zM5R4YEosU/OdxZyK8tIYNSC584YSnwhfMyz/K4yYC6m53f/zhXpvPVurtrKxciObKzeHrgDvq9sXOic6IrrXreqbEJ3AJQWXcP6Q8ynMKuzUBSak+9lAwAmRTU3YhgYCjU3Yxo7vA40N2IZGAk2Nx7gPflboM52gar0nd4HHREVh4uJIu+ZqhVoJLyULIWssDJ7qdiW9ij/gZ3ft7lB43Vy5mY2VG9lWvQ1fwAdAhIkgNymXgtQCzh9yPgVpBQxPHU5ucm6njSSy1mLr6/HX1hKorSVQU4O/ppZAXS3+mhoCtXUEamqCz4Pn1NbgbzleW4u/tvaIK/q2ExmJJzGRiMREIpKS8KQlEJGURERiIp6kRCKSkjvsUfWkpGA8aqclvCnUHqfdlQ2hvWGXbatg/b5qAsGhxKOzk7myKIeiIekU56UxMDXO7XIda5+CulKYqm18ulqTv8kJrxVOeG1pPPfU7gkNDYqOiCY/JZ+i/kUUpBYwLGUYBWkFDEoc1O3L40vfYP1+/FVVwTlPbTaWLyvHX1lBoK6uNYh2eN+IbWjANjefXAGRkUTExmLiYomIiSUiLhYTE0tEbCye5GQi+mc5z1uOt72Pdc470v2hx/QHmYSl3Stg7yq48DdaIOokWWvZV7fPCa/BXteNFRvZUrWFJn9T6LxBiYMoSC3gjJwzKEgtYHjacPJT8o96UTnQ3NwaRGtrCRwaRmtr2gVTf13wnNpa/LUt59R2PNf0EBEJCcEwmognMQlPcjJRgwYGQ2oSEYkJeJKSiEhoOSexNbAGH5uYGA31lT5LobYDPn+AT/bVBAOsE2T3VjlDiROiPUzOTeMrZw+neEgak3PTSOyuocQnomo3vPRd5+rv0LPdrqbX8Pq9bKveFupxbel93Vmzs3XYkolkSMoQxvcbz6UFl1KQWkBBagE5STkaeiunxAYCTkht2a/vkKAa2lC+IrjZfGXlEf+Y8qSkOCtHxsURERODiYvDk5iE6ZcZDIoxRMTGHfH+mEE0JgYTpTnTfZH2qT0BJQshKgEmXO12JT2etZayxjI2VrQJr8FRUHXeutB5WXFZFKQVcNWAqxieOty5iJw6jPio+CN+dqCpiaYNG2hcu47GdWtpXLeOpo2bsEeYU9qWiY4OhssEPIlOyIzKHUxsQmLr8aSk9sE02JsaCqbx8bowJ3KK9Bd20M7yehYt38Xy7RWs3FFBXXAocXZKLEV5aRTnpVE8JJ1RA5KI9PTwXjVfM/zzOvA1wpUPQUQPr7cH8gV87KjZERqy1BJid1TvwGedYUse42Fw0mBGpI3ggvwLGJY6rNOHLUnvZq0lUFPjLMpRUeHctwTVDjaU91dUgN/f4WdFpKQQmZaGJyOD6CFDiCssar+hfEZG6559qamYSP36l66hfWqPU2OVM592/BUQm+x2NT1KZWNlqN1te6tqqgqdkxqTyvC04cwbOo/haa3hNSUm5aif7a+upnHdJzSuW0vTunU0rl1H05Ytod+tEYmJxI4aRepllxHZL+OYvaQRWq1XpEfQXzVB+6sb+cPrGxk1IJnLCnOcVYmHpDOopwwlPhGL74Jdy+DKv0LmCLer6dHazrlpG163VW0LLd5kMOQk5VCQWsA5uecwLHUYBakF5KfkE+1RYyatnA3l61q3OGgJo217T9sG1YoKOMIc0oikpFAojRo8mLiJE535TxntN5T3pKc72yCoV1QkvHz0BHjr+/QCUXXeunZDhlva4IMNB0PnJEYlUpBawJy8OaGRT8NSh5ERm3HUobbWWnylB9qF18Z16/Du2hU6JzIzk5gxo0k852xiR48hdsxoonJyNIRXJAz1uFBrjBkK3AWkWGuv6K6fO2lwKh/98DySYsP8D8OP/glLF8CM27UvbRsBG2Bv3d52Kx1uqtx02JybgQkDKUgrYNagWaHGMz8ln7jIMLy4IZ0i0NyM/8ABfMEN4zsKqr7yYA9rWdkRFzqKSEgIrS4ZlZ1N7LixRKY5K01Gpgf360t3elo9aWm6+i/Sm1nrDD0eONm59XKNvka2VG0JDRluGQW1p25P6JxYTyzDUocxc+BMp/1Nc9rg/vH9jxkybSCAd8cOGtuE18Z16/CXlYXOicrLJXbcOFKvvJLYMaOJHT1ai8qJ9CLdEmqNMQuBi4FSa+24NsfnAv8HeID7rbW/tNZuAW4yxizqjtpaRHoiSOrpw4qPpXQdPPtVyJ0B5/7I7WpcYa1lf/3+dsF1U8UmNldtpsHXunJgVnwWw1OHM2XAlFB4HZo6lISoBBerl+5kvV4nqJaW4jtwAG9pqfO49EDw3jnur6jo8P0mLi4YRtOJyswiduSoUC9qu6Aa7FGNiOldK1yLyCnY+QGUroVP3e12JZ2qZd2JduG1anO7dSeiIqLIT8lnUtYkrki9IhRgj3fRRNvcTNPmze3Ca9MnnxCoC86rjYwkpqCAxDPOIHb0aGLHjCZm1Cg8iYld+dVFxGXd1VP7EPBH4OGWA8YYD3APMAfYBSwzxjxjrV3bTTX1Lo3V8PjnIDrRmUfrCfMe52Nou2BE22HDmys3U+utDZ2XEZtBQVoBlw2/LDTndWjqUJKjNX+pt7I+n9N7WlqK70BpKKB6gyG1JbT6y8ud3pK2PB5nw/isLGfIb1EhUVlZRGZm4unXr3VeanoaEfFHXnREROSoShZCTDKMu9ztSk5awAbYXLmZlaUrWVG6gvXl69lWta3duhO5ybmMSBvBhfkXhsJrblLucS+aGKiro3H9+sMWcGqZtmHi44kdOZKUSy5xwuvo0cQMH66RLiJ9ULeEWmvtW8aYIYccngpsCvbMYoz5B3AJoFB7oqyFp78M5VvhumchaYDbFXW66uZqSvaVsGzfMtaWrWVz1ebDFowoSC3goqEXhXpeC1ILSI1Nda9o6VTW78dfXt7ao9omoIaC64FS/GXlh6/4GxHhbByflUVU//7EjR9PZFYWkVmZzrGW4JqerhUoRaRr1ZfDmqeg8FqIDp/RQc3+ZtaUrWHF/hWsLF3JytKVVDdXA5AZl8nYjLHMHjy73dSdE1l3wldeHgqvLXNgm7dvD1189KSlETt6NInXXUvM6NHEjh5DdF6ufmeLCODunNpBwM42z3cB04wxGcDPgcnGmO9aa3/R0ZuNMbcCtwLk5uZ2da0925J7YN0zMOenMGSm29V0inpvPStKV7B071KW7lvKuvJ1BGyAWE8so9JHnfCCEdJz2UAAf2Vl+3AaCq4HW48fPNjhyr+eYFiNzMokduwYIjNbw6rzOIvIjHSt9isiPcPS+8DfBMU3uF3JUdU017CqdBUrS1eyfP9yVh9cTXPA2bM6PyWfOXlzmJw1mcL+heQkHv/iStZavLv3HLaAk2///tA5UQMHEjNmNMnzLg4t4BTZ/9hza0Wk7+pxf+VZa8uA247jvAXAAoDi4mJ7jNN7r23vwis/gNHz4LSvuF3NSWvyN/Fh6Ycs3eeE2I8PfIzP+oiMiGRi5kS+MOELTB0wlQmZE7TicJiw1gbD6oE2PauloWHB3pb5qwcOgM932Ps9aWlEZjrhNGb48MN6VSOzsojs10+r/opI+PjkefjPL2H0p6D/WLeraWd/3f5QgF1ZupINFRuwWCJNJKMzRnPNqGso7F/I5KzJpMemH9dnWp+P5q1b2y/g9MknBKqCI60iIogemk/81Kmh+a+xo0bhSU3tui8qIr2Sm6F2NzC4zfOc4DE5XjX7YNENkDYELvkThNEVTG/Ay5qDa5wQu3cpK0tX0hxoJsJEMDZjLNeNvY6p2VOZnDVZKw+7KNDURKC6Gn9NLYHaGvw1NQRaHlfXBI/VEqipwV/rvOavqSZQWYXvwIEOVwKOSEkhKiuTyMwsYqbmB3tZs0K9rVFZWXgyMzUnSkTwlZXRsGoV1u8Hvx/rD0Cg7b0fAoHg6wFs4Aj3ft/RXw/et//s4L3Pf/jrfj82cBz3Pl/rc28TNFRiyYbo9Zj/m4mJjHRGkURFYiKjQs9NVNQRjkdCZPBY1BGOB5877285J6r1WGQkxhPJ/uaDbKjZwvqazayr3si+pgP4IiAqOobh/UZzYc7nGNd/IqOzxhEfn9Ja6xGG+wYaG2lav75dgG3asAHb5OwwYGJiiBkxguTzzw+tPhwzYgQRcWrjReTUuRlqlwHDjTH5OGH2GuAzJ/IBxph5wLyCgoIuKK+H83vhnzdAUw18/qkev3G7P+BnfcV6lu5dygf7PmDF/hXU++oBGJk2kqtHXc20AdMo7F9IUnSSy9X2DoHmZgI1NU7gbBdKDwmnLaG05bVgOA3U1Bxxe5oQY4hISHD2VE1KIiIpiajMLCIKCpwe1ba9qsHHEbGx3fMPQETCXuPatez68u2n/kEREeDxYE703uM54usmOrqD8yMgInjviWw97q2DDS9gIqNg3KUQFY/1ebE+H3h9WF/bmxe8XqzXh21sIuCrc87zBY/5fFivN3Q+bR8fup7AMeQHb3PbHa0Hlgdvztyww7QN1MGbr6ws9PMjkpOJHTWKtGuuaV3AaehQTQMRkS7TXVv6PAbMBvoZY3YBP7TWPmCMuR14GWdLn4XW2jUn8rnW2meBZ4uLi2/p7Jp7vFd/BDveg8vug/5j3K7mMNZaNldu5oN9H7B071JK9peEFpTIT8ln3rB5TMueRnH/YtJi01yutuexXi/+2kPCZrsg2qZX9NBwWltLoLoa29x8zJ8TER/vBNLkJCISk/BkpBOdl9fuWERSohNYE4P3Scl4khKJSEoiIiHB+aNNRCSoMy84x02ezJAnF2GOEjCJiHBeC923hMo257o5kqm+HB6YA/2a4abnod/wLvtRNhDA+nzU1Veyeu8qPt63itWlH7Kx9BP83iY8ARgcl82Y5BGMSh7OiOShZEVlgL8lLLcNzV6np9l7SOBued3besz6fERmZgaHEI8hatAgzX8VkW7VXasfzz/C8ReAF7qjhl5lzVOw5I8w5RaYcJXb1QBOiN1Zs5MP9n3Asr3L+GDfB5Q3lgMwKHEQ5+ady5QBU5g6YCpZ8VkuV+u+QEMDDR99TP3yEho/+hhfRXm7kGobG4/5GSYuLtQ76klMxJOSQlTOIDyJSUQkB3tOE5NaA2hiIp7k5NZjiYlaNVJEOl1nXnD2JCYSN7ZnzT09Id5GeGw+VO6Ea5/uskB7sOEgK/avYEXpClbsX8H6ivUEbIAIE8Go9FHMnnZVaD5sv7h+XVKDiIibwnocSJ8cfnxwIzx9OwwqhvP/29VS9tXtY+m+pXyw9wOW7lvKvrp9gLO0/4yBM5g2YBpTBkwhJynH1Tp7Al9FBQ0rV1K/fDkNJctpWLvW2WfPGKKHDSUqqz9RA7KDvaLJzn1iMLAmJTpBNDmpNZwmJmqBJBGRniwQgH9/AXa+D1c8CHkzOuVjrbVsr94eCrArS1eyo2YHALGeWCZkTuDWCbcyOWsyEzMnkhAVPtsGiYicrLAOtX1u+HFTLTz+OYiMhqv+6tx3o4MNBynZVxIaUtzSiKbFpFE8oJibx93M1OypDEke0ueHHXl376Z++XLql6+gfnkJzZs2Oy9ERRE3bhwZ119HXFER8ZMn40lJcbdYERHpfK/+ANY+Bef9DMZddtIf4wv4+KT8k1CAXVG6IjQSKi0mjclZk7lq5FVMzprM6IzRREXogqeI9D1hHWr7FGvh2a/BwQ3wuX9BStf3flY1VVGyvyS0V+ymyk0AJEYlUty/mGtGXcPUAVMZnjacCNN351XaQICmjZtoWLGc+pLl1K9YgW/vXgAiEhOJmzyZlIvnEV9USOz48d2yUJLX62XXrl00HscwZpG+LDY2lpycHKI08kE60wcL4L27YeqtMOPEFrqq99bz0cGPWLl/JctLl/PRgY9o8DUAkJOYw6xBs0L7w+Yn5/f5i8jhRG2zyPE5mbZZoTZcLL0PVi+Cs78Pw87qkh9R761n+f7lob1i15Wtw2KJ9cRS2L+Qi4dezLTsaYxKH0VkRN/9TyfQ3Ezj6tWhocT1K1cSqHYWwYrMzCSuuIj4m24ivqiQmBEjXJm3umvXLpKSkhgyRL3mIkdiraWsrIxdu3aRn5/vdjnSW3zyPLz0bRh5Icz95TG32ytrKGNV6SqWly5n5f6VrCtfh9/6MRhGpo/k0oJLnfmwmZPpn9C/m76EdAW1zSLHdrJtc1gnkz4zp3bnUnj5v2DEXJj1jU772EZfIx8e+DC0V+zqg6vxWR9REVFMzJzIFyd9kakDpjKh3wSiPH23F8NfU0PDqlVOL2xwYaeWlYWj8/NJPv884gqLiC8uIionp0c0VI2NjWo0RY7BGENGRgYHDhxwuxTpLXYth0U3wcDJcPkDEHH4Rc3q5mpe3/G6M5R4/wq2VW8DIDoimvGZ47lx3I0U9i9kYuZEbXHXy6htFjm2k22bwzrU9ok5tXUH4YnrIGUQfPpeZ6+9k+QNeFlzcE1oYadVpatoDjTjMR7G9hvL9eOuZ+qAqUzKmkRcZN/dDN27v7TdUOKm9eudBT88HmLHjCHtM58hrqiQ+KIiItPT3S73iNRoihyb/j+RTlO+FR69ChKzYP7jEB1/2ClrytZwxxt3sLduL8nRyUzOmsynh3+awqxCxmSMIdrTvWtlSPfT7xyRYzuZ/0/COtT2egE/LLoRGsrhplcg7sT2cw3YAOvK17F071I+2PcBK/avoMHXEBrSdM2oa5iWPY3CrEISoxO76Ev0bNZamrdubR1KvGIF3p07AWfLnLhJE+n3pS8RX1RI3MSJRMQf/keKHNn+/fu54447eP/990lLSyM6OppvfetbfPrTn+62Gt5//32+9rWv0dTURFNTE1dffTU/+tGPjnh+SUkJDz/8MH/4wx8Oe23IkCGUlJTQr9+Rt8Q4nnNEpJepL4e/XwHWD597EhIzDzvl3xv/zc/e/xnpcek8eP6DFPYv7NPrUYh71DZLb6RQ25O98XPY+h+45B7InnBCb6331vPN/3yTt3e/DcDQlKFcMuwSpmVPo7h/MamxqV1QcM9nvV4aP/mE+pLlTm/s8hX4y51VJD3p6cQXFZL2mc8QX1xE7KhR2jbnFFhrufTSS7nuuut49NFHAdi+fTvPPPNMt9Zx3XXX8cQTTzBx4kT8fj/r168/6vnFxcUUFxd3U3UiEva8DUfdi7bZ38wvlv6CRRsWMS17Gr8641ekx/bcUT7Su6ltlt4qrC8RGmPmGWMWVFVVuV1K51v/Irz9v1B4LUz+3Am9tbyxnJtevol397zL/yv6f7x+5es8fenT3DX9Ls7NO7dPBdpAXR11S5Zw4O4/sv2GG1g/bTrbrryK0v/5Hxo/WU/iGWcw4Kc/YegLLzD83XfIuftuMm64nrjx4xVoT9Hrr79OdHQ0t912W+hYXl4eX/nKV9i2bRunn346hYWFFBYW8t577wHw5ptvcvHFF4fOv/3223nooYcA+M53vsOYMWOYMGEC3/zmNwH45z//ybhx45g4cSJnnHFGh3WUlpaSnZ0NgMfjYcyYMQDU1dVx4403MnXqVCZPnszTTz99WA1lZWWcd955jB07lptvvhlrbehz//a3vzF16lQmTZrEF77wBfx+f2f8YxORcNJ2L9rL/nLYXrT76vZx/UvXs2jDIm4cdyP3nnuvAq24Sm2z9FZh3VPba+fUlm+Bf30BsifCBb8+obfurt3NF175Avvq9vH72b/nrNyuWSm5p/KVl7cbSty4di34/WAMMaNGkXrZZc5Q4sIiovpnuV1ur7ZmzRoKCws7fC0rK4tXXnmF2NhYNm7cyPz58ykpKTniZ5WVlfHvf/+bTz75BGMMlZWVAPzkJz/h5ZdfZtCgQaFjh7rjjjsYOXIks2fPZu7cuVx33XXExsby85//nLPPPpuFCxdSWVnJ1KlTOffcc9u998c//jGzZs3iBz/4Ac8//zwPPPAAAOvWrePxxx/n3XffJSoqii996Uv8/e9/59prrz3xf1AiEr5e+T6sfdrZi3Zs+6GbS/cu5c637qTR18jvZv+Oc/POPcKHiHQftc3SW4V1qO2VvA3w+LXOFgBXPQxRx7+n6fry9dz26m00+5u577z7mJw1uQsLdZ+1Fu+uXa1DiUuW07x1KwAmOpq4CRPIuPlm4ouLiJs0CU9S31xF8sfPrmHtnupO/cwxA5P54byxJ/SeL3/5y7zzzjtER0fz6quvcvvtt7Nq1So8Hg8bNmw46ntTUlKIjY3lpptu4uKLLw5drZ05cybXX389V111FZdddlmH7/3BD37AZz/7WRYvXsyjjz7KY489xptvvsnixYt55pln+M1vfgM4q1Lu2LGj3Xvfeust/vWvfwFw0UUXkZbmzGt/7bXXWL58OVOmTAGgoaGBrCxdJBHpUz74Cyz542F70Vpr+euav/K7Fb8jLzmPB89/kKGpQ10sVHoitc1qm6VzKdQG+SsrafzkE3eLsBbe/QNs3gBn/wDW7wP2Hddb15ev5+6VdzM6Mo47iv4fA7c0Ubfl/a6t1w2BAE1btlK/vISG5SvwlZYCEJGcTHxhISmXfZr4omJix40lIlqrSLpp7NixPPnkk6Hn99xzDwcPHqS4uJjf/e539O/fnw8//JBAIEBsrHPxJjIykkAgEHpPywb1kZGRLF26lNdee41Fixbxxz/+kddff517772XDz74gOeff56ioiKWL1/ON7/5TVauXMnAgQN54YUXABg2bBhf/OIXueWWW8jMzKSsrAxrLU8++SQjR45sV/f+/fuP+d2stVx33XX84he/OOV/TiIShj55Hl48fC/aOm8dP3j3Byzevpg5eXP46cyfkhCV4HKxIq3UNktvFdahtjP3qW385BN2XH/DqRfVKTLgjXtO6B1xwLcAqMP34F3sOPrpYS9ywADip0xxemELi4gZXoA5he2OerMTvWrbWc4++2z+67/+iz//+c988YtfBKC+vh6AqqoqcnJyiIiI4K9//WtozkteXh5r166lqamJhoYGXnvtNWbNmkVtbS319fVceOGFzJw5k6FDnV6PzZs3M23aNKZNm8aLL77Izp07efDBB9vV8fzzz3PhhRdijGHjxo14PB5SU1M5//zzufvuu7n77rsxxrBy5UomT24/uuGMM87g0Ucf5Xvf+x4vvvgiFRUVAJxzzjlccskl3HHHHWRlZVFeXk5NTQ15eXld+s9URHqAXSXOXrSDCtvtRbu1aitff+PrbKvexh1Fd3DD2Bu0fYsckdpmtc3SucI61HbmnNrYMWPIe+ThTqjqJB3cAM9/01nleM5P4DiX+X91x2v8fe3fGJ42nK8Wfo3EPnBFOGrgQKIGDXK7DDkGYwxPPfUUd9xxB7/61a/IzMwkISGB//mf/6GwsJDLL7+chx9+mLlz55KQ4Px3O3jwYK666irGjRtHfn5+qCGrqanhkksuobGxEWstv/3tbwG488472bhxI9ZazjnnHCZOnHhYHY888gh33HEH8fHxREZG8ve//x2Px8P3v/99vv71rzNhwgQCgQD5+fk899xz7d77wx/+kPnz5zN27FhOO+00cnNzARgzZgw/+9nPOO+88wgEAkRFRXHPPfeo4ZRezRiTAPwH+JG19rljnd8rlW+BR6+GpP7t9qJ9bftr3PXuXURHRPOXOX9hevZ0lwsV6ZjaZumtTNsVw8JVcXGxPdpE9h6vvhz+ciZg4db/QELGMd9ireXulXdz38f3cdbgs/jVGb8iNvL4599K77Zu3TpGjx7tdhkiYaGj/1+MMcuttb1i/whjzELgYqDUWjuuzfG5wP8BHuB+a+0vj/E5PwFqgbXHE2rDvm0+VF0ZPDAnuHf8q9CvAH/Az90r7+aB1Q8wLmMcv539W7ITs92uVHootc0ix+9E2+aw7qntFQIB+NctULsPbnzpuAKtL+DjJ0t+wr83/ZsrRlzBXdPuIjJC/ypFRKRDDwF/BELDkYwxHuAeYA6wC1hmjHkGJ+AeOiHtRmAisBbom1dPvQ3wj/lQtQuuewb6FVDRWMG33voW7+99n8uHX853p32XGE+M25WKiPRJSkJue+vXsOlVuOi3MKjomKc3+Bq48z938p9d/+G2ibfxpYlf0pwdERE5ImvtW8aYIYccngpsstZuATDG/AO4xFr7C5xe3XaMMbOBBGAM0GCMecFaGzj0vF4ptBftUrjyIcidzpqDa7jjzTsoayjjx6f9mMuGd7zCq4iIdA+FWjdtfBXe/AVMnA/FNx7z9MrGSm5//XY+OvAR35/+fa4aeVU3FCkiIr3QIGBnm+e7gGlHOtlaexeAMeZ64OCRAq0x5lbgViA0zy3shfai/TmMvZR/b/w3P3v/Z2TEZfDwBQ8ztp87C/6IiEirsA61nbn6cber3AH/uhmyxji9tMfobd1bu5cvvPoFdtfs5rezf6tN3EVEpNtZax86xusLgAXgzKntjpq6VGgv2i/QPPUWfrHkxyzasIjp2dP51Rm/Ii02ze0KRUQECOs9UKy1z1prb01JSXG7lBPja4InroWAH65+JLR64pFsqNjA5174HAfrD/KXOX9RoBURkVO1Gxjc5nlO8Ji0WPdccC/ai9h3+te4/uUbWLRhETeNu4l7z71XgVZEpAcJ657asPXit2HPSrjmUcgYdtRTS/aV8NXXv0pcZBwPXfAQI9JGdFORIiLSiy0Dhhtj8nHC7DXAZzrjg8N6FFWLncvgSWcv2g9m3sa3XvgMTf4mfj/795yTd47b1YmIyCHCuqc2LK16FJY/CDO/DqMuOuqpr21/jS+88gX6xffjkQsfUaCVsPLzn/+csWPHMmHCBCZNmsQHH3wAwO9///vQRu9uqays5E9/+tMRX/d4PEyaNImxY8cyceJE/vd//5dAwN01cf77v//b1Z8v4csY8xiwBBhpjNlljLnJWusDbgdeBtYBT1hr13TGzwvbUVQtyrfAY1djk/rzYOGnufXNr5Aak8qjFz2qQCthT21z51Lb3HMo1HanfR/Dc3fAkNPh7O8f9dQn1j/BN/7zDUZljOLhuQ8zMHFgNxUpcuqWLFnCc889x4oVK/joo4949dVXGTzYGekYDg1nXFwcq1atYs2aNbzyyiu8+OKL/PjHP+7GCg+nhlNOlrV2vrU221obZa3NsdY+EDz+grV2hLV2mLX2527X2SPUlcHfrqCOAP9v1HR+u3oB5+Sew6MXPcrQlKFuVydyStQ2dz61zT2HQm13aaiExz8PcWlwxULwdDzy21rLn1b9iZ++/1NOH3Q69593P6mxqd1aqsip2rt3L/369SMmxtmzsV+/fgwcOJA//OEP7Nmzh7POOouzzjoLgMWLFzNjxgwKCwu58sorqa2tBWD58uWceeaZFBUVcf7557N3714AZs+ezde+9jUmTZrEuHHjWLp0KQB1dXXceOONTJ06lcmTJ/P0008DsGbNGqZOncqkSZOYMGECGzdu5Dvf+Q6bN29m0qRJ3HnnnUf9LllZWSxYsIA//vGPWGvx+/3ceeedTJkyhQkTJvCXv/wl9J3POOOMUF1vv/02AC+99BKFhYVMnDiRc84556i1PvTQQ1x22WXMnTuX4cOH861vfQuA73znOzQ0NDBp0iQ++9nPds6/JJEuYoyZZ4xZUFVV5XYpJya4F+2Wur18Jn8Er+3/gP9X9P/43zP/l4SoBLerEzllapvVNvdq1tqwvxUVFdkeze+39tFrrP1xurXblxzxNJ/fZ3/03o/suIfG2e+/833r9Xu7sUjpTdauXevqz6+pqbETJ060w4cPt1/84hftm2++GXotLy/PHjhwwFpr7YEDB+zpp59ua2trrbXW/vKXv7Q//vGPbXNzs50xY4YtLS211lr7j3/8w95www3WWmvPPPNMe/PNN1trrf3Pf/5jx44da6219rvf/a595JFHrLXWVlRU2OHDh9va2lp7++2327/97W/WWmubmppsfX293bp1a+h9HUlISDjsWEpKit23b5/9y1/+Yn/6059aa61tbGy0RUVFdsuWLfY3v/mN/dnPfmattdbn89nq6mpbWlpqc3Jy7JYtW6y11paVlR211gcffNDm5+fbyspK29DQYHNzc+2OHTuOWJN0jo7+fwFKbA9o38L51uPb5rb8fmv/8Tn7yv8MsNMeLrKnP3a6fX/P+25XJb2M2ma1zXL8TrRt1kJR3eHd38P6F2DuLyF3eoenNPoa+fZb3+b1na9zy/hb+Mrkr2COsc2PyHF58TvO0PfONGA8XPDLI76cmJjI8uXLefvtt3njjTe4+uqr+eUvf8n111/f7rz333+ftWvXMnPmTACam5uZMWMG69evZ/Xq1cyZMwcAv99PdnZ26H3z588H4IwzzqC6uprKykoWL17MM888w29+8xsAGhsb2bFjBzNmzODnP/85u3bt4rLLLmP48OGn9NUXL17MRx99xKJFiwCoqqpi48aNTJkyhRtvvBGv18ull17KpEmTePPNNznjjDPIz88HID09PfQZHdUKcM4559AyF3HMmDFs3749NDxMRLqGf/H3uHvff3igfz/Gp4/gt7N/y4CEAW6XJb2Z2ma1zdKpwjrUhsUKi1v+A6//FMZeBtNu6/CUqqYqvvr6V1lZupLvTv0unxndKQtQirjK4/Ewe/ZsZs+ezfjx4/nrX/96WMNprWXOnDk89thj7Y5//PHHjB07liVLlnT42Yde8DHGYK3lySefZOTIke1eGz16NNOmTeP555/nwgsv5C9/+QtDh57Y3LgtW7bg8XjIysrCWsvdd9/N+eeff9h5b731Fs8//zzXX3893/jGN0hL63jLjyPV+sEHH4SGhYHzz9Dn851QrSJyYire+R3f2vYE76emcMXwK/jutO8S7Yl2uyyRLqG2WW1zbxXWodZa+yzwbHFx8S1u19Kh6j2w6EbIGA6fuhs66HndV7ePL776RbZXb+dXZ/6KuUPmulCo9GpHuWrbVdavX09EREToyuuqVavIy8sDICkpiZqaGvr168f06dP58pe/zKZNmygoKKCuro7du3czcuRIDhw4wJIlS5gxYwZer5cNGzYwduxYAB5//HHOOuss3nnnHVJSUkhJSeH888/n7rvv5u6778YYw8qVK5k8eTJbtmxh6NChfPWrX2XHjh189NFHTJw4kZqamuP6LgcOHOC2227j9ttvxxjD+eefz5///GfOPvtsoqKi2LBhA4MGDeLgwYPk5ORwyy230NTUxIoVK7jrrrv40pe+xNatW8nPz6e8vJz09PQj1no0UVFReL1eoqKiTuHfjEjXC4sLzkFrlt3LHesXUBYXz09m/JBPj7jc7ZKkr1DbrLZZOlVYh9oezdcMT1wHvka4+m8Qk3jYKZsrN3Pbq7dR01zDvefey9TsqS4UKtL5amtr+cpXvkJlZSWRkZEUFBSwYMECAG699Vbmzp3LwIEDeeONN3jooYeYP38+TU1NAPzsZz9jxIgRLFq0iK9+9atUVVXh8/n4+te/Hmo4Y2NjmTx5Ml6vl4ULFwLw/e9/n69//etMmDCBQCBAfn4+zz33HE888QSPPPIIUVFRDBgwgP/6r/8iPT2dmTNnMm7cOC644AJ+/etft6u/ZeEHr9dLZGQkn//85/nGN74BwM0338y2bdsoLCzEWktmZiZPPfUUb775Jr/+9a+JiooiMTGRhx9+mMzMTBYsWMBll11GIBAgKyuLV1555Yi1Hs2tt97KhAkTKCws5O9//3un/vsS6Uw9/oJz0L8++F9+vu5BMjzRPHzeQsZmF7ldkkiXUtustrk3M86c2/BWXFxsS0pK3C6jvRe/DR/cC1c8COMuO+zlVaWr+PJrXybaE82fz/0zo9JHuVCk9Fbr1q1j9OjRbpfRJWbPns1vfvMbiouL3S5FeomO/n8xxiy31uo/slPQI9tmoNnfzH+/9V2e3LGY6V741WVPkZY+zO2ypA9Q2yxy/E60bVZPbVf4eJETaKd/qcNA+8aON7jzrTsZkDCAe8+9l5ykHBeKFBER6Vv21e3jjte+wuqKT7i5tpnbr3oGjwKtiEjYU6jtbKXr4JmvwODpMOcnh7385IYn+cn7P2FM+hjuOfce0mPTXShSJHy9+eabbpcgIkfRU+fUfrD3A+78zzdpbqzi9wcrOeeqRZA58thvFJFjUtssbotwu4BepakGHv88RCfAlQ+Bp3XSuLWWv3z4F3605EfMGDiDB85/QIFWRER6HWvts9baW1u2v3CbtZaFqxdy6yu3kt7cwGO793DOhX+E3GlulyYiIp1EPbWdxVp4+stQvgWufRqSW/ft8gf8/GLpL3h8/eN8atin+NFpPyIqQqukiYiIdKU6bx3ff/f7vLL9Fc6Lyean65cSf97PYcwlbpcmIiKdSKG2s7z/J1j7tDPkOP/00OEmfxPfffu7vLL9FW4YdwN3FN5x2D5eIiIi0rm2VG3h6298nR3VO/hm5mlcu/QfmGm3OetdiIhIr6JQ2xm2vweLvw+jLobTvho6XN1czdde/xol+0v41pRv8fkxn3exSBERkb7h1e2vctc7dxEbGcuCkdcz9cUfOW30+f/d4Z7xIiIS3jSn9lTV7Id/3gBpQ+DSP4Uay9L6Uq5/6XpWHVjFr874lQKt9Dn79u3jmmuuYdiwYRQVFXHhhReyYcOGk/qsm2++mbVr13ZqfZWVlfzpT3864ffV1tbyxS9+kWHDhlFYWEhRURH33XffKdXy0EMPcfvttwNw77338vDDD5/U52zbto1HH330lGoROVXGmHnGmAVVVVXd/rN9AR+/W/477njzDoalDuPxwv9i6iu/gEFFcNl9EOHp9ppEehK1zcdPbXN4CetQ62bDCYDfB4tugMYquPoRiHUWxdhatZXPv/B5dtfs5k/n/IkL8i9wpz4Rl1hr+fSnP83s2bPZvHkzy5cv5xe/+AX79+8/qc+7//77GTNmTKfWeLIN580330xaWhobN25kxYoVvPTSS5SXlx92ns/nO6m6brvtNq699tqTeq8aTukJ3FooqryxnNtevY2Fqxdy5YgreWjKDxjw7y9CUjZ85nGIju/WekR6GrXNapt7s7AOta6vsPjaj2H7uzDv/6D/WAA+OvAR1754LY3+Rh6c+yAzBs5wpzYRF73xxhtERUVx2223hY5NnDiR008/HWstd955J+PGjWP8+PE8/vjjgLMdwOzZs7niiisYNWoUn/3sZ7HWAs6m7iUlJQAkJiaGPnPRokVcf/31AGzevJnp06czfvx4vve974XOq62t5ZxzzqGwsJDx48fz9NNPA/Cd73yHzZs3M2nSJO68804Afv3rXzNlyhQmTJjAD3/4w8O+1+bNm1m6dCk/+9nPiIhwfn1mZmby7W9/O/QdTj/9dD71qU+FGvpLL72UoqIixo4dy4IFC0Kf9eCDDzJixAimTp3Ku+++Gzr+ox/9iN/85jehnzd37lyKioo4/fTT+eSTTwC4/vrr+epXv8ppp53G0KFDWbRoUeg7vf3220yaNInf/e53J/hvTeT/t3fn4VHV59/H39+kkU0iIGuLCoYiJCQBE2IAAxENUEF+oFAEfooPgrjg0qdQWewD/QmC6IOIdUOWgkWxLLFoW0UsGBAQiApFQAIWiiySBJuK7Mn9+yNxZMkyCZNMZvJ5XddcV+bMOd9z3+fMOXe+Z84SuLZlbWPAuwP47JvP+J+O/8P/i3mIyxYNzL+J438vhVr1/R2iiN+pNqs2BzNdU1tW25fDupkQfy/EDgAg7es0Rn00ivo16vPqLa9yVfhVfg5SxD+2bdtGXFxcoZ8tW7aMzz//nC1btpCVlUX79u3p3LkzAJ999hlffPEFP/3pT+nUqRMff/wxN954o1fzfPTRR3n00UcZOHAgr7zyimd49erVSU1NJTw8nKysLBITE+nduzdTp05l27ZtfP755wCsWLGCjIwMNm7ciJnRu3dv0tLSPLEBfPHFF8TGxnqKZmE+/fRTtm3bRvPmzQGYO3cu9erV48SJE7Rv35477riD06dPM2HCBNLT07niiiu46aabaNeu3UVt3Xfffbzyyiv8/Oc/55NPPuHBBx/k73//OwCHDh1i7dq17Ny5k969e9OvXz+mTp3Ks88+y7vvvuvVMhMJBkt3LWXyJ5NpUKMBC25dQFT4tTC/N/znIAx5B66M8HeIIpWCarNqczBTp7YssnbD2w/mX6PTYwoAb+9+m4nrJtKybkteuuUl6tfQUWGpHJ7e+DQ7j+70aZut6rXi8YTHyzTt2rVrGThwIKGhoTRq1IguXbqwadMmwsPDSUhIoGnTpgC0bduWvXv3el04169fz9tvvw3AoEGDGDVqFJB/utW4ceNIS0sjJCSEAwcOFHqq1YoVK1ixYoWngB07doyMjIzzCueFJk+ezOLFizly5AgHDx4EICEhwVM0AWbOnElqaioA+/fvJyMjg8OHD5OcnEyDBg0AGDBgwEXXNB07dox169bRv39/z7BTp055/u7Tpw8hISFERkaW+dQxkUB2KvcUUz6ZwtKMpXRo0oGnOz9N3cvCYfEQ+HoT/HI+XJXg7zBFCqXarNosvqVObWmd/h7+dBeEhkH/+VjoZcz5x2ye//R5EpskMuOmGdQKq+XvKEX8KioqynPaTWlUq1bN83doaGih176c+0iskydPltjmwoULyczMJD09nbCwMJo1a1bodGbG2LFjGTFiRJFtRUZGsmXLFvLy8ggJCWH8+PGMHz/+vNOuatX6cftfvXo1K1euZP369dSsWZPk5GSvYgbIy8ujTp06nqPVFzp3Wf1wKphIVXHo2CF+tfpXfJH9BcOjh/NQ24cIDQmF98bBjneg+xQ9i1bkAqrN+VSbg5M6taVhBu88Ckd2wF3LyLviZ0zb9DQLdyzkF81/weROkwkLDfN3lCLnKetR20vRtWtXxo0bx6xZs7jvvvsA2Lp1Kzk5OSQlJfHqq68yZMgQjh49SlpaGs8884znmpSSNGrUiB07dnDdddeRmppK7dq1AUhMTGTp0qUMGDCARYsWecbPycmhYcOGhIWFsWrVKvbt2wdA7dq1+e677zzjde/end/+9rcMHjyYyy+/nAMHDhAWFkbDhg0947Ro0YL4+HieeOIJnnzySUJDQzl58mSRhSsnJ4e6detSs2ZNdu7cyYYNGwC44YYbePTRR8nOziY8PJzFixcTGxt73rTh4eE0b96cxYsX079/f8yMrVu3XjTeuS7MScQfnHO3Abe1aNGiXNrfcGgDv/noN5zOO82Mm2Zw89U3F3zwMmx4EW54ADroWbRSuak2qzaLbwX0jaIq3KbZ8I/FcNN4Tje7kd+k/YaFOxZyV+RdTE2aqg6tSAHnHKmpqaxcuZKIiAiioqIYO3YsjRs3pm/fvsTExBAbG0vXrl2ZNm0ajRs39qpNgKlTp9KrVy86duxIkyZNPJ/PmDGD6dOnExMTw+7du/nhBnKDBw9m8+bNREdHs2DBAlq1agXAlVdeSadOnWjTpg2jR4+mW7duDBo0iA4dOhAdHU2/fv0KLUKzZ88mOzvbU0RTUlKYNm1aoTH36NGDs2fP0rp1a8aMGUNiYiIATZo0YeLEiXTo0IFOnTrRunXrQqdfuHAhc+bMITY2lqioKM+NNIoSExNDaGgosbGxuhmF+E153cTRzJjzjzmM+GAE9arX482eb/7Yod3xDrw3tuBZtJN9Ol+RYKHanE+1OTi5YPhpPD4+3n64+1q52b8J5v0CIrpy7I7ZPPbR/+WTw5/w67hfc0+be8p33iKltGPHjiJ3xoEoOjqa5cuXn3c9zIWOHz9OjRo1cM6xaNEi3nzzzRILjQgUvr0459LNLN5PIQUFX9bm7898z28//i0f7PuAbtd048lOT1IzrOARPfs3wvzboHE03L1cj+6RSku1WbVZvFfa2qzTj73xfVb+jSfCm5D1i6d44IN72f3tbp668Slui7jN39GJBLWUlBSio6OLLZoA6enpjBw5EjOjTp06zJ07t4IiFJHy9FXOVzy26jH+9Z9/MSp+FHdH3v3j9XvZe+CNAfnPoh24SB1akQqi2iyVjTq1JcnLhaX3wvdZ7LtzASNWjeToyaO8cPML3Pgz7+78JiJl98EHH3g1XlJSElu2bCnnaESkIm08tJGH//4w1X9Snde6vUb7xu1//PD7LFjYL/9vPYtWpEKpNktlo05tSVY9BV+tZtst43gw/SkA5nafS5v6bfwcmIiISHC7ts61JDZJZOwNY2lc65zr+86cgDfv1LNoRUQEUKe2eF++B2ue5eM2PfnVvqXUq16PV1Ne5Zrwa/wdmYiISNCrX6M+z3d9/vyBebmwdBh8vVnPohUREUB3Py7a0X9C6n2887NWjDy+natrX80fb/2jOrQiIiL+tOIJ2PkudH9Kz6IVERFAndrCnTkBf7qL+bUuY9xlx7m+URzzesyjfg1dryMiIlIc59xtzrlZOTk5vm98/Uuw4SU9i1ZERM5T6Tq1zrlazrn5zrnXnHOD/RFD3l9+zTOn9/NseA26N+vOy7e8TO3LavsjFJGAdfjwYe68804iIiKIi4vj1ltvZdeuXWVqa9iwYWzfvt2n8f373//mpZdeKvV0x44d44EHHiAiIoLrr7+euLg4XnvttUuK5Q9/+AMjR44E4JVXXmHBggVlamfv3r288cYbRX5+6NAhevXqVaa2vdWsWTOysrIA6NixY6Fxbd68mUceeaRc4zh3mf7+97/XHTcrUHk9p5bty+H9cdD6Nj2LVqSMVJu9p9rse+VZmyukU+ucm+ucO+Kc23bB8B7OuS+dc7udc2MKBt8OLDGz4UDviojvXGc2zWXcgfdYcEU4g1oNYlrnaVwWellFhyES0MyMvn37kpyczJ49e0hPT2fKlCl88803ZWpv9uzZREZG+jTGshbOYcOGUbduXTIyMvj000957733OHr06EXjnT17tkxx3X///dx9991lmrakwjl9+nSGDx9eprbLYt26dcDFccXHxzNz5swKi2Po0KG88MILFTY/KQf7N8Ky4dA0Hm5/DUJC/R2RSMBRbVZthuCtzRX1S+0fgB7nDnDOhQIvAr8AIoGBzrlIoCmwv2C03AqKD4Dv/7WOkelT+cvltXi07cOMSRhDiKt0P2aLVHqrVq0iLCyM+++/3zMsNjaWpKQkzIzRo0fTpk0boqOjeeuttwBYvXo1ycnJ9OvXj1atWjF48GDMDIDk5GQ2b94MwOWXX+5pc8mSJdxzzz0A7Nmzh8TERKKjo3niiSc84x07doybb76Z66+/nujoaM9D38eMGcOePXto27Yto0ePBuCZZ56hffv2xMTEMGHChIvy2rNnDxs3bmTSpEmEhOTvGxo0aMDjjz/uySEpKYnevXt7Cn2fPn2Ii4sjKiqKWbNmedqaN28eLVu2JCEhgY8//tgzfOLEiTz77LOe+fXo0YO4uDiSkpLYuXMnAPfccw+PPPIIHTt25Nprr2XJkiWenNasWUPbtm157rnnLop/6dKl9OiRvyvOzc1l1KhRtGnThpiYGE9h+fDDD2nXrh3R0dEMHTqUU6dOAflHeSdMmOBZjj/Ekp2dTbdu3YiKimLYsGGedXbuurowrtWrV3uOSh89epQ+ffoQExNDYmIiW7du9SyHoUOHkpyczLXXXnteof3jH/9IQkICbdu2ZcSIEeTm5ha7TGvWrEmzZs3YuHHjRctEAsAPz6IN/2n+s2jDavg7IpGApNqs2nzuugq62mxmFfICmgHbznnfAXj/nPdjC153Ab0Khi3ypu24uDi7VFnZu+2Xs9tY7LwoS932+iW3J+JP27dv9+v8n3/+eXvssccK/WzJkiV2yy232NmzZ+3w4cN21VVX2cGDB23VqlUWHh5u+/fvt9zcXEtMTLQ1a9aYmVmXLl1s06ZNZmZWq1YtT1uLFy+2IUOGmJlZz5497Y033jAzs5dfftkz3pkzZywnJ8fMzDIzMy0iIsLy8vLsn//8p0VFRXnaev/992348OGWl5dnubm51rNnT/voo4/Oi/3Pf/6z9enTp8i8V61aZTVr1rSvvvrKMyw7O9vMzI4fP25RUVGWlZVlBw8etKuuusqOHDlip06dso4dO9pDDz1kZmYTJkywZ555xszMunbtart27TIzsw0bNthNN91kZmZDhgyxfv36WW5urn3xxRcWERHhmX/Pnj0Lje2rr76y66+/3vP+pZdesjvuuMPOnDnjifPEiRPWtGlT+/LLL83M7K677rLnnnvOzMyuueYamzlzppmZvfjii3bvvfeamdnDDz9sv/vd78zM7N133zXAMjMzzezHdXVhXOe+HzlypE2cONHMzD788EOLjY31LIcOHTrYyZMnLTMz0+rVq2enT5+27du3W69evez06dNmZvbAAw/Y/Pnzi12mZmaTJk2yZ599ttBlU9j2Amy2CqqPwfryRW22Y5lmz7c1e7q5WdbuS29PxI9Um1WbL6Ta7Lva7M9H+vyMH3+RBfgauAGYCfzeOdcTeKeoiZ1z9wH3AVx99dWXHMw/961if4gxs+2v6Bz135fcnkhlcfippzi1Y6dP26zWuhWNx40r07Rr165l4MCBhIaG0qhRI7p06cKmTZsIDw8nISGBpk2bAtC2bVv27t3LjTfe6FW769ev5+233wZg0KBBjBo1Csg/cDdu3DjS0tIICQnhwIEDhZ5qtWLFClasWEG7du2A/KPIGRkZdO7cuch5Tp48mcWLF3PkyBEOHjwIQEJCAs2bN/eMM3PmTFJTUwHYv38/GRkZHD58mOTkZBo0aADAgAEDLrqm6dixY6xbt47+/ft7hv1wZBbyjzKHhIQQGRnp1aljhw4d8swPYOXKldx///385Cf5ZaBevXps2bKF5s2b07JlSwCGDBnCiy++yGOPPQbA7bffDkBcXBzLli0DIC0tzfN3z549qVu3bomxnGvt2rUsXboUgK5du5Kdnc1//vMfT3vVqlWjWrVqNGzYkG+++YYPP/yQ9PR02rdvD8CJEydo2LAhn3zySbHLtGHDhp4j2BJADm+F49kweImeRStBRbVZtRlUm31Zmyvdc2rN7Hvg/3gx3ixgFkB8fLyVMHqJ4tsN4/2IW6kd/tNLbUqkyouKivKcdlMa1apV8/wdGhpa6LUvzjnP3ydPniyxzYULF5KZmUl6ejphYWE0a9as0OnMjLFjxzJixIgi24qMjGTLli3k5eUREhLC+PHjGT9+/HmnXdWqVcvz9+rVq1m5ciXr16+nZs2aJCcnexUzQF5eHnXq1OHzzz8v9PNzl1X+wcvi1ahRw+t5F+WHeRa1bnytsO+DmTFkyBCmTJly3rg//NNUlJMnT1Kjhk5bDTgRXeGxf0B1H990SqQKUm3Op9p8aSprbfbnBaMHgKvOed+0YJjfqEMrwajxuHFc8/oCn75KOhLctWtXTp06dd51Klu3bmXNmjUkJSXx1ltvkZubS2ZmJmlpaSQkJHidT6NGjdixYwd5eXmeo6wAiYmJnqOKixYt8gzPycmhYcOGhIWFsWrVKvbt2wdA7dq1+e677zzjde/enblz53Ls2DEADhw4wJEjR86bd4sWLYiPj+eJJ57wXCty8uTJIgtXTk4OdevWpWbNmuzcuZMNGzYAcMMNN/DRRx+RnZ3NmTNnWLx48UXThoeH07x5c89nZsaWLVuKXTYX5nSuli1bsnfvXs/7lJQUXn31VU8BPHr0KNdddx179+5l9+7dALz++ut06dKl2Hl27tzZc6OJv/3tb3z77beliispKYmFCxcC+f9o1K9fn/Dw8CLnd/PNN7NkyRLPujl69Cj79u0rcZnu2rWLNm3aFJuLVFLq0EoQUm1WbQbVZl/WZn92ajcBP3fONXfOXQbcCSwvTQPl+iw8ESkz5xypqamsXLmSiIgIoqKiGDt2LI0bN6Zv377ExMQQGxtL165dmTZtGo0bN/aqTYCpU6fSq1cvOnbsSJMmTTyfz5gxg+nTpxMTE8Pu3bv54XEigwcPZvPmzURHR7NgwQJatWoFwJVXXkmnTp1o06YNo0ePplu3bgwaNIgOHToQHR1Nv379Ct3Zz549m+zsbE8RTUlJYdq0aYXG3KNHD86ePUvr1q0ZM2YMiYmJADRp0oSJEyfSoUMHOnXqROvWrQudfuHChcyZM4fY2FiioqI8N9IoSkxMDKGhocTGxl50M4patWoRERHhKYrDhg3j6quv9qyLN954g+rVqzNv3jz69+9PdHQ0ISEh591QpDATJkwgLS2NqKgoli1bVujlIMXFNXHiRNLT04mJiWHMmDHMnz+/2PlFRkYyadIkunXrRkxMDCkpKRw6dKjEZfrxxx+TkpJSbNsiIsFMtTmfanPJcQVibXbe/DR+yTNx7k0gGagPfANMMLM5zrlbgRlAKDDXzMr04Ln4+Hj74e5rIgI7duwocmcciKKjo1m+fPl518Nc6Pjx49SoUQPnHIsWLeLNN98ssdBUNampqaSnpzNp0iR/h1KhPvvsM6ZPn87rr79e6OeFbS/OuXQzi6+I+IKVarPI+VSbVZsLo9rsm9pcIdfUmtnAIob/FfhrRcQgIoEpJSWF6OjoYosmQHp6OiNHjsTMqFOnjk8f6B0s+vbtS3Z2tr/DqHBZWVk8+eST/g5DRCRoqDb7jmqzb1TIL7XlxTl3G3BbixYthmdkZPg7HJFKI9iOBouUJ/1S61uqzSKFU20W8V5pa7M/r6m9ZGb2jpnd98P5+SIiIuJfqs0iIlLRArpTKyJFC+SzMEQqirYTEalI2ueIlKws24k6tSJBqHr16mRnZ6t4ihTDzMjOzqZ69er+DkVEqgDVZpGSlbU2V8iNosrLOdft+DsUkUqladOmfP3112RmZvo7FJFKrXr16jRt2tTfYYhIFaDaLOKdstTmgO7Umtk7wDvx8fHD/R2LSGUSFhZW4h0JRUREpOKoNouUH51+LCIiIiIiIgFLnVoREREREREJWAHdqXXO3eacm5WTk+PvUERERERERMQPXDDcgc05lwnsK3h7BVDWXm59IKuQdrz5u6h5FxePN7FeSj5QeE6lyaEy51TW9VKanLyNtSJz8iafsgwv7TjF0fcu+HLyNn5/bku+3IdfY2YNyhiHoNpcAu0jVZvLEqu+dxdTbVZtzmdmQfUCZl3CtJsLa8ebv4uad3HxeBPrpeRTVE6lyaEy51TW9VKanLyNtSJz8iafS81J3zv/fu8qY07exu/Pbak89uF6+eZV2dZrZdn2zm2rKu1PLjWnitiflDYnb/cnl5KTvnf+/d5Vxpy8jd+f21J57MMLewX06cdFeKcc2vHm76LmXVw83sTqq3zObas0OVTmnMq6XkqTk7exVrbvXVmGl3Ycb+l7V/T8Ayknb+MPxm1JLl1lW6+VZds7t62qtD/x5n1Rw4obXtbxStOOanPJw7z5rDTjeEO1ufTDihte1vFK045Pa3NQnH7sK865zWYW7+84fEk5BYZgyynY8gHlFAiCLR/JF4zrVTkFhmDLKdjyAeUUCCoqn2D8pfZSzPJ3AOVAOQWGYMsp2PIB5RQIgi0fyReM61U5BYZgyynY8gHlFAgqJB/9UisiIiIiIiIBS7/UioiIiIiISMBSp1ZEREREREQCljq1IiIiIiIiErDUqfWSc+5a59wc59wSf8fiK865Ps6515xzbznnuvk7Hl9wzrV2zr3inFvinHvA3/H4gnOulnNus3Oul79j8QXnXLJzbk3Bekr2dzy+4JwLcc5Nds694Jwb4u94LpVzLqlg/cx2zq3zdzy+4Jy72jn3tnNurnNujL/jEd9QbQ4Mqs2Vn2pz5afaXLwq0aktWFBHnHPbLhjewzn3pXNud0kL0sy+MrN7yzdS7/kop7fNbDhwPzCgPOP1ho9y2mFm9wO/BDqVZ7wl8UU+BR4H/lQ+UZaOj3Iy4BhQHfi6vGL1lo9y+i+gKXAGP+fko+1oTcF29C4wvzzj9YaP1lE0sMTMhgLtyi1Y8Zpqc+FUm8uXanORVJvLkWpzkXxWm6vE3Y+dc53J31AXmFmbgmGhwC4ghfwv+iZgIBAKTLmgiaFmdqRguiVm1q+iYi+Kj3P6/8BCM/u0gsIvlK9ycs71Bh4AXjezNyoq/gv5Ih8gFriS/CKTZWbvVkz0hfNRTllmluecawRMN7PBFRV/YXyU01DgWzN71d/7CB/vG/4E3Gtm31VQ+IXy0TrKBZaQ/4/b62Y2r2Kil6KoNqs2+4Nqs2qzP6g2l39t/klZJwwkZpbmnGt2weAEYLeZfQXgnFsE/JeZTQEq/akkvsjJOeeAqcDf/F00wXfrycyWA8udc38B/FY4fbSOkoFaQCRwwjn3VzPLK8+4i+PjbelboFq5BFoKPlpPXwOnC97mlmO4JfLVOnLOXQ3k+Ltogs/W0ShgQkFbSwB1av1MtVm12R9Um1Wb/UG1ufxrc5Xo1BbhZ8D+c95/DdxQ1MjOuSuByUA759zYgpVT2ZQqJ+Bh4BbgCudcCzN7pTyDK6PSrqdk4Hbyd8h/Lc/AyqhU+ZjZeADn3D0UHEUt1+jKprTr6HagO1AH+H25RlZ2pd2WlgEvOOeSgLTyDKyMSpsPwL1U7o5faXN6D5jonBsE7C3HuOTSqDarNvuDarNqsz+oNvuwNlflTm2pmFk2+de3BA0zmwnM9HccvmRmq4HVfg7D58zsD/6OwVfMbBn5hSZomNlx8gtN0DCzCf6OwZfMbBvg99NTxbdUmwODanPlp9ocGFSbi1YlbhRVhAPAVee8b1owLJApp8ov2PIB5RQIgi0fCM6cJDjXq3Kq/IItH1BOgSDY8gE/5lSVO7WbgJ8755o75y4D7gSW+zmmS6WcKr9gyweUUyAItnwgOHOS4FyvyqnyC7Z8QDkFgmDLB/yZk5kF/Qt4EzjEj7f0vrdg+K3k36FrDzDe33Eqp+DKKdjyUU7+j7Uq5hOsOekVnOtVOVX+V7Dlo5z8H2tVzKcy5lQlHukjIiIiIiIiwakqn34sIiIiIiIiAU6dWhEREREREQlY6tSKiIiIiIhIwFKnVkRERERERAKWOrUiIiIiIiISsNSpFRERERERkYClTq2IiIiIiIgELHVqRUREREREJGCpUytSBTjnejjnPi94feKc07YvIiLiR6rNIr7jzMzfMYhIOXPOZQCdzeyQv2MRERER1WYRX9IRIZGq4a/AVufcDH8HIiIiIoBqs4jP/MTfAYhI+XLOdQQc0MTMzvo7HhERkapOtVnEt/RLrUjw6w/sMrOzLl+4vwMSERGp4lSbRXxI19SKBDnnXAIwBzDgBPCgmaX7NyoREZGqS7VZxLfUqRUREREREZGApdOPRUREREREJGCpUysiIiIiIiIBS51aERERERERCVjq1IqIiIiIiEjAUqdWREREREREApY6tSIiIiIiIhKw1KkVERERERGRgKVOrYiIiIiIiASs/wXkkerFzkYjzAAAAABJRU5ErkJggg==\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 (conditioned)\")]:\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[1].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
|
|
}
|