{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Systems of linear equations\n", "\n", "> **⚠️ Warning**: The following Python examples are intended for educating concepts rather than for practical use. For practical implementations, consider using standard libraries such as `numpy`, `scipy` or `linalg`.\n", "\n", "Systems of linear equations appear in multitude of applications and solving them efficiently is important.\n", "Linear systems are ubiquitous in computational physics, such as solving circuit equations, performing data fitting, applying numerical solvers, and modeling physical systems.\n", "\n", "Generally a system of $N$ linear equations corresponds to the following construction\n", "\\begin{align}\n", "\\sum_{j=1}^N a_{1j} x_j & = v_1,\\\\\n", "\\ldots \\\\\n", "\\sum_{j=1}^N a_{Nj} x_j & = v_N.\n", "\\end{align}\n", "which can be written in a matrix form\n", "$$\n", "\\mathbf{A} \\mathbf{x} = \\mathbf{v}.\n", "$$\n", "where\n", "$$\n", "\\mathbf{A} = \\left(\\begin{array}{ccc} \n", "a_{11} & \\ldots & a_{1N} \\\\\n", "\\ldots& \\ldots & \\ldots\\\\\n", "a_{N1} & \\ldots & a_{NN} \\\\\n", "\\end{array}\\right)\n", "$$ \n", "\n", "A unique solution to the system of equations exists if they are all linearly independent.\n", "Alternatively, this is the case if the determinant of matrix $\\mathbf{A}$ is non-zero, $\\det \\mathbf{A} \\neq 0$.\n", "\n", "There are many efficient libraries that exist to solve systems of linear equations efficiently, and these should be used. Nevertheless, going over the basic methods for solving linear equations is important to understand the concepts, limitations, and possible issues and solutions to them." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gaussian elimination and backsubstition\n", "\n", "Gaussian elimination is the most basic approach of solving systems of linear equations.\n", "\n", "It is based on the fact that the following two operations leave the system of equations equivalent.\n", "1. One can multiply any equation by a constant (non-zero) factor, and it will still be the same system of equations with the same solution.\n", "2. One can take any linear combination of two (or more) equations to get another correct equation. This new equation can replace any of the equations entering this linear combinations. In other words, we can subtract from any of the equations any other equation, and the resuling system of linear equation will stay equivalent to what we had before.\n", "\n", "### Gaussian elimination\n", "\n", "Gaussian elimination simplifies a system of linear equations by systematically eliminating variables from equations, reducing the system to an upper triangular form. Let us take a system of 4 equations as an example\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "a_{11} & a_{12} & a_{13} & a_{14} \\\\\n", "a_{21} & a_{22} & a_{23} & a_{24} \\\\\n", "a_{31} & a_{32} & a_{33} & a_{34} \\\\\n", "a_{41} & a_{42} & a_{43} & a_{44}\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "v_1 \\\\\n", "v_2 \\\\\n", "v_3 \\\\\n", "v_4\n", "\\end{array}\\right)\n", "$$\n", "\n", "\n", "1. We start from the first row. We divide the row by $a_{11}$ such that its diagonal element is equal to unity\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & a_{12}/a_{11} & a_{13}/a_{11} & a_{14}/a_{11} \\\\\n", "a_{21} & a_{22} & a_{23} & a_{24} \\\\\n", "a_{31} & a_{32} & a_{33} & a_{34} \\\\\n", "a_{41} & a_{42} & a_{43} & a_{44}\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "v_1/a_{11} \\\\\n", "v_2 \\\\\n", "v_3 \\\\\n", "v_4\n", "\\end{array}\\right)\n", "$$\n", "\n", "2. Then, we make all entries in the first column below the main diagonal to go to zero. To achieve that we subtract the first equation multiplied by $a_{j2}$ from the $j$th equation:\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & a_{12}/a_{11} & a_{13}/a_{11} & a_{14}/a_{11} \\\\\n", "0 & a_{22} - a_{21} a_{12}/a_{11} & a_{23} - a_{21} a_{13}/a_{11} & a_{24} - a_{21} a_{14}/a_{11} \\\\\n", "0 & a_{32} - a_{31} a_{12}/a_{11} & a_{33} - a_{31} a_{13}/a_{11} & a_{34} - a_{31} a_{14}/a_{11} \\\\\n", "0 & a_{42} - a_{41} a_{12}/a_{11} & a_{43} - a_{41} a_{13}/a_{11} & a_{44} - a_{41} a_{14}/a_{11}\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "v_1/a_{11} \\\\\n", "v_2 - a_{21} v_1/a_{11} \\\\\n", "v_3 - a_{31} v_1/a_{11}\\\\\n", "v_4 - a_{41} v_1/a_{11}\n", "\\end{array}\\right)\n", "$$\n", "which we can denote as\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & a_{12}^{'} & a_{13}^{'} & a_{14}^{'} \\\\\n", "0 & a_{22}^{'} & a_{23}^{'} & a_{24}^{'} \\\\\n", "0 & a_{32}^{'} & a_{33}^{'} & a_{34}^{'} \\\\\n", "0 & a_{42}^{'} & a_{43}^{'} & a_{44}^{'}\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "v_1^{'} \\\\\n", "v_2^{'} \\\\\n", "v_3^{'} \\\\\n", "v_4^{'}\n", "\\end{array}\\right)\n", "$$\n", "\n", "3. Repeat steps 1-2 to make all elements below the main diagonal in the 2nd column go to zero\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & a_{12}^{''} & a_{13}^{''} & a_{14}^{''} \\\\\n", "0 & 1 & a_{23}^{''} & a_{24}^{''} \\\\\n", "0 & 0 & a_{33}^{''} & a_{34}^{''} \\\\\n", "0 & 0 & a_{43}^{''} & a_{44}^{''}\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "v_1^{''} \\\\\n", "v_2^{''} \\\\\n", "v_3^{''} \\\\\n", "v_4^{''}\n", "\\end{array}\\right)\n", "$$\n", "\n", "4. Repeat until all elements below the main diagonal are zero and all diagonal elements are equal to unity\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & \\tilde a_{12} & \\tilde a_{13} & \\tilde a_{14} \\\\\n", "0 & 1 & \\tilde a_{23} & \\tilde a_{24} \\\\\n", "0 & 0 & 1 & \\tilde a_{34} \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "\\tilde v_1 \\\\\n", "\\tilde v_2 \\\\\n", "\\tilde v_3 \\\\\n", "\\tilde v_4\n", "\\end{array}\\right)\n", "$$\n", "\n", "*Example:* System of equations \n", "$$\n", "\\left(\\begin{array}{cccc} \n", "2 & 1 & 4 & 1 \\\\\n", "3 & 4 & -1 & -1 \\\\\n", "1 & -4 & 1 & 5 \\\\\n", "2 & -2 & 1 & 3\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "-4 \\\\\n", "3 \\\\\n", "9 \\\\\n", "7\n", "\\end{array}\\right)\n", "$$\n", "becomes\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "1 & 0.5 & 2 & 0.5 \\\\\n", "0 & 1 & -2.8 & -1 \\\\\n", "0 & 0 & 1 & 0 \\\\\n", "0 & 0 & 0 & 1\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "-2 \\\\\n", "3.6 \\\\\n", "-2 \\\\\n", "1\n", "\\end{array}\\right)\n", "$$\n", "\n", "## Backsubstitution\n", "\n", "In the end we have the following system of equations\n", "\\begin{align}\n", "x_1 + \\tilde a_{12} x_2 + \\tilde a_{13} x_3 + \\tilde a_{14} x_4 & = \\tilde v_1, \\\\\n", "x_2 + \\tilde a_{23} x_3 + \\tilde a_{24} x_4 & = \\tilde v_2, \\\\\n", "x_3 + \\tilde a_{34} x_4 & = \\tilde v_3, \\\\\n", "x_4 & = \\tilde v_4.\n", "\\end{align}\n", "\n", "The solution now is trivial and proceeds through backsubstition, starting from \n", "$$\n", "x_4 = \\tilde v_4\n", "$$\n", "to\n", "$$\n", "x_3 = \\tilde v_3 - \\tilde a_{34} x_4\n", "$$\n", "to\n", "$$\n", "x_2 = \\tilde v_4 - \\tilde a_{23} x_3 - \\tilde a_{24} x_4\n", "$$\n", "and finally to\n", "$$\n", "x_1 = \\tilde v_4 - \\tilde a_{12} x_2 - \\tilde a_{13} x_3 - \\tilde a_{14} x_4.\n", "$$\n", "\n", "The algorithm complexity is $O(n^3)$.\n", "One limitation of simple Gaussian elimination is its susceptibility to numerical instability, especially for ill-conditioned matrices." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def linsolve_gaussian(A0, v0):\n", " # Initialization\n", " A = A0.copy()\n", " v = v0.copy()\n", "# A = A0\n", "# v = v0\n", " N = len(v)\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " # Divide the current row by the pivot (diagonal element)\n", " # Ensure diagonal element is non-zero to proceed\n", " div = A[r,r]\n", " if (div == 0.):\n", " print(\"Diagonal element is zero! Cannot solve the system with simple Gaussian elimination\")\n", " return None\n", " A[r,:] /= div\n", " v[r] /= div\n", " \n", " # Now subtract this row from the lower rows\n", " for r2 in range(r+1,N):\n", " mult = A[r2,r]\n", " A[r2,:] -= mult * A[r,:]\n", " v[r2] -= mult * v[r]\n", " \n", " # Backsubstitution\n", " x = np.empty(N,float)\n", " for r in range(N-1,-1,-1):\n", " x[r] = v[r]\n", " for c in range(r+1,N):\n", " x[r] -= A[r][c] * x[c]\n", " \n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we solve a system of 4 linear equations using the Gaussian elimination algorithm. The solution is verified by checking if $\\mathbf{A} \\mathbf{x} = \\mathbf{v}$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2. 1. 4. 1.]\n", " [ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. -2. 1. 3.]]\n", "x = [ 2. -1. -2. 1.]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "\n", "x = linsolve_gaussian(A,v)\n", "print(A)\n", "print('x =',x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pivoting\n", "\n", "Simple Gaussian elimination relies on the diagonal element of the present row being non-zero.\n", "This is not always the case: the diagonal element could be zero even in non-singular systems that have a perfectly valid solution. Simple Gaussian elimination will fail.\n", "\n", "Consider our example where we set the very first element to zero\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "0 & 1 & 4 & 1 \\\\\n", "3 & 4 & -1 & -1 \\\\\n", "1 & -4 & 1 & 5 \\\\\n", "2 & -2 & 1 & 3\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "-4 \\\\\n", "3 \\\\\n", "9 \\\\\n", "7\n", "\\end{array}\\right)\n", "$$\n", "\n", "The system has a solution but the solver will fail" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Diagonal element is zero! Cannot solve the system with simple Gaussian elimination\n", "x = None\n" ] }, { "ename": "TypeError", "evalue": "unsupported operand type(s) for *: 'float' and 'NoneType'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[3], line 8\u001b[0m\n\u001b[1;32m 6\u001b[0m x \u001b[38;5;241m=\u001b[39m linsolve_gaussian(A,v)\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mx =\u001b[39m\u001b[38;5;124m'\u001b[39m,x)\n\u001b[0;32m----> 8\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mAx =\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[43mA\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mv =\u001b[39m\u001b[38;5;124m'\u001b[39m, v)\n", "\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for *: 'float' and 'NoneType'" ] } ], "source": [ "A = np.array([[ 0, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = linsolve_gaussian(A,v)\n", "print('x =',x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution here is to apply the so-called partial pivoting. Pivoting involves swapping rows in the matrix to ensure that the largest possible pivot element is used at each step of Gaussian elimination. This reduces round-off errors and improves stability.\n", "\n", "Recall that the system of equations does not change if one swaps any two equations (rows). The solution thus is to swap the current row, if its diagonal element is zero, with one of the lower rows where the corresponding element is non-zero. (If all elements in all rows below are zero, we are dealing with a singular matrix that has no solutions). \n", "\n", "In pratice, at each step one chooses the row where the pivot element is largest in magnitude, even the pivot of the present row is non-zero. In this way one minimizes round-off error associated with the subtraction of large numbers.\n", "In our example from above,\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "0 & 1 & 4 & 1 \\\\\n", "3 & 4 & -1 & -1 \\\\\n", "1 & -4 & 1 & 5 \\\\\n", "2 & -2 & 1 & 3\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "-4 \\\\\n", "3 \\\\\n", "9 \\\\\n", "7\n", "\\end{array}\\right)\n", "$$\n", "we will swap the first and second rows to obtain\n", "$$\n", "\\left(\\begin{array}{cccc} \n", "3 & 4 & -1 & -1 \\\\\n", "0 & 1 & 4 & 1 \\\\\n", "1 & -4 & 1 & 5 \\\\\n", "2 & -2 & 1 & 3\n", "\\end{array}\\right)\n", "\\left(\\begin{array}{c} \n", "x_1 \\\\\n", "x_2 \\\\\n", "x_3 \\\\\n", "x_4\n", "\\end{array}\\right)\n", "=\n", "\\left(\\begin{array}{c} \n", "3 \\\\\n", "-4 \\\\\n", "9 \\\\\n", "7\n", "\\end{array}\\right)\n", "$$\n", "\n", "This procedure is performed at each step.\n", "To implement partial pivoting we will keep track of all the row swaps by keeping a map between the original row numbers and the current ones." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def linsolve_gaussian_partialpivot(A0, v0):\n", " # Initialization\n", " A = A0.copy()\n", " v = v0.copy()\n", " N = len(v)\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " # Find the pivot element (largest in magnitude)\n", " r_pivot = r\n", " for i in range(r + 1, N):\n", " if (abs(A[i][r]) > abs(A[r_pivot][r])):\n", " r_pivot = i\n", " \n", " # Swap rows to move the largest pivot element to the diagonal\n", " A[[r,r_pivot]] = A[[r_pivot,r]]\n", " v[[r,r_pivot]] = v[[r_pivot,r]]\n", " \n", " # Divide row r by the pivot element\n", " div = A[r,r]\n", " if (div == 0.):\n", " print(\"Diagonal element is zero! The system appears to be singular\")\n", " return None\n", " A[r,:] /= div\n", " v[r] /= div\n", " \n", " \n", " # Now subtract this row from the lower rows\n", " for r2 in range(r+1,N):\n", " mult = A[r2,r]\n", " A[r2,:] -= mult * A[r,:]\n", " v[r2] -= mult * v[r]\n", " \n", " \n", " # Backsubstitution\n", " x = np.empty(N,float)\n", " for r in range(N-1,-1,-1):\n", " x[r] = v[r]\n", " for c in range(r+1,N):\n", " x[r] -= A[r][c] * x[c]\n", " \n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we solve a system of linear equations using Gaussian elimination with partial pivoting. The results show that the solution is accurate, even when a diagonal element starts at zero." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 2. -1. -2. 1.]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = linsolve_gaussian_partialpivot(A,v)\n", "print('x =',x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 1.61904762 -0.42857143 -1.23809524 1.38095238]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 0, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = linsolve_gaussian_partialpivot(A,v)\n", "print('x =',x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LU decomposition\n", "\n", "Now that we have explored Gaussian elimination and pivoting, we move to LU decomposition, which builds upon these concepts for efficient solutions.\n", "\n", "LU decomposition factors $\\mathbf{A}$ into a lower triangular matrix $\\mathbf{L}$ and an upper triangular matrix $\\mathbf{U}$, simplifying the solution of systems of equations. This is particularly useful when the coefficient matrix $\\mathbf{A}$ remains constant while the right-hand side $\\mathbf{v}$ changes.\n", "\n", "LU decomposition is the following represetation of matrix $\\mathbf{A}$:\n", "\\begin{align}\n", "\\mathbf{A} = \\mathbf{L} \\mathbf{U},\n", "\\end{align}\n", "where $\\mathbf{L}$ and $\\mathbf{U}$ are lower and upper triangular matrices, respectively.\n", "\n", "At the end of Gaussian elimination our matrix became upper triangular\n", "$$\n", "\\mathbf{A x} = \\mathbf{v}, \\qquad \\Rightarrow \\qquad \\mathbf{U x} = \\mathbf{\\tilde v},\n", "$$\n", "\n", "Discarding pivoting for a moment, all steps of the Gaussian elimination can be represented by matrix multiplication, i.e.\n", "$$\n", "\\mathbf{U} = \\mathbf{L_{N-1}} \\ldots \\mathbf{L_{0}} \\mathbf{A}\n", "$$\n", "where e.g.\n", "$$\n", "\\mathbf{L}_0 = \n", "\\frac{1}{a_{11}}\n", "\\left(\\begin{array}{cccc} \n", "1 & 0 & 0 & 0 \\\\\n", "-a_{21} & a_{11} & 0 & 0 \\\\\n", "-a_{31} & 0 & a_{11} & 0 \\\\\n", "-a_{41} & 0 & 0 & a_{11}\n", "\\end{array}\\right),\n", "$$\n", "$$\n", "\\mathbf{L}_1 = \n", "\\frac{1}{a_{22}'}\n", "\\left(\\begin{array}{cccc} \n", "a_{22}' & 0 & 0 & 0 \\\\\n", "0 & 1 & 0 & 0 \\\\\n", "0 & -a_{32}' & a_{22}' & 0 \\\\\n", "0 & -a_{42}' & 0 & a_{22}'\n", "\\end{array}\\right),\n", "$$\n", "and so on.\n", "\n", "These are lower triangular matrices. Their inverses are also lower triangular matrices\n", "$$\n", "\\mathbf{L}_0^{-1} = \n", "\\left(\\begin{array}{cccc} \n", "a_{11} & 0 & 0 & 0 \\\\\n", "a_{21} & 1 & 0 & 0 \\\\\n", "a_{31} & 0 & 1 & 0 \\\\\n", "a_{41} & 0 & 0 & 1\n", "\\end{array}\\right),\n", "$$\n", "$$\n", "\\mathbf{L}_1^{-1} = \n", "\\left(\\begin{array}{cccc} \n", "1 & 0 & 0 & 0 \\\\\n", "0 & a_{22}' & 0 & 0 \\\\\n", "0 & a_{32}' & 1 & 0 \\\\\n", "0 & a_{42}' & 0 & 1\n", "\\end{array}\\right),\n", "$$\n", "and so on.\n", "\n", "Matrix $\\mathbf{A}$ can therefore be represented as\n", "\\begin{align}\n", "\\mathbf{A} & = \\mathbf{L}_0^{-1} \\ldots \\mathbf{L}_{N-1}^{-1} \\mathbf{U}\\\\\n", "& = \\mathbf{L} \\mathbf{U},\n", "\\end{align}\n", "where\n", "$$\n", "\\mathbf{L} = \\mathbf{L}_0^{-1} \\ldots \\mathbf{L}_{N-1}^{-1} =\n", "\\left(\\begin{array}{cccc} \n", "a_{11} & 0 & 0 & 0 \\\\\n", "a_{21} & a_{22}' & 0 & 0 \\\\\n", "a_{31} & a_{32}' & a_{33}'' & 0 \\\\\n", "a_{41} & a_{42}' & a_{43}'' & a_{44}'''\n", "\\end{array}\\right).\n", "$$\n", "\n", "This is the $LU$-decomposition of our matrix into a product of lower and upper triangular matrices. It is trivial to modify the Gaussian elimination code to calculate $L$ and $U$ matrices." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def lu_decomp(A):\n", " # Initialization\n", " U = A.copy()\n", " N = len(A[0])\n", " L = np.zeros((N,N), float)\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " # Record elements of the lower triangular matrix L\n", " for r2 in range(r,N):\n", " L[r2][r] = U[r2][r]\n", " \n", " # Divide row r by diagonal element\n", " div = U[r,r]\n", " if (div == 0.):\n", " print(\"Diagonal element is zero! LU decomposition without pivoting is not possible!\")\n", " return None\n", " U[r,:] /= div\n", " \n", " # Now subtract this row from the lower rows\n", " for r2 in range(r+1,N):\n", " mult = U[r2,r]\n", " U[r2,:] -= mult * U[r,:]\n", " \n", " return L, U" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "L = [[ 2. 0. 0. 0. ]\n", " [ 3. 2.5 0. 0. ]\n", " [ 1. -4.5 -13.6 0. ]\n", " [ 2. -3. -11.4 -1. ]]\n", "U = [[ 1. 0.5 2. 0.5]\n", " [ 0. 1. -2.8 -1. ]\n", " [-0. -0. 1. -0. ]\n", " [-0. -0. -0. 1. ]]\n", "LU = [[ 2. 1. 4. 1.]\n", " [ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. -2. 1. 3.]]\n", "A = [[ 2. 1. 4. 1.]\n", " [ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. -2. 1. 3.]]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "L, U = lu_decomp(A)\n", "print('L =', L)\n", "print('U =', U)\n", "print('LU =', np.dot(L,U))\n", "print('A =', A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "LU decomposition is useful for repeated solution of systems of linear equations\n", "$$\n", "\\mathbf{A x} = \\mathbf{v},\n", "$$\n", "when the matrix $\\mathbf{A}$ stays the same but where the vector $\\mathbf{v}$ can change.\n", "\n", "Indeed, the system of equations becomes\n", "$$\n", "\\mathbf{LUx} = \\mathbf{v}.\n", "$$\n", "Let us define \n", "$$\n", "\\mathbf{Ux} = \\mathbf{y},\n", "$$\n", "then\n", "$$\n", "\\mathbf{L y} = \\mathbf{v}.\n", "$$\n", "\n", "We can solve the system for $\\mathbf{x}$ in two steps.\n", "1. First we solve the equation $\\mathbf{L y} = \\mathbf{v}$ using forward substitution, in analogy to backsubstitution we used before.\n", "2. Once we have $\\mathbf{y}$, we can solve $\\mathbf{Ux} = \\mathbf{y}$ for $\\mathbf{x}$ using backsubstitution." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def solve_using_lu(L,U,v):\n", " # L*U*x = v\n", " # First solve L*y = v with forward substitution\n", " # Then solve U*x = y with backsubstitution\n", " # Initialization\n", " \n", " N = len(v)\n", " # Forward substitution for L*y = v\n", " y = np.empty(N,float)\n", " for r in range(N):\n", " y[r] = v[r]\n", " for c in range(r):\n", " y[r] -= L[r][r - 1 - c] * y[r - 1 - c]\n", " y[r] /= L[r][r]\n", " \n", " # Backsubstitution for U*x = y\n", " x = np.empty(N,float)\n", " for r in range(N-1,-1,-1):\n", " x[r] = y[r]\n", " for c in range(r+1,N):\n", " x[r] -= U[r][c] * x[c]\n", " return x" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 2. -1. -2. 1.]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "L, U = lu_decomp(A)\n", "\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = solve_using_lu(L,U,v)\n", "print('x =', x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The computational complexity of LU decomposition is $O(n^3)$, the same as Gaussian elimination. However, for repeated solutions with different $\\mathbf{v}$, the decomposition step is done only once, making it more efficient overall." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $LU$ decomposition with pivoting\n", "\n", "Not every non-singular matrix allows for $LU$-decomposition because diagonal elements can be zero.\n", "In general case we need to allow the possibility to perform partial pivoting by exchanging the rows of our matrix.\n", "If we do that, what we get $LU$-decomposition with pivoting which can be written as\n", "$$\n", "\\mathbf{P A} = \\mathbf{L U}.\n", "$$\n", "The permutation matrix $\\mathbf{P}$ records the row swaps performed during partial pivoting, ensuring that the LU decomposition is valid for singular and near-singular matrices.\n", "\n", "Solving the system of equations \n", "$$\n", "\\mathbf{A x} = \\mathbf{v},\n", "$$\n", "is also straighforward using forward and backsubstituion passes, except that we have to exchange the rows in the vector $\\mathbf{v}$ to account for the row swaps that we did. \n", "\n", "In cases where diagonal elements are zero or very small, partial pivoting improves numerical stability by selecting the largest pivot element, reducing round-off errors." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def lu_decomp_partialpivot(A):\n", " # Initialization\n", " U = A.copy()\n", " N = len(A[0])\n", " L = np.zeros((N,N), float)\n", " \n", " # Keep track of all row swaps\n", " row_map = [i for i in range(N)]\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " # Find the pivot element (largest in magnitude)\n", " r_pivot = r\n", " for i in range(r + 1, N):\n", " if (abs(U[i][r]) > abs(U[r_pivot][r])):\n", " r_pivot = i\n", " \n", " row_map[r], row_map[r_pivot] = row_map[r_pivot], row_map[r]\n", " U[[r,r_pivot]] = U[[r_pivot,r]]\n", " L[[r,r_pivot]] = L[[r_pivot,r]]\n", " \n", " # Record the elements of L\n", " for r2 in range(r,N):\n", " L[r2][r] = U[r2][r]\n", " \n", " # Divide row r by the pivot element\n", " div = U[r,r]\n", " if (div == 0.):\n", " print(\"Diagonal element is zero! The system appears to be singular\")\n", " return None\n", " U[r,:] /= div\n", " \n", " \n", " # Now subtract this row from the lower rows\n", " for r2 in range(r+1,N):\n", " mult = U[r2,r]\n", " U[r2,:] -= mult * U[r,:]\n", " \n", " return L, U, row_map" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "L = [[ 3. 0. 0. 0. ]\n", " [ 1. -5.33333333 0. 0. ]\n", " [ 2. -1.66666667 4.25 0. ]\n", " [ 2. -4.66666667 0.5 -1. ]]\n", "U = [[ 1. 1.33333333 -0.33333333 -0.33333333]\n", " [-0. 1. -0.25 -1. ]\n", " [ 0. 0. 1. 0. ]\n", " [-0. -0. -0. 1. ]]\n", "LU = [[ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. 1. 4. 1.]\n", " [ 2. -2. 1. 3.]]\n", "A = [[ 2. 1. 4. 1.]\n", " [ 3. 4. -1. -1.]\n", " [ 1. -4. 1. 5.]\n", " [ 2. -2. 1. 3.]]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "L, U, row_map = lu_decomp_partialpivot(A)\n", "print('L =', L)\n", "print('U =', U)\n", "print('LU =', np.dot(L,U))\n", "print('A =', A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The matrices A and L*U coincide up to a permutation of rows, as they should." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def solve_using_lu_partialpivot(L,U,row_map,v):\n", " # L*U*x = v\n", " # First solve L*y = v with forward substitution\n", " # Then solve U*x = y with backsubstitution\n", " # Initialization\n", " \n", " N = len(v)\n", " # Backsubstitution for L*y = v\n", " y = np.empty(N,float)\n", " for rr in range(N):\n", " r = row_map[rr]\n", " y[rr] = v[r]\n", " for c in range(rr):\n", " y[rr] -= L[rr][rr - 1 - c] * y[rr - 1 - c]\n", " y[rr] /= L[rr][rr]\n", " \n", " # Backsubstitution for U*x = y\n", " x = np.empty(N,float)\n", " for rr in range(N-1,-1,-1):\n", " x[rr] = y[rr]\n", " for c in range(rr+1,N):\n", " x[rr] -= U[rr][c] * x[c]\n", " return x" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 2. -1. -2. 1.]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "L, U, row_map = lu_decomp_partialpivot(A)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = solve_using_lu_partialpivot(L,U,row_map,v)\n", "print('x =', x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 1.61904762 -0.42857143 -1.23809524 1.38095238]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 0, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "L, U, row_map = lu_decomp_partialpivot(A)\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = solve_using_lu_partialpivot(L,U,row_map,v)\n", "print('x =', x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculating the matrix inverse\n", "\n", "The inverse $\\mathbf{A}^{-1}$ of matrix $\\mathbf{A}$ satisfies\n", "$$\n", "\\mathbf{A} \\mathbf{A}^{-1} = \\mathbf{I}.\n", "$$\n", "\n", "We can therefore evaluate $\\mathbf{A}^{-1}$ by solving $N$ systems of linear equations\n", "$$\n", "\\mathbf{A} \\mathbf{x}_k = \\mathbf{v}_k, \\qquad k = 1 \\ldots N\n", "$$\n", "where \n", "$$\n", "\\mathbf{v}_{k,j} = \\delta_{kj}\n", "$$\n", "and $\\mathbf{x}_k$ is the $k$th column of the inverse matrix $\\mathbf{A}^{-1}$.\n", "Note that each equation has the same matrix $\\mathbf{A}$, therefore, we can reuse the $LU$-decomposition $\\mathbf{A} = \\mathbf{L U}$.\n", "\n", "The computational complexity of inverting a matrix using LU decomposition is $O(n^3)$, the same as direct methods. However, LU decomposition is more efficient for reuse when solving multiple systems of equations." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def matrix_inverse_with_ludecomp(A):\n", " # First step: LU decomposition of matrix A\n", " L, U, row_map = lu_decomp_partialpivot(A)\n", " N = len(row_map)\n", " \n", " Ainv = A.copy()\n", " for c in range(N):\n", " v = np.zeros(N, float)\n", " v[c] = 1.\n", " x = solve_using_lu_partialpivot(L,U,row_map,v)\n", " Ainv[:,c] = x\n", " \n", " return Ainv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this section, we compute the matrix inverse using LU decomposition and verify the result by multiplying the matrix with its inverse. The result should be approximately the identity matrix.\n", "\n", "Test it with our matrix" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A*A^{-1} = \n", " ------------ ------------ ------------ ------------\n", " 1 -5.55112e-17 1.11022e-16 -1.11022e-16\n", "-2.77556e-17 1 -1.11022e-16 3.33067e-16\n", " 2.77556e-17 5.55112e-17 1 4.44089e-16\n", "-2.77556e-17 -5.55112e-17 -5.55112e-16 1\n", "------------ ------------ ------------ ------------\n" ] } ], "source": [ "# if unavailable run\n", "# pip3 install tabulate \n", "from tabulate import tabulate \n", "\n", "A = np.array([[ 0, 1, 4, 1 ],\n", " [ 3, 4, -1, -1 ],\n", " [ 1, -4, 1, 5 ],\n", " [ 2, -2, 1, 3 ]],float)\n", "\n", "Ainv = matrix_inverse_with_ludecomp(A)\n", "print(\"A*A^{-1} = \\n\", tabulate(A.dot(Ainv)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we test the matrix inversion with a random matrix of a larger size to evaluate the generality and scalability of the method." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A*A^{-1} =\n", " ------------ ------------ ------------ ------------ ------------ ------------ ------------ ------------\n", " 1 -5.99656e-17 -1.82773e-16 -1.52031e-16 -2.43055e-16 -2.73188e-16 -3.02291e-17 2.31419e-16\n", "-5.23166e-17 1 1.02429e-16 -6.06064e-17 4.62205e-16 3.66022e-16 -3.18536e-16 1.90081e-16\n", " 1.69258e-16 -6.81865e-16 1 -2.46065e-16 -2.12089e-16 -2.22237e-16 1.43684e-17 1.87233e-16\n", "-2.26756e-17 -3.71079e-16 3.77089e-17 1 -6.92619e-16 1.71737e-16 -9.67012e-17 -4.82271e-16\n", "-1.903e-17 -1.97257e-16 8.19161e-17 6.49946e-17 1 2.75247e-16 -2.0933e-16 -1.87141e-16\n", " 1.06784e-16 2.22085e-16 -4.1101e-16 -2.51684e-16 -5.81629e-16 1 8.40635e-17 8.99321e-16\n", " 8.02874e-17 5.08112e-16 2.91742e-18 -1.04034e-16 7.83533e-17 -4.42833e-18 1 3.41899e-16\n", " 2.53577e-17 -8.82038e-16 3.72835e-17 -7.87716e-17 -7.67222e-16 2.734e-16 -4.19679e-16 1\n", "------------ ------------ ------------ ------------ ------------ ------------ ------------ ------------\n" ] } ], "source": [ "# A is a square random matrix of size n\n", "n = 8\n", "A = np.random.rand(n, n)\n", "Ainv = matrix_inverse_with_ludecomp(A)\n", "print(\"A*A^{-1} =\\n\", tabulate(A.dot(Ainv)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tridiagonal systems\n", "\n", "A particular case of linear systems of equations is tridiagonal systems.\n", "In this case the matrix $\\mathbf{A}$ has non-zero elements only at the main diagonal, as well lower and upper subdiagonal:\n", "$$\n", "\\mathbf{A} = \n", "\\left(\\begin{array}{ccccc} \n", "d_1 & u_1 & & & \\\\\n", "l_2 & d_2 & u_2 & & \\\\\n", " & l_3 & \\ddots & \\ddots & \\\\\n", " & & \\ddots & \\ddots & u_{n-1} \\\\\n", " & & & l_n & d_n\n", "\\end{array}\\right)\n", "~.\n", "$$\n", "\n", "Tridiagonal systems of linear equations often appear in physics, e.g.\n", "- Nearest-neighbor interaction (linear chain of springs)\n", "- Finite differences applied to partial differential equations (heat equation)\n", "\n", "A tridiagonal system of linear equations can efficiently be solved using Gaussian elimination in linear time.\n", "This is because one only needs to process one row at a time (the rest already have zeros in the column), and at most two elements during the elimination process.\n", "Similarly, during the backsubstitution step, one only has to subtract a single element from the upper superdiagonal, i.e.\n", "\\begin{align}\n", "x_n & = \\tilde{v}_n, \\\\\n", "x_k & = \\tilde{v}_k - \\tilde u_{k} x_{k+1}, \\qquad k = 1\\ldots N-1.\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Solve tridiagonal system of linear equations\n", "# d: vector of diagonal elements\n", "# l: vector of elements on the lower subdiagonal\n", "# u: vector of elements on the upper superdiagonal\n", "# v0: right-hand-side vector\n", "def linsolve_tridiagonal(d, l, u, v0):\n", " # Initialization\n", " N = len(v0)\n", " a = d.copy() # Current diagonal elements\n", " b = u.copy() # Current upper diagonal elements\n", " v = v0.copy()\n", " \n", " # Gaussian elimination\n", " for r in range(N):\n", " if (a[r] == 0.):\n", " print(\"Diagonal element is zero! Cannot solve the tridiagonal system with simple Gaussian elimination\")\n", " return None\n", " b[r] /= a[r]\n", " v[r] /= a[r]\n", " a[r] = 1.\n", " if (r < N - 1):\n", " a[r + 1] -= l[r+1] * b[r]\n", " v[r + 1] -= l[r+1] * v[r]\n", " \n", " # Backsubstitution\n", " x = np.empty(N,float)\n", " \n", " x[N - 1] = v[N - 1]\n", " for r in range(N-2,-1,-1):\n", " x[r] = v[r] - b[r] * x[r + 1]\n", " \n", " return x\n", "\n", "# Returns tridiagonal vectors d,l,u of matrix A\n", "# See also scipy.sparse.diags at https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.diags.html\n", "def get_tridiagonal(A):\n", " n = len(A[0])\n", " d = np.zeros(n, float)\n", " l = np.zeros(n, float)\n", " u = np.zeros(n, float)\n", " for r in range(n):\n", " d[r] = A[r][r]\n", " if (r > 0):\n", " l[r] = A[r][r-1]\n", " if (r < n - 1):\n", " u[r] = A[r][r+1]\n", " return d, l, u" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [-2.04 0.08 -1.76 2.92]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 0, 0 ],\n", " [ 3, 4, -5, 0 ],\n", " [ 0, -4, 3, 5 ],\n", " [ 0, 0, 1, 3 ]],float)\n", "\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "\n", "x = linsolve_tridiagonal(*(get_tridiagonal(A)), v)\n", "print('x =', x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test for a random tridiagonal matrix" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def random_tridiagonal(n):\n", " A = np.random.rand(n, n)\n", " for r in range(n):\n", " for c in range(0,r-1):\n", " A[r][c] = 0.\n", " for c in range(r + 2, n):\n", " A[r][c] = 0.\n", " return A" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " x = [ 0.50650253 0.72652636 -3.37949485 ... 9.6947616 4.14125389\n", " -8.35797187]\n", "Ax = [0.53846347 0.584372 0.49716602 ... 0.76883773 0.4259966 0.02554295]\n", " v = [0.53846347 0.584372 0.49716602 ... 0.76883773 0.4259966 0.02554295]\n", "CPU times: user 2.57 s, sys: 37.2 ms, total: 2.6 s\n", "Wall time: 2.59 s\n" ] } ], "source": [ "%%time\n", "\n", "n = 5000\n", "A = random_tridiagonal(n)\n", "# print(\"A =\\n\", tabulate(A))\n", "v = np.random.rand(n)\n", "x = linsolve_tridiagonal(*(get_tridiagonal(A)),v)\n", "print(\" x = \", x)\n", "print(\"Ax = \", A.dot(x))\n", "print(\" v = \", v)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " x = [ 0.50650253 0.72652636 -3.37949485 ... 9.6947616 4.14125389\n", " -8.35797187]\n", "Ax = [0.53846347 0.584372 0.49716602 ... 0.76883773 0.4259966 0.02554295]\n", " v = [0.53846347 0.584372 0.49716602 ... 0.76883773 0.4259966 0.02554295]\n", "CPU times: user 57.4 s, sys: 847 ms, total: 58.3 s\n", "Wall time: 58.6 s\n" ] } ], "source": [ "%%time\n", "\n", "L, U, row_map = lu_decomp_partialpivot(A)\n", "x = solve_using_lu_partialpivot(L,U,row_map,v)\n", "print(\" x = \", x)\n", "print(\"Ax = \", A.dot(x))\n", "print(\" v = \", v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us apply it to the springs example from Section 6.1 of M. Newman *Computational Physics*\n", "\n", "We have linear chain of strings. The amplitudes $x_i$ of the vibration obey\n", "\n", "\\begin{align}\n", "(\\alpha - k) x_1 - k x_2 & = 0, \\\\\n", "\\alpha x_i - k x_{i-1} - k x_{i+1} & = 0, \\quad i = 2,\\ldots, N-1, \\\\\n", "(\\alpha - k) x_N - k x_{N-1} & = 0. \n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbqdJREFUeJzt3Xd4XOWZN/7vmTNNddTrqLkbjI0xAcxGYENwMIljUEwoWQK7gXcJIbHjZXkD+b3BySawmxBiZwlhScg6hLLZGJFGdcCyxWICBgtTbGNjyRrJ6m3Upp05vz9mnjOjPpLm9PtzXbqueDSSHpTRmfs8z104URRFEEIIIYTohEXtBRBCCCGEzAYFL4QQQgjRFQpeCCGEEKIrFLwQQgghRFcoeCGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdMWq9gKSLRwO48yZM8jIyADHcWovhxBCCCEJEEURg4ODKCkpgcUy/d6K4YKXM2fOoKysTO1lEEIIIWQOPB4P3G73tM8xXPCSkZEBIPIfn5mZqfJqCCGEEJIIr9eLsrIy6X18OoYLXthRUWZmJgUvhBBCiM4kkvJBCbuEEEII0RUKXgghhBCiKxS8EEIIIURXKHghhBBCiK5Q8EIIIYQQXaHghRBCCCG6QsELIYQQQnSFghdCCCGE6IrhmtQRkkyCIKC+vh5tbW0oLi5GdXU1eJ5Xe1mEEGJqFLwQMoXa2lps3boVLS0t0mNutxu7du1CTU2NiisjhBBzo2MjQiZRW1uLLVu2jAlcAKC1tRVbtmxBbW2tSisjhOiBIAioq6vDM888g7q6OgiCoPaSDIWCF0LGEQQBW7duhSiKEz7HHtu2bRtdjAghk6qtrUVlZSXWr1+PG2+8EevXr0dlZSXd9CQRBS+EjFNfXz9hxyWeKIrweDyor69XcFVEaXTnTOaCdm2VQcELIeO0tbUl9XlEf+jOmcwF7doqh4IXQsYpLi5O6vOIvtCdM5kr2rVVDgUvhIxTXV0Nt9sNjuMm/TzHcSgrK0N1dbXCKyNyoztnMh+0a6scCl4IGYfneezatQsT374gBTQ7d+6kfi8GRHfOZD5o11Y5FLwQMomrr74G53zle+Az8sY8nldYjD179lCfF4OiO2cyH9XV1SgtdU/5edq1TR4KXgiZxEsftsNbdB6Wb/sNnn95L67e/m8ovOF+fO0XL1LgYmB050zmg+d5XP4Pd0/6Odq1TS4KXggZRxRF/MdrJwEA/1i9CFdt+Az+4Ss3wVm+EoeaB1ReHZET5TuR+WjqHsbB8CLkX30v8grHBrhut5t2bZOIghdCxnn1aCeOtnmRZufxD39XCQC4sCoHAPBRmxdeX1DF1RE5TZfvBACiCPz0pz+lO2cyqX/9y0cICGF89vNfQFtLM17+619RdPXdKLzhfux7+wMKXJKIghdC4oiiiP/YF9l1+crFlchKtQMACjOdqMxNhSgC7zT1qblEIrOamhrceM/OCflOfEYe8q++B0PFa1RaGdGy14514NVjnbBaONy36WxYrVZsuPxyXLxhM5zlK9HQ4lV7iYZCgxkJiVN/ohvvefrhtFnw1U9XjfncBVU5aOoZwZuNPVi/rEClFRIlhMo/hdLbH8ctC31YlimguLgYJzg3fvjicfzg+Y9wblkWVpVlqb1MohH+kIDv//kjAMBXP12FRQXp0ufWVGTjndN9eKe5D19cM3UyL5kdCl4IifNwNNflxgsqkJfuGPO5C6py8T+HWvBWY68aSyMKEcIi3m8dAGfh8eWrN2JJYQYA4FJRxDvNA3jpw3bc8dS7eP6bn5Z25oi5/aq+EU09IyjIcOAbly8e87nzyrMBAO+eph3bZJL12OjAgQPYtGkTSkpKwHEc/vCHP0z7/Lq6OnAcN+Hj2LFjci6TEADAm6d68FZTL+y8Bf906YIJn2d5L++3DGAkEFJ6eUQhJzuHMBIQkGbnsTA/dgfNcRx+dO1KlOekorV/FP/8P+8hHJ4qO4aYxZn+Uemm596rliPdMXZP4LyKLADA8Y5BDFK+XNLIGrwMDw9j1apVePjhh2f1dcePH0dbW5v0sXjx4pm/iJB5YhegL33KjcJM54TPu7NTUOxyIhQWcbi5X+HVEaW85+kHAJzjdoG3jK06ynTa8MiXz4PdasGrxzrxWP0pFVZItOT+F45iNCjgU5XZ2HxuyYTPF2Q4UZaTAlEEGqKvLTJ/sgYvGzduxA9+8INZZ1gXFBSgqKhI+qDMfiK3d5v78PrJblgtHG6/dOGkz+E4Ttp9+RsdHRlWQ0s/AEyZ07Ki1IUdm84GAPz45eN0jGhib3zSjb8caYOFA3Z84ewpS+zXRI+O3qGjo6TRZLXR6tWrUVxcjMsvvxz79u2b9rl+vx9er3fMByGzxXZdas4rhTs7dcrnXVCVCwB4q7FHkXUR5bGdl3PdWVM+54YLynD1uSUQwiK+8cy76B7yK7M4ohlBIYzv/SmSpPvlCytwdolryueuqYjmvdCObdJoKngpLi7GY489hmeffRa1tbVYunQpLr/8chw4cGDKr3nggQfgcrmkj7KyMgVXTIzgg9YBvHasExYOuGPdommfe0F05+Vwcz/8IRrOZzS+oIBj7YMApt55ASK7cD+85hwsKkhHh9ePbf/dAIHyX0zltwdP43jHILJTbfjnDUumfe7q6M7L4dN9lCeVJJoKXpYuXYrbbrsN5513HtauXYtHHnkEn/vc5/Dggw9O+TX33HMPBgYGpA+Px6PgiokRsF2XL6wqQWVe2rTPXZifhtw0O/yhMI60ULddo/nwzACEsIj8DAeKXRPznuKlOaz4xZfPQ4qNx+snu/GzV08otEqitq5BP36692MAwN1XLpux6mxZUQZS7TwG/SGc6BxSYomGp6ngZTIXXXQRTpyY+qLgcDiQmZk55oOQRH3cMYiXPmwHxwFfXz/9rgsQueNmuy+U62A8DZ5IQLrKnTVl/kK8xYUZuL9mBQDgZ6+dQP2JLlnXR7ThRy8dw6A/hHNKXfjS+TPv9lt5C86N7uS920x5L8mg+eDl8OHDmhiCJggC6urq8Mwzz6Curg6CQEcGRsB2XTauKMLiaD+PmVxASbuGJeW7lE2dvzDeNavduOGCMogisO2/G9A+4JNpdUQL3m3uw+/faQEAfG/z2RMq0qbC8l4oaTc5ZG1SNzQ0hJMnT0r/bmxsRENDA3JyclBeXo577rkHra2teOKJJwBEpm1WVlbi7LPPRiAQwJNPPolnn30Wzz77rJzLnFFtbS22bt2KlpYW6TG3241du3bRrAodO9U1hL8cOQMgsV0X5sJo0u47Tb0ICWFYec3fA5AEvTdDpdFU7tt0Nt7zDOCjNi++8cy7ePq2i2Cj14XhCGER9/3xQwDAljVuqQFdIqhZXXLJ+td16NAhrF69GqtXrwYAbN++HatXr8Z3v/tdAEBbWxuam5ul5wcCAdx1111YuXIlqqur8frrr+P5559XNUCora3Fli1bxgQuANDa2ootW7agtrZWpZWR+Xqk7hOEReAzywumrRQYb2lRBjKdVgwHBHzURtVtRtE3HMDpnhEAwMrSrFl9rdPG45Evn4d0hxVvN/XhwVeOy7BCorb/OeTB+60DyHBY8X+vXDarr11dngUAONU9jN7hgAyrMxdZg5d169ZBFMUJH7t37wYA7N69G3V1ddLz7777bpw8eRKjo6Po7e1FfX09rrrqKjmXOC1BELB161aI4sTscPbYtm3b6AhJhzy9I3jucCsA4M7LZtcEkbdw+FQl5b0YDdt1WZCXBleqbdZfX5mXhh9vWQkA+M/9p7D3o45kLo+orH8kgB+9FOn2vu2KJcjPcMzwFWNlpdqlmUeHKe9l3mhfcxr19fUTdlziiaIIj8eD+vp6BVdFkuHR/Z9ACIuoXpwnJdLNBuW9GM97LFl3HgMXN55TjH/4u0oAwD//TwM8vSNJWBnRgof2foy+kSCWFKbjK2sr5vQ9qFld8lDwMo22trakPo9oQ/uAD78/FAlKvzHLXReGBS9vN/VS3waDkPJd3IkfIU7mno3LsaosC15fCF9/+l2M+AOU7K9zH53x4sk3TwOIdNKdaz4Tm3NEwcv8UfAyjUSrnLRQDUUS958HPkFACOPCqhwpCJmtFaUupNp59I8E8XHnYJJXSJQmiqJUaTSfnRcAsFst+PmNq+FKseHNV19EcWk51q9fjxtvvBHr169HZWUl5crpiCiKuO9PHyAsAp9bWYyLF+bN+XuxiqP3WvoRFMLJWqIpUfAyjerqarjd7in7PXAch7KyMlRXVyu8MjJXXYN+PP23SJL4XHddAMDGW6QLEeW96F9L3yh6hgOw8RyWF8+/V5Q7OxU12S3o+sP98PaMzX2hZH99+WPDGbzd1IcUG4/vXLV8Xt9rQV46XCk2+IJhHGujm575oOBlGjzPY9euXQAwIYBh/965cycNjtSRX71+Cv5QGOeWZeHvFuXO63tdUEl5L0bBjoyWF2fCaZv/37MgCPjVj++b9HOU7K8fQ/4Q7n/hKADgzssWoSQrZV7fz2LhcF606uid03TdmA8KXmZQU1ODPXv2oLS0dMzjbrcbe/bsoT4vOtI3HMBvD0bOrb95+aKEOqhOJ77T7mQVaUQ/pCOjaYYxzgYl++tXfEPSb/30KXQMjKAyNxW3Vlcl5fuzfi/v0JDGeZG1SZ1R1NTUYPPmzXj6jy9j++46pGfn4cTj22G10q9PT/7rfxsxEhBwdkkm1i8tmPf3W1WWBTtvQdegH009I6iaYS4S0a5kVBrFo2R/fZqsISmfkYd//MG/w2Fdn5SfIU2YpqTdeaGdlwTxPI/rvnAlMlesg1h8NrqHQ2ovicyC1xfEf73RBAD4xmXz33UBIo3JWJn1W4098/5+RB0hIYz3WyPBy2zGAkyHkv31Z6qGpMJgN7637dak5SitKsuChQNa+0dplMQ8UPAyC3arBRU5qQCAkzQZVFeeeKMJg74QlhSmY8NZRUn7vlK/l1N0fq1XJzqHMBoUkO6wYkFeelK+JyX768t0DUmZZOUopTmsUlI4DWmcOwpeZmlBfuTi9kkXBS96MewP4fHXGwFEZhhZEhyklogLF1DSrt6xfJeVblfSXhuU7K8vSuconUfN6uaNgpdZYu2daedFP57622n0jQRRlZeGz68sSer3Pq88G7yFQ2v/KFr6qJuqHjUkqb/LeJTsrx9K5yjRhOn5o4zTWaLgRR8EQUB9fT1Oe1rxk7p2iLmL8bV1CxMeX5+oNIcVK0pdeM/Tj7ebeuHOTk3q9yfya0hypVE8luz/jw88gRfeOoqrLz4bj/7L39OOi8YonaPEgpcPzwzAFxSSUp5vNhS8zJIUvNCxkWZNVjFgd+UDFzwCnF+W9J93YVUO3vP0463GXlyz2p3070/kMxII4eOOSLOwucy4SgTP87j88vXYP1QAsbiQAhcNYjlKra2tk+a9cBwHt9udtBwld3YK8tId6B7y44PWAZxfObdO32ZGx0aztCA/Ug7bNejHwGhQ5dWQ8aaqGAgMdOP6L31Jlq6m1KxOvz5o9SIsAoWZDhS5nLL9HJYIfKqbbnq0KD5HaTw5cpQ4jsMamnM0LxS8zFKm04bCzMgodEra1ZbpKwbk62r6qcoccBxwqmsYXYP+pH5vIq9kN6ebysKCyE1Pc88IzbTRKJajZM8cO7tIrhwlynuZHwpe5oDyXrRJra6mrlQblhVFSh9pzpG+NLBJ0jIdGTFFmU6k2nmEwiKaeymxW6s++7kvoOifHkfhDffjl//1BPbt24fGxkZZkqulZnXN/dShew4oeJmDRaxcmoIXTVGzq+mF0qgAalanJ2znRa58F4bjOKkDM103tKuxexichUfJ8vNx6y03Yd26dbLlKJ1d4oKN59A95Iend1SWn2FkFLzMwcIC6vWiRWp2NZWa1dHOi250D/nR0jcKjgPOcSens+50FuazvJdh2X8WmZvG6P83Soz6cNp4rCiNvO7eaabrxmxR8DIHbOeFjo20Rc2upp+KJu0e7xhE/0gg6d+fJN+R6JHRwvx0ZDptsv88lux/im56NEvJ4AUA1lCzujmj4GUOWM5Lc+8IfEEaaa8VYysGlO1qmp/hwIL8NIgicKiJLkR60MCGMcqcrMsslLpz086LVrHghQWacosNaexX5OcZCQUvc5Cf4UCGw4qwCDT10IVIS1jFQGbe2KnRSnQ1lfJemmgLWA9i+S7yHxkBtPOiB+z/mwUK7bycFw1ejrV7MeSnYb+zQU3q5oDjOCwsSEeDpx+fdA5LlSZEG2pqavDiYCn+/PJruHpJKq5ftwrV1dWyNwe7oCoHz7zlobwXHRBFEe8pVGnEsKOIvpEgeocDyEmzK/JzSWJEUZTykaqSNKBzJoWZTpRmpaC1fxTvefrxd4vyZv4iAoB2XuaMyqW1rbnXD2f5Stxw442yVgzEu7AqFwDwQesA3UVpXHPvCPpHgrDzFsVuPlLtVpRmpQCg3Rct6hkOYNAXAscBFbnKjfmgfi9zQ8HLHNGYAO0Kh0XpOK8qV5ntXwAoyUqBOzsFQljEu3Qh0jQ2z+iskkzYrcpdBtnREVUqas+paC5SaVaKorOGYv1e6JoxGxS8zNFCqjjSrDavD/5QGDaeQ0mWfC3fJ3OB1O+Fjo607L1osq7c/V3Gk8qlKWlXcxqjoxuUqjRizitnSbt9CIepWV2iKHiZI7bzcqpriF5wGtMUPbcuy0mFlVf2JX4hBS+6EMt3USZZl4ntvFDwojUs30WpZF1mWXEGUmw8vL4Q7cjNAgUvc1SWnQI7b4E/FEZrP3VH1BK1LkIAcEE076XB009l9BoVFML4oFXZMmkmtvNCb1Ja09jFyqSVSdZlbLxFCqIp7yVxFLzMkZW3SNuLdHSkLWznpVLBfBemMjcV+RkOBISwVIpLtOV4+yD8oTAynVbFXyNs56W5lwY0as0phRvUxZOOjijvJWEUvMwDmxRLW33aIgUvKlyEOI6joyONiy+Rtlgm78Ysl/gBjad7aECjVghhEad71AteqOJo9ih4mQcaE6BNjSoeGwGxvBfq96JNbEdM6SMjIBLcUrM67WntG0VQEGG3WlASLWdX0urozssnXcPoG6bxIomg4GUeFlKvF80JCWE090buaNXYeQFieS/vnO6jowENYpVGSjWnG29BHo0J0JpTrNIoNw28wrtxAJCTZpeC2sMe2n1JBAUv8xDf60UUqeJIC1r7RxEKi3BYLSjKVLZMmllckI6sVBtGg4KUGEq0YcgfwsedgwCAVQpMkp4MJe1qDytdV+PIiImVTPertgY9oeBlHhbkpYPjgP5ou2+ivvikO6XzGRiLhZOmTFPei7Z80DoAUQRKXE4UqBTcSsdG3bTzohXSNGmFBjJOhvJeZoeCl3lIsfNSu286OtIGNSuN4lHSrjZJ+S4qHRkB8dOl6ZqhFWrnyQGx4KXB048QHTfPiIKXeaIxAdqiZqVRPDbn6K2mXgjUxFAzlB7GOBl2NEE7ttohTZNWcedlUX46MpxWjAYFHGsfVG0dekHByzxRxZG2qNmgLt7y4gykO6wY9IVwnC5EmtHQ3A9AnUojJn7HlnZf1DcaEHBmwAdAuWnSk7FYOKnqiPq9zIyCl3liFUdUOaANbCCj2jsvVt4ibQP/rbFH1bWQiE6vD2cGfOA44ByVknUZKpfWDnbNcKXYkJ1qU3Uta8op7yVRFLzMEzs2+oR2XlQXCIXR2hcZ1VCZp9xI+6nQkEZtea8lUvm1uCAd6Q6rqmuhAY3aIeW75KeB49RJ8mcoaTdxFLzMEzs2au0fxbA/pPJqzK25dwRhEUiz88hPd6i9nDFJu1RKrz41m9ONtzCfunNrBdv9UrNMmllV5oKFA1r6RtHp9am9HE2j4GWestPsyEmzA6C7KLXFlzuqfQcFRI4mHFYLeoYDdKyoAVpI1mUW0M6LZmglTw4AMpw2LCnMAEB5LzORNXg5cOAANm3ahJKSEnAchz/84Q8zfs3+/fuxZs0aOJ1OLFiwAI8++qicS0yKRVT6qAlaKZNmHFZeajxFR0fqCodFaeflXA0EL+zY6HTvCAIhKotVU+zYSL1k3Xh0dJQYWYOX4eFhrFq1Cg8//HBCz29sbMRVV12F6upqHD58GPfeey+++c1v4tlnn5VzmfNGYwK0oVHFwWpTieW9UNKumpp6huH1hWC3WrC0KEPt5aAw04E0Ow8hLErjLIjyRFHURHfdeBS8JEbWrLWNGzdi48aNCT//0UcfRXl5OXbu3AkAWL58OQ4dOoQHH3wQX/ziF2Va5fwtouBFExo1dhECxg5pFEVRE8dZZsSOjFaUZMLGq39aznEcqvLT8EGrF590DUnXEKKsvpEgBkaDALSzY8t2az9o9cIfEuCw8iqvSJvU/yuOc/DgQWzYsGHMY5/97Gdx6NAhBIPBSb/G7/fD6/WO+VAaJd9pg1bKpOOtLs+G1cKhbcCHlmglFFGe2sMYJ0MVR+prjA5kLHE5kWLXRpBQkZuK3DQ7AkIYH7Qq/36mF5oKXtrb21FYWDjmscLCQoRCIXR3d0/6NQ888ABcLpf0UVZWpsRSx2B3TU09w9TWWSWjAQFtrNGURu6ggEhDspXRniJ/o7wX1TRoKN+FYdOlqdeLeljgqJV8FyCyK3deBRvSSEdHU9FU8AJgwrY6KzGdarv9nnvuwcDAgPTh8XhkX+N4Ja4UpNh4BAURp+n8WhVs1yUr1YbsaPWXVlzARgVQ3osqAqEwPjoTuYPVQpk0s7CAdmzVFj/IVUso72VmmgpeioqK0N7ePuaxzs5OWK1W5ObmTvo1DocDmZmZYz6UZrFw0oWI8l7UobVKo3gXLqBmdWo61u5FQAjDlWJDRa76zQsZtvPySdcw9QFSiRbz5IBY3ss7zX302piCpoKXtWvXYu/evWMee+WVV3D++efDZlO3bfNMaFKsurRYacSsqciGhQOaekbQQY2nFBc/SVpLCdNVeWngOGBglAY0qiW+u66WrHS7YLVw6Br0U67cFGQNXoaGhtDQ0ICGhgYAkVLohoYGNDc3A4gc+XzlK1+Rnn/77bfj9OnT2L59O44ePYpf//rXePzxx3HXXXfJucykoAGN6mJ3UFrcecl02rC8KA2+5iP48SO/Rl1dHQRBUHtZptEQTdY9V+V5RuOl2HmUuCIDGtnxBVGOEBalm54FKg5knIzTxuPs0sjrlZrVTU7W4OXQoUNYvXo1Vq9eDQDYvn07Vq9eje9+97sAgLa2NimQAYCqqiq88MILqKurw7nnnot//dd/xc9+9jNNl0kzNONIXSznpUpjd1AAUFtbi9f/9Xp0PHMvfnLP17F+/XpUVlaitrZW7aWZgpY6645HAxrVc6Z/FIFQGDaeQ2l2itrLmYCGNE5P1j4v69atm/a8bvfu3RMeu/TSS/Huu+/KuCp5LCoYe36tpe1pM2jsjiRKa6nSCIgELlu2bJnwd9Da2ootW7Zgz549qKmpUWl1xuf1BaWj3JUaStZlFuano/5EN42PUAE7MqrITQNv0d71+ryKLPz6f2nnZSqaynnRM/YHMOQPocPrV3s5pjLoC6J7KPI718I0aUYQBGzdunXSAJ49tm3bNjpCktEHLQMQRaA0KwX5GeoP6xxvIe28qKZRQzONJsMqjo62DdLQ30lQ8JIkdqsFFTmRN07Ke1FWU3TXJS/dgQyndhK76+vr0dLSMuXnRVGEx+NBfX29gqsyl4bokZGW+rvEiyX6086L0qRp0ho8agaAYlcKSlxOCGFROvokMRS8JFFsxtGgyisxl1ilkXZ2XYBITlcyn0dmL1ZppK1kXYY1R2umAY2K09I06alQs7qpUfCSRNKMI9oCVpRWe7wUFxcn9Xlk9qSxABrMdwHGD2ik3RclaW2a9GRYv5d3m/vVXYgGUfCSRNIWcCddhJTELkJa2/6trq6G2+2eMnmb4ziUlZWhurpa4ZWZQ/uAD+1eHywcsKJUmzsvHMdJb550dKQcX1BAa3+kf4oWe0MxayqyIYYFHNhfh6eeepraLMSh4CWJaOdFHVLworGdF57nsWvXLgATx1uwf+/cuRM8r42BcEbD8gSWFGYgzSFrYeW8xMqlKXhRyumeEYgikOG0Ildj40TiHf/bX3Hm0a/ik9134+///svUZiEOBS9JxCoHugb90ph1Ij8tTpNmampqsGfPHpSWlo55vNTtpjJpmb2nwWGMk6Hu3Mpj06QX5KVptq1FbW0trv/SlxAaHDuUmLVZMHsAQ8FLEmU4bSjMjJRj0oVIGX3DAfSPRAJFreW8MDU1NWhqasJfX30NhZv/BYU33I//PfwRBS4y03JzunjUqE55pzSe70JtFmZGwUuSSUdHVC6tCFZpVOxyIsWu3eMXnudx+WXrsfzTV8FZvhJnBqgXkJzCYRFHNJ6sy8SXS9MQPmWc0uhARobaLMyMgpckW5RPYwKUpNVKo6mURXsBNfeOqLwSYzvVPYRBfwhOmwVLCrV5d83QgEblSXlyGg1eqM3CzCh4STLaeVEWuwhpMd9lMuU5kRkqHgpeZMWGMZ5T6oKV1/ZlzmnjUZoVeV1QxZEytDpNmqE2CzPT9l+1DlHynbK03uJ7vHLaeZGVIAioq6vDU08/DV/zEZxTnKH2khLCci8o70V+/SMBaYdLqzu21GZhZhS8JBnbeWnuHYEvaN5kKqVoudJoMhS8yKe2thaVlZVYv349/vDQ/0XHM/fiZ//0WV1UZbDgmyWSEvmw33FRplOzJfTUZmFmFLwkWX6GAxlOK8Ji7I2VyEMURTR2aXM0wFRYzgsdGyUXm949Psmxt7NdF2WlbLQI5crJr1HjybrMVG0W3NRmAQAFL0nHcRzlvSika8iP4YAACxcLCrSOrbN7KECTYpPECGWlC2nnRTFaz3eJx9osfOunTyJv07/gi//vMTQ2Npo+cAEoeJEFjQlQBpsmXZqdAodVH9unmU4bslIjk689fbT7kgxGKCtdWEADGpVyKtqgTus7LwzP8/jM5Zch7axLES46y9RHRfEoeJEBjQlQht7KpBkp76WHgpdkMEJZaUEGDWhUCuvxooedF6YsO3rc3Deq8kq0g4IXGbBeL3RsJC89jLSfDPV6SS4jlJVyHBfLe6FyadmEw6KUi7ggT9v9f+KxG56uQT9GA9o9/lQSBS8yYDsvp7qGIISpY6ZcmnTW44Upp6TdpDJKWSkLwqnNgnzavD74gmFYLRzc2SlqLydhrlQbMp2Ryig6bo6g4EUG7uwU2HkL/KEwzvTTNp9c9FYmzUjBC20BJ4VRykpjvV5o50UurNKoPDdV880Lx6NKxbH09f+eTlh5i5QMRkdH8giHxViLb73mvNBFKGlYWWlRccmYx/VUVkoNLuUXP01ab+i6MZY2O/QYwKKCdBzvGMTJziGsX1ag9nIMp93rgz+kv+1fYOyxUTgswmKZ/LiDzE5NTQ3yV/wdrt3xa+RwI3jktstRXV2t+R0XJjZdOjKgcapjMDJ3Wp8mPR0KXsai4EUmC/Np50VOLN+lPEd/27/FLid4Cwd/KIyuIT8KM51qL8kwzgz44SxfifOX5GPdugvUXs6sxA9o7BkOIC/dofaSDEfr06Sn46ZjozH0ddXXkVjlAAUvcjil02RdIHKsyAbx0V1Ucnl6I3lEZTrbjQPGDmikvBd5aH2a9HRiO7aUKwdQ8CKb+F4vk3X+JPPTpOOLEEC9XuTCgsFynXRcHm8hDWiUjT8koCVaqaOnHi9M/LERvadQ8CKbBXnp4DigfySyBUySS6+VRgz1epEH+33qZVzEeOxNlXZsk6+5ZwRhEUh3WJGvwyO5kiwnOA4YDQroHqL3FApeZJJij20B07C15Dul00ojpiwn8tqg8+vkYnfWet15oXJp+ZyK263VYzK0w8qjOJofR71eKHiRFY0JkEdICEtv+lU63P4FqHJADsP+kHRHytqp681C2nmRjZ7zXRjq9RJDwYuMaEyAPM70+xAURDisFulORG8oeEm+lmjTv0ynFa7o8Eu9YTkvnr5RGtCYZI06nGk0XhnlykkoeJGRtPNCwUtSNUbzXSpyU3XbI4UFL52DfviCNKskGaRk3Vx97roAkQGN6Q4rDWiUgd6mSU8m1p2bghcKXmS0sIDOr+XQ2KX/i5ArxYaM6KySFroQJQXbStfrkREQGWewQOoRRdeNZGrs1t9AxvFoxzaGghcZsWOj1v5RDPtDKq/GOJqiW6Z6rTQCIm9SdCFKLr2XSTNSuXQ37dgmy8BoUMqH0mueHBCf6E+9Xih4kVF2mh25aXYAtPuSTHqdaTQe9XpJLraD5dZ58CJNl6adl6Rh1wx2LKdXLOelbYByoih4kdlCqeJoUOWVGEejjrvrxovtvNBdVDIYZedlAe28JF2jAfJdACA/3QGnzYKwCJzpN/d1g4IXmUmTYukuKikCoXCsS6bOL0TUqC55RFHU9WiAeAsL2M4LdedOFiNUGgGR42aW02X26wYFLzKjiqPk8vRFumSm2XnkZ+ivS2a8curZkDTdQwGMBgVwHFCq8+ClMjcyoNHrC1F37iQ5ZYBkXYYqjiIoeJEZNapLLnYHVZGrzy6Z8WhWSfKwC3lRphMOK6/yaubHaePhzqYBjcmk52nS49GObQQFLzJjwUtT9zCCgrkTrJKBzTTSc8UAU5KVQrNKksSj85lG47EdAuq0O3+iKMaS/A1w3aAd2wgKXmRWnOlEio1HKCyaPlJOBqNUGgGA3WpBiStyh02vjfkxQo+XeDRdOnk6vH6MBgXwFk73ydxA/IgAStiV3SOPPIKqqio4nU6sWbMG9fX1Uz63rq4OHMdN+Dh27JgSS006i4WTEvAo72X+jFJpxNCAxuQwSqURE5suTcdG88UCwPKcVNh4/d+vU3+oCNn/n/zd736Hbdu24Tvf+Q4OHz6M6upqbNy4Ec3NzdN+3fHjx9HW1iZ9LF68WO6lyoZmHCVPkwGGq8WjC1FySJVGOfpO1mVY8EI7L/N3ymDXDPYaHxgNYmA0qPJq1CN78PLQQw/hq1/9Km699VYsX74cO3fuRFlZGX7xi19M+3UFBQUoKiqSPnhev0l4Urk0XYjmxRcUcGbAB8A4FyIKXpLDaDsv7IanuXcE/hDNvpoPI0yTjpdqtyIvPdL81Mw7trIGL4FAAO+88w42bNgw5vENGzbgjTfemPZrV69ejeLiYlx++eXYt2+fnMuUHUva/YR2XuaFJetmOq3I1unU4PFoxP38BYUw2gbYzosxgpf8aCfYsEgdmOdLmmlkgGRdhq4bMgcv3d3dEAQBhYWFYx4vLCxEe3v7pF9TXFyMxx57DM8++yxqa2uxdOlSXH755Thw4MCkz/f7/fB6vWM+tEYKXrqGqSR2HuKPjPReJs1Q5cD8nekfRVgEHFYL8tP13fuH4TgOCynvJSlOGWCQ63i0YwsoMuRh/BuNKIpTvvksXboUS5culf69du1aeDwePPjgg7jkkksmPP+BBx7A9773veQuOMkqctPAWzgM+UNo9/pQ7DLGubzSGrsjf6hGvAi1eX3whwTd9yhRQ3NcmbTFYoygFoiMCXivZYCOm+chEArD0xfZlTNCgzqGuuzKvPOSl5cHnucn7LJ0dnZO2I2ZzkUXXYQTJ05M+rl77rkHAwMD0ofH45nXmuVgt1pQEX2TojEBc9dksEojAMhJsyPNzkMUgdY+c5c+zpVRxgKMt1BK2qVrxlx5+kYghEWk2nkUZhpjVw6I77Jr3muGrMGL3W7HmjVrsHfv3jGP7927FxdffHHC3+fw4cMoLi6e9HMOhwOZmZljPrRIGtDYSQMa58poiXdAdFYJbQHPi9GSdZkFlOg/b/GddY1y1AxQzgugwLHR9u3bcdNNN+H888/H2rVr8dhjj6G5uRm33347gMjOSWtrK5544gkAwM6dO1FZWYmzzz4bgUAATz75JJ599lk8++yzci9VVosK0rH3ow4aEzAPjT3GC16AyJvusfZBU1+I5oONBjBKsi4TXy493VE7mZpRpkmPx8qlW6I7S7yBjksTJXvwct1116Gnpwff//730dbWhhUrVuCFF15ARUUFAKCtrW1Mz5dAIIC77roLra2tSElJwdlnn43nn38eV111ldxLlRX1epmfIX8IXYN+AMY6NgJoVsl8GW00ADN+QGOeQZKRlSRVGhnsmlHsSoHVwiEoiOjw+lCSZawj00QokrB7xx134I477pj0c7t37x7z77vvvht33323AqtS1sK4iiMyeyzfJTfNjkynMcqkGaocmB+jjQZg2IBGT+8oPukcouBlDtixETuCMwrewsGdnYKmnhE0946YMnjRf69knWDJd12DflN3RZwrI+a7MLHgxbzJd3M16AuibyTy92SU7rrxpBlH3XTTMxdG664bz+w7thS8KCTDaUNRphMAHR3NhRErjZj45DvqAzQ7rNIoO9WGDIPtyAFx06XpmjFrg76gdNRshGnS47HrRgsFL0Ru1Gl37oy88+KOlvgO+UPSLgJJjFErjRg21JV2XmavKdoXKi/dYbijZoCOmyl4UVCsYyYFL7Nl1EojIJLbwHblqOJodlqilUZugwYv0s4LXTNm7VS00shoybqM2RvVUfCioEUFVHE0V9KxUa4xL0Rmv4uaK8PvvERveDw0oHHW4nu8GJHZG9VR8KIgqVEd3UXNSv9IQDpOqcwz5puU2ZPv5sqolUZMfoYDGTSgcU6ko2YD5rsAseCla9CP0YD5AlsKXhTEdl48vSPwBc33YpsrdhEqzHQg1a5Idb/iaEDj3Bh954XjOKlZHR0dzY5Re7wwrlQbMpyR6yFr1GgmFLwoKD/dgQxn5C6qqYcS8BLVZOB8F6Y8N5K0SzsviQuHRbREt8yNGrwAsXJp6hGVOFEUpWnSCwy68wKY+6aHghcFcRyHhXkp8DUfweO/eRJ1dXUQBNqBmYkRp0mPRzkvs9c15Ic/FIaFA4qznGovRza08zJ7XYN+DAcEWDjjdV6OZ+brhjH34DWqtrYWr+64A0M9HXjwGeBBAG63G7t27UJNTY3ay9OsRoMn6wKxnI0z/aMICmHYeLqvmAm7YJdkpRj69yU1qqOdl4Sx0vKynFQ4rLzKq5GPmXPljPsXrzG1tbXYsmULhno6xjze2tqKLVu2oLa2VqWVaV+TgXu8MPkZDjisFoTFSABDZmb0ZF1mgRS8DFETwwQZvdKIiTW4NN81g4IXBQiCgK1bt0564WGPbdu2jY6QJiGKoimCF47jTL0FPBdGT9ZlKnJTpQGN3UMBtZejC0adJj0e5bwQWdXX16OlpWXKz4uiCI/Hg/r6egVXpQ/dQwEM+kPgDH52DZj7/Hou2N2mEWcaxXPaeGl36RTlvSTE6JVGTFl2LNHfbLtyFLwooK2tLanPMxNWaVSalQKnzbhn14C5z6/nQjo2MnhQC8Qn7VLeSyJYzovRpkmPV5qdAo4DRoOC6XblKHhRQHFxcVKfZyZGnmk0Xrk0aM1859dzwXpbmCF4WRiX90KmFxTCUkM/o183HFYexWy0iMl6vVDwooDq6mq43W5wHDfp5zmOQ1lZGaqrqxVemfaZodKIoWOjxPlDAtq9PgDGz3kBqFx6Nlr6RhEKi3DaLNLMMCNzmzTvhYIXBfA8j127dgHAhACG/Xvnzp3geWMfi8yFNNPI4HdQAFCeS8FLolr7RiGKQIqNR26aXe3lyE7aeaHp0jOKJeumw2KZ/IbRSKSbHpONj6DgRSE1NTXYs2cPSktLxzzudruxZ88e6vMyBbMk3gGxkt+B0SAGorOcyOTiK42m2tE0kgX5aRDDAj4+/CaeePIpanA5DVYmbYZrBhA/oNFcwQs1qVNQTU0NNm/ejIee+CN+9NybWFRZhtd/egftuEwhHBalhF0z7Lyk2HnkZzjQNeiHp28ErlSX2kvSLDZJ1+iVRkz9K8/jzH/ejpC3Gzc/E3mMGlxO7pSJ8uSA2N+A2XZsaedFYTzPY9PGzyDtrEsxkruUApdpdAz64AuGwVs4uLPN8SZFeS+JMVOlUW1tLa699lqEvN1jHqcGl5NrNEmDOqbcpI3qKHhRAXux0fHA9NiRUXlOqqHbv8eL79tApuYxSYM6anA5e9JRs4EHMsZjAXzbwCgCobDKq1GOOd4RNCbVbkVeugMAvUlNJ1ZpZOw3qHi085KYZpOMBqAGl7Mz7A9JVWhm2XnJT3fAaTPfaBEKXlRSQZUlMzJTpRFTZtKyx9mSdl4MHthSg8vZYTc8OWl2ZKUavwoNiLbayDbf+wkFLyphd9ine6n0cSqN3ZE/RLNUDQC085KIgZEgvL4QABg+F4oaXM6OmaoT45mx4oiCF5XQHfbMzFRpxLCdhNa+UQhhc80qSRS7QOelO5BqN3bBJDW4nB2zTJMez4yjRSh4UUkF23kxWWOhRAlhUWq6ZIbuukxhhhN23oJQWETbgHnOr2dDyncxQZk0NbhMnCAIqD9Qh+GP9iPY+oGpkpjNeDNMwYtKqJvq9M70jyIghGG3WlCSZfw3KcZi4eA2ad+GRDWbpNKIoQaXM6utrUVlZSWe/df/g+4//xi7vvX3qKysNE0ZuRnLpSl4UQnbeTnTP4qgYJ7ytkSxs+uKnFTwJmjxHa/chHdRs+ExSaVRvJqaGjQ1NeHq7/wn8jb9C+7a9RQaGxspcEEkcNmyZcuEqiwz9cExY64cBS8qyc+Ilbe19pknWk5UowkrjRgzXohmw2w7LwzP87i4+hKknXUpnOUr6agI1AeHYYnrA6NBDIyao3cYBS8q4TiO3qSm0WiyFt/xYq8LCmon0xIN9t0myHkZryIn8vdwmq4ZAKgPDpPmsCIvPVIabpYdWwpeVBQrlzbHi202WKWRGYMXM1YOJEoIi2jpM+fOCxCXK9dDLRYA6oMTz2xJuxS8qMhsL7bZiHXXNWHwkk2vi6l0eH0ICiKsFg7FLhPuvESDl5a+UYQoV4764MQxW6M6Cl5UxJJ2m6lceoygEJaOBsy58xJ5U+4dDmDQZ47z60SxC3NpdorpErmBaCm9lZXS+9RejuqoD06M2RrVUfCiIrYFTMdGMYIg4Pd/fhneD+oQbv0QeWnGbkI2mQynDTlp7Pya8l7imbHSKJ7FwknDO6lH1Ng+OOOZrQ+O2XLlKHhRUXk0+c7TOzJptrzZsF4NX77mKnT/+cfwPPl/UVVVZYpSx/Eo72VyUvBiwnwXpiKXJe1S3gsQKSP//e9/D2tm3pjHzdYHhyWwm+W42Xy3tRrCytuG/CH0DgeQG500bUasV8P4II71ajDTRQiI3EW95+k3zYUoUZ7ocaIZk3UZqlKc6LIrN6Hkn2zwt3yIf99YjoqyUlRXV5tix4Vhr4uWvhEIYdHwx6q086Iip41HUaYTgLkvRNSrYaJy6rI7KTONBpiKNJGejo0kp3tHwFl4VKz4FG6+6ctYt26dqQIXACh2pcBq4RAURHR4jZ8PRcGLymhMAPVqmIzZku8S5TFpg7p4LHihnJeY0z2sI7f5EvwZ3sKhNNs8Nz0UvKisnCqOqFfDJCjnZSJfUEDnoB+AeRN2gViuXDPlyklYIMcCO7My05EiBS8qqzDRi20q1KthIun8uncU4TC9QQGQmtNlOKzISrWpvBr1uLNTwHGxXDlCwQtTJl03jP9+okjw8sgjj6CqqgpOpxNr1qyZcft///79WLNmDZxOJxYsWIBHH31UiWWqgsqlqVfDZNj5dUAIo2PQ+OfXiWABvjsndcrXihnE58qZ+boRTzo2MmFTy3hmalQne/Dyu9/9Dtu2bcN3vvMdHD58GNXV1di4cSOam5snfX5jYyOuuuoqVFdX4/Dhw7j33nvxzW9+E88++6zcS1UFTRAe16th3JuS2Xo1MLyFk6rRzHykGI/1vCk3cbIuQ8fNY7Egzuw7L3RslEQPPfQQvvrVr+LWW2/F8uXLsXPnTpSVleEXv/jFpM9/9NFHUV5ejp07d2L58uW49dZb8Y//+I948MEH5V6qKtiLrd3rgy9onmqa8WpqarBnzx5k5xWOedxsvRriUd7LWM0mb1AXj5J2Y4b9IXRFc6HMnLALxCf6G79Rnax9XgKBAN555x18+9vfHvP4hg0b8MYbb0z6NQcPHsSGDRvGPPbZz34Wjz/+OILBIGw2Y51156TZkWbnMRwQ0NI3ikUF6WovSTU1NTV4R1yIx599AZeV2/F/rlxjul4N8Wj21VgseCk3+d01EDseocA29jvISrXBZeJcKCDWQqBr0I/RgIAUu3GvnbLuvHR3d0MQBBQWjr2bLiwsRHt7+6Rf097ePunzQ6EQuru7Jzzf7/fD6/WO+dATjuNQnhvrtGt2zf0+OMtX4uovfsmUvRrimWkLOBFmHw0QL/baoC67UrKuicvnGVeKDRnOyJ6E0dssKJKwOz65ThTFaRPuJnv+ZI8DwAMPPACXyyV9lJWVJWHFymJn+KdpzD1VDcSh4CVGFEUaDRCHjo1iKFk3huM40+RRyhq85OXlgef5CbssnZ2dE3ZXmKKiokmfb7VakZubO+H599xzDwYGBqQPj8eTvP8AhcS2gI1/TjmdQCgslcNWmnCa9HhmG7Q2nb6RIIYDkZwwlshsZiy3ozN6PGBmlKw7llkqjmQNXux2O9asWYO9e/eOeXzv3r24+OKLJ/2atWvXTnj+K6+8gvPPP3/SfBeHw4HMzMwxH3pTRlvAAIDW/lGERcBps6Agw7xznhj2uuge8mMkEFJ5NepiF+LCTAecNvMeJTKuVBtcKZHrodHfpGZCOy9jmaVru+zHRtu3b8evfvUr/PrXv8bRo0fxrW99C83Nzbj99tsBRHZOvvKVr0jPv/3223H69Gls374dR48exa9//Ws8/vjjuOuuu+ReqmroeCCiKXoRqsxNM3UfD8aVEnuDajFB9cB0aCzAROx3YfbjZjpqHiuW6G/sa4bsU6Wvu+469PT04Pvf/z7a2tqwYsUKvPDCC6ioqAAQafke3/OlqqoKL7zwAr71rW/h5z//OUpKSvCzn/0MX/ziF+Veqmriu+zOlA9kZKe72R0UXYSY8pxUvN86gOaeESwpzFB7OaqhMumJynOjrw0T3/QEQmGc6Y+8SdN1I6Iseqxq9JwX2YMXALjjjjtwxx13TPq53bt3T3js0ksvxbvvvivzqrSjJCsFFg7wBcPoGvSjINo902yaeijfZTwpeDH4hWgmLBeKknVjKnIoabelbwRhEUi188hPp6NmYOxOvpFvhmm2kQbYrRaUZJlnGuhUTscdG5EIalQX0UyVRhNUmCS3YTqn444TjfomPVul0dlXo0EBPQaefUXBi0aU010UnV1PwixljzOJjQag1wYTP13arOioeSKHlUdxdPfeyK8NCl40wuxJuyEhLP23085LjNlfF0DktdEazWsoo7lGEvaG3dI3AsGkk8djZdJ0zYjnNsFNDwUvGmGW8rapnOn3IRQWYbdapIm5JPZmzc6vzahtwAchLMLOW1CYQa8NpijTCTtvQVAQpaRVs6Hd2smZYceWgheNMPsdNiuTrshJhcVCZ9cMS+b2h8LS8DmzYRdgd3YKvTbiWCwc3DnmzpWTeryYfCDjeGZ4P6HgRSMqTH5+TY2mJmfjKZnbQ5VGUzJzxZEQFqVcKNp5GavMBEEtBS8awSLlrkFzdlOVyqTpIjSBGe6iphOrNKJ8l/FYsH/ahN25270+BIQwbDwnBfgkotwEjeooeNGI+HbfRn7BTUXaeaEeLxOYPXihSqOpSa8NE+68sGuGOzsVPB0njsF2KdsGRhEIhVVejTwoeNEQM7f7pp2XqZml3fdUqLvu1Mzc64WSdaeWn+6A02ZBWIRhk7kpeNEQs95hC2FRunOkMumJzFA5MB3qrjs1KXjpMV81mhS80OtiAo7jDD9dmoIXDWHl0mZ7k6Kz6+mZNagFgGF/CN1DkS6hFLxM5M5OBccBg/4Q+kaCai9HUZTkPz3ppqfPmNcNCl40RDo2MtmbFOuSWZZDZ9eTYa+Ldq8PvqCg8mqUxS688RO2SYzTxkt9kcx23EzHRtMz+mgRCl40pMLgL7apNNJMo2llpdqQ4YjMUG3pM+b59VQoWXdmRn+TmowoirTzMoMygx83U/CiIezF1tI7aqp233QHNT2O4wx/IZoKlUnPzIy9XnqGAxgOCOA4em1Mxejl0hS8aEixywmrhUNACKPD61N7OYpp6qadl5mYoenUZDy9lKw7Exb0myl4Yf+txZlOOKy8yqvRJqNfMyh40RArb4E7O/KCM+OFiHZepmbWpF0PlUnPqDyXdec2T84LHRnNjP3NDIwGMTBqvGRuCl40xmzHA+GwKHUHpZ2XqZk2eIkm7FLOy9TMmCtHNzwzS3NYkZduB2DM9xMKXjTGbE2nOgf98AXD4C0cSrPp7HoqZgtqgUhSJjuvp2OjqbFrRofXb5pqNNp5SYw727jXDQpeNMZs5dJNUovvFNh4ejlOJX7nxSzNyLqHAhgNRpIyS6n/z5SyUu3IdEaq0cxy08Ouj7TzMj0j79jSu4XGlJtsujTdQSWmNDsFHAeMBAT0DAfUXo4i2N9AcaYTditdqqYjDWg0Sa4cHRslxsiN6uiKoDGxQWvmSL5jM42q6CI0LYeVR3G0GZkRt4AnQ2MBEmemuWiDviB6owE83fRML1ZxZLxyaQpeNIaNCOgbCcLrM16G+Hi085I4szUjY/OuKHiZWbmJcuXYrktumh3p0eaNZHJGzpWj4EVj0h1W5KYZN0N8vKbu6EDGPHqDmonZBjRSpVHizNSojo6MEsf+dlr7jNf4lIIXDZLusA1+IRJFUUrYpZ2XmRk5+W4y1F03cWYa6spaK9A1Y2bFrhTDNj6l4EWDzFIu3TXkx0hAgIWD1JyPTM1MRwMAzTWaDfZG7ukbMdwd9ninu2nnJVHxLSiMdt2g4EWDzFIuzbZ/S7JSqMV3AmI9G4yXfDdeUAijbYB6vCSqKNMJO29BUBCl35tRxXZe6HWRCKPu2FLwokFGTrKKRzONZoddhM4MjCIQCqu8Gnmd6R9FWAScNgvy0x1qL0fzeAsHN6ssMfhxcyznha4biWA3PS0Gez+h4EWDzNLumxLvZicv3Y4UGw9RBFr7jX133Rw304jjOJVXow8VJtix9QUFtEdzNypoRy4htPNCFMNyG1r7RhESjHuHzZJ1aeclMRzHGfZCNF4sWZfeoBJVboKKo5a+EYhipCozJ1qVSaZn1GsGBS8aVJgR6SgaCotoGzBWhng82nmZPbP0eqFk3dkzw3TpprhkXdqRSwyr1vP0GWu3loIXDbJYOJRFM8SNehcVXyZdmUc7L4kyS68X9t9HVWiJM0OvF5ppNHvsmtE16MdowDiDOyl40aiKXGPPOOobCWLQFwLH0d31bJSzuyiDvi4YalA3e1KLhR7jDu+kjtyz50qxISM6uNNIM44oeNGoWLm0MbeA2a5LcaYTThuVSSeqNMsBX/MRvLn3T6irq4MgGOdOKh7lvMwe+10N+kPoHzHmaBHpqJleFwmLz5Uz0k0PBS8aZfRyaVYmTXdQiautrcVXrvgUOp65Fw1PfB/r169HZWUlamtr1V5aUnl9QenNl4KXxDltPIqiwzuNWnFEOy9zU5ZtvFw5Cl40yujl0myaNM00SkxtbS22bNmC9rYzYx5vbW3Fli1bDBXAsIA9hwbvzRqrVDTidOmQEEZLNOmUcl5mx4jduSl40ajYRciY59d0B5U4QRCwdevWSV8H7LFt27YZ5giJVRrRrsvsVRh4LtqZfh9CYRF2q0XaYSKJie3kG6fiiIIXjWLbfIO+EAZGjXd+Le280B3UjOrr69HS0jLl50VRhMfjQX19vYKrko9HalBHlUazZeTRIiz/rzwnFRYLlUnPBvtbMlIaAgUvGpVi51GQEWmLbsTSR9p5SVxbW1tSn6d1VGk0d+W5xt15oWTduYtvVGeUnXwKXjTMqJ0R+0cCUkImnV3PrLi4OKnP0zr2eqfgZfbYzYARqxTphmfuSrNTwHHAaFBAz3BA7eUkBQUvGmbEJCsgdgdVkOFAqp0SMmdSXV0Nt9s9ZUdRjuNQVlaG6upqhVcmDw+VSc8Z25Xo8PrhCxojB4qhjtxz57DGKtGM8n4ia/DS19eHm266CS6XCy6XCzfddBP6+/un/ZpbbrkFHMeN+bjooovkXKZmlRs0+Y5mGs0Oz/PYtWsXAEwIYNi/d+7cCZ7Xf7+ccFiU2pjTzsvsZaXGNSQzyJsUQ8HL/LijPaKeeuppQ/SIkjV4ufHGG9HQ0ICXXnoJL730EhoaGnDTTTfN+HVXXnkl2trapI8XXnhBzmVqVoXBd17oIpS4mpoa7NmzB6WlpWMeLywuwZ49e1BTU6PSypKra8iPQCgM3sKh2EUVJbPFcZz0d2WkXDlRFKWjMDo2mr3a2lq89P+2oOOZe/HQvXcaokeUbHv2R48exUsvvYQ333wTF154IQDgl7/8JdauXYvjx49j6dKlU36tw+FAUVGRXEvTDaPmvNBMo7mpqanB5s2bUV9fj395og4enwOP/svfY/PqMrWXljTstV6S5YSVp1PtuajIScMHrV5DVRx1DvrhC4Zh4YDSLKpCmw3WI2p8oi7rEaXXmx/Zrg4HDx6Ey+WSAhcAuOiii+ByufDGG29M+7V1dXUoKCjAkiVLcNttt6Gzs1OuZWpaeU7kzf3MwCgCobDKq0me01KZNAUvs8XzPNatW4fLrqqBs3wlPuk2Tt8GIL5Mmnbl5ipWcWScpF12zSjNToHdSkFtoozcI0q2V0F7ezsKCgomPF5QUID29vYpv27jxo146qmn8Nprr+EnP/kJ3n77bVx22WXw+/2TPt/v98Pr9Y75MIq8dDtSbDxEEWjtN86bVGw0AL1BzdXiwnQAwImOQZVXklxUaTR/Ruz1IlUa5dANz2wYuUfUrIOXHTt2TEioHf9x6NAhABOTC4HIL2uqqgkAuO666/C5z30OK1aswKZNm/Diiy/i448/xvPPPz/p8x944AEpIdjlcqGszDhb6PEDtYzS7tvrC0qlehS8zN2SwgwAwMcGDV6o0mjujNhll/Lk5sbIPaJmnfNy55134vrrr5/2OZWVlThy5Ag6OjomfK6rqwuFhYUJ/7zi4mJUVFTgxIkTk37+nnvuwfbt26V/e71eQwUw5bmpON4xaJjKAXZBzUu3I8NpU3k1+sWCl6aeEfhDAhxW/VcaAUALjQaYN3Zs5OkbgRAWwRugGy3bRaLgZXaM3CNq1sFLXl4e8vLyZnze2rVrMTAwgLfeegsXXHABAOBvf/sbBgYGcPHFFyf883p6euDxeKb85TocDjgcjoS/n97Edl6MEbw0UaOppCjMdCDDacWgL4TG7mEsK8pUe0lJ0UyjAeat2JUCG88hKIho9/oMkeBKDermhvWIam1tnTTvheM4uN1uXfaIki3nZfny5bjyyitx22234c0338Sbb76J2267DZ///OfHVBotW7YMzz33HABgaGgId911Fw4ePIimpibU1dVh06ZNyMvLwzXXXCPXUjXNaOXStP2bHBzHxR0dDam8mvkTBAGv/PVVfPK3l+FrPoJSl3FvSOTGWzgp4dkox8103ZgbI/eIkjVt+6mnnsI555yDDRs2YMOGDVi5ciV++9vfjnnO8ePHMTAwACDyi37//fexefNmLFmyBDfffDOWLFmCgwcPIiMjQ86lalaZwcqlWbIuVRrN3xKDJO3W1taisrISn73iM+j+04/R8cy9WLNiqa57UKjNSDOO+kcC0nBaSuSeval6RLndbt2WSQMy9nkBgJycHDz55JPTPid+KyslJQUvv/yynEvSnfEDtaZLdtYDuoNKnsUFkYD+eLt+gxej9qBQW4WBKo5onMj8sR5Ru2tfxLd/ewDZeQU4+stv6XLHhaGCeY1zRwdqjQSMMVCLRgMkDzs2OtGpz2MjI/egUFuZgSqOKFk3OXiex42bNyL97EsRKFiO3pGQ2kuaFwpeNM5h5VEcHail96TdkUAInYORfj0UvMwfOzY63TOsyyF8Ru5BoTYjTZc+3U3JusmSYuelXTk979gCFLzoglT6qPMtYBZ8Zafa4EqlMun5ys9wwJViQ1gEPunS3+6LkXtQqC1+vtFkO1t6Iu28UL5LUiwtiuzYHmvXd0NXCl50wCjl0k10B5VUkYojlrSrv+DFyD0o1MauGYO+EPpHgiqvZn6kMmmahZYUS6NtFWjnhcjOKAMam6SZRnQHlSx67rTLelBMlYTOcRzKysp02YNCbU4bj8LMSLm53q8bUpI/7bwkxbLozstxHV4z4lHwogPl0Z0K/R8b0c5Lsum510t8DwoYrAeFFrA5QHquOIrPk6OE3eRgx0YfdwxCCOv3SJGCFx2IDVrTd/KdVGmURxehZJEGNHbq8y6K9aDIzhs7MkTvPSi0wAjTpdmukSvFhqxUu8qrMYbK3DQ4rBb4gmFd78pR0bwOsO3SDq8fvqAAp02fd6KxHi+085IsbOeluXcEowEBKXb9vTZqampwwF+B//7zK/hslRP/cMVqVFdX047LPFUYIFeO+kIlH2/hsLgwHR+0enG83YsqneYS0c6LDmSl2pDhiMSZej068gUFtA34AFCZdDLlpTuQk2aHqNOKI+ZYxzCc5Svx91++EevWraPAJQnYzouej43oqFkeSwsjSbvHdJy0S8GLDnAcF9sC1umFiK07w2lFNpVJJ9XigsjRkR6TdgHAHxJwMtpob3mxOceAyKHcAI3qKFlXHsuK9Jvoz1DwohN6rziKn2mk9xEHWqPnpF0AONk5hFBYRKbTaogJyFrBdivavT5dNjEE6NhILrFeLxS8EJnpvdcLXYTko/cBjUfbIuteXpxJgW0SZRvguJkVKdCxUXKxnZembn125wYoeNENvXfZZZVGek0O07LFbOdFpxVHR9sinT6XF2eqvBJj0ftxcyAURmvfKAC66Um2/AwHslIj3blP6nQ2GgUvOlGu8ymxVGkkH3Zs5OkdxUhAf8PWWPByFgUvSRc/JkBvWvtHERYBp82CggyH2ssxFI7jsLRQ30dHFLzoBGs45ekdQViHjYUapZwXuoNKtpw0O/LSIz0w9HYXJYoi7bzIqDx63dDjzotUaZRDeXJykDrt6nTGEQUvOlGc5QRv4eAPhaWOk3rhDwk4M8C2f2nnRQ6LC/SZtNvh9aNvJCj1niDJFdt50V+jOsqTkxebcUQ7L0RWNt6CkiwnAP3dRXl6RyGKQJqdl3YISHLpNWmX7bosyEvTbfNFLdPzcTMFL/JaKu286OuawVDwoiMVOt0Cjm80Rdu/8lis0wGNH9GRkaxY8NLSO6q7OTbUoE5eLHjpHPSjbzig8mpmj4IXHSnL0eesEmmaNM00ko1ee71Qvou8SrJSYOM5BIQw2r0+tZczK2y3iHZe5JHusMKdHemrpMejIwpedKRCp2WPdAclP3Zs1No/imG/fiqOYsELddaVA2/h4M7WX95LOCxK1zm240yST89JuxS86Ihez6+lnRe6g5JNVqod+dFy0hM6qTjyBQWpCo3KpOXDrht66hHV7vUhEArDauGkXD+SfFLei852bAEKXnRFjxchgHZelMJ2Xz7WyRbw8fZBhEUgNy0WeJHk02OvF7ZWd3YKrDy9TcmFVRzRzguRFeuW2T0U0M3RQFAIoyXaJZOmScsrVi6tj+AlPt+FErnlo8cdW7rhUUZsQOMQRFFfCd0UvOhIptOGrOhEZr3kvbT2RaocnDYLCjPp7lpOUtKuTo6NKN9FGSwA0NN0aUrWVUZVXhpsPIchf0i6ydQLCl50pkJn06Ube2iatFL01uslfiAjkU9sqKt+EnZp50UZNt6ChfmR64be+r1Q8KIzsXJpfQQvp7vZRYjuoOTGer20Dfjg9QVVXs30aCyAcljw4vWF0D+ij34eUoO6HLpuyC2WtEvBC5GR3sqlY5VGdAclN1eKDUWZkcqMExqvHmjpG8WgPwQbz0l3fkQeKXZeGmyoh6RdURSpu66CWPCit14vFLzoTLnOjo1o+1dZi3VydMQ66y4uyIDdSpchuUkVRzq4bvQOBzDkD4HjYjvNRD567fVCVw2dKdNd8EI9XpSkl067dGSkrPK4qfRaxwKsokwnzbtSACuXPtU1jEAorPJqEkfBi86wHYyWvhHNzyoJCWF4+qLbv3m086IEKWm3U9s7L1RppCw9TZeO7dbSDY8SSlxOZDitCIVFfNKl7ZueeBS86ExRphM2nkNQENE2oO3StrYBH4KCCLvVguJM6pKpBL0MaGSVRtRZVxl6alQXS9alGx4lcByHpTq5bsSj4EVneAuHsmx9HB01Re+gynNSYbFQmbQSFhdEdl46vH4MjGqz4mjQF5Reu3RspAw95cpJwQsNclWMHpN2KXjRIb2US9NMI+VlOG0ocbGKI21eiFg/iaJMJ7LT7CqvxhxY8NLu9cEXFFRezfSkYyPaeVFMLGlXm9eMyVDwokN6KZdmPV6oTFpZizWetEv5LsrLSbMj3WGFKEby5bSMyqSVF5txRMELkZFetoCbeihZVw3SgEaN7rx8RJ11FcdxXFynXe1eN4b8IfQMRxrplVPwohiW89LaP6r5BpcMBS86pJdy6SZpNABdhJTEdl60WnFEZdLq0EPSLjsyykmzI9NpU3k15uFKjTW41MtUegpedEgPx0ZCWJRycujYSFla7vUihEVpa5qCF2WV6+C6QUdG6tFb0i4FLzrEqo36R4KarShp9/oQEMKw8RyKXVQmrSRWcdQ16NfcLJvTPcMYDQpw2iyoouNERbEEWF0EL9RZV3F6S9ql4EWH0hxW5KVHZpVotWMmS9Yty06FlaeXmZLSHFaUZqUA0N7uC+vvsrQwAzyVzytKD43qaJyIepZS8EKUUJ4TeXPS6l1UE23/qkqrSbuU76IelrDr6RtFWKPduenYSD2xYyMvRFGbr494sgYvP/zhD3HxxRcjNTUVWVlZCX2NKIrYsWMHSkpKkJKSgnXr1uHDDz+Uc5m6pPXKAbqDUhfLe9FarxcKXtRT7HLCauEQCIXR7vWpvZxJ0WgA9SwqSAdv4eD1hdDh9au9nBnJGrwEAgFce+21+NrXvpbw1/zoRz/CQw89hIcffhhvv/02ioqKcMUVV2BwUFsXYbW5s5zwNR/BS3/cg7q6OgiCthpPUaWRurTa6+UjCl5UY+UtcGdHdmy1eNPjDwloiwZVdNOjPIeVl/LQjulgwrSswcv3vvc9fOtb38I555yT0PNFUcTOnTvxne98BzU1NVixYgV+85vfYGRkBE8//bScS9WV2tpa/NvNl6PjmXvxp53fxvr161FZWYna2lq1lyY5TT1eVKXFY6P+kQDaBiJvTsuoQZ0qynNZ0q728l48vaMQRSDNziOXOi+rQk95L5rKeWlsbER7ezs2bNggPeZwOHDppZfijTfemPRr/H4/vF7vmA8jq62txZYtW9DT2Tbm8dbWVmzZskUTAYwoinE7LxS8qGFRtOKoZziAniFtbAGzXRd3dgr18FBJhYaPm+OPmjmOkrnVsKyQgpc5aW9vBwAUFhaOebywsFD63HgPPPAAXC6X9FFWVib7OtUiCAK2bt06aTIVe2zbtm2qHyF1DvrhC4bBWzhpm5ooK9VuRVmOtiqOjlJnXdVpuUcUJeuqT0+9XmYdvOzYsQMcx037cejQoXktanzULYrilJH4Pffcg4GBAenD4/HM62drWX19PVpaWqb8vCiK8Hg8qK+vV3BVEzVGy6Td2SmwUZm0apYUaKvTLkvWPYuCF9VoebQIJfmrjwUvJ7uGEBLCKq9metbZfsGdd96J66+/ftrnVFZWzmkxRUVFACI7MMXFxdLjnZ2dE3ZjGIfDAYfDMaefpzdtbW0zP2kWz5MLXYS0YUlRBl491qmZvBeqNFIf+5vU5LFRL+28qK0sOxWpdh4jAQFNPcNYVKDd3LRZBy95eXnIy8uTYy2oqqpCUVER9u7di9WrVwOIVCzt378f//7v/y7Lz9ST+IAuGc+TS5M0FoAuQmqKJe2qf2wUFMI4EV0H7byopywnBWJYQMfxI3h8dwcWVpahuroaPM+rvTTqrqsBFguHxYUZeM/Tj2Ptg5oOXmTd029ubkZDQwOam5shCAIaGhrQ0NCAoaHYxXTZsmV47rnnAESOi7Zt24b7778fzz33HD744APccsstSE1NxY033ijnUnWhuroabrd7yiM0juNQVha5GKmJdl60YXFBrNeL2k2nTnUNIyCEke6wUh6Uil76y5/Q9p9fRccz9+LWf/iKZioVhbCIlj6qUNQCvSTtznrnZTa++93v4je/+Y30b7absm/fPqxbtw4AcPz4cQwMDEjPufvuuzE6Ooo77rgDfX19uPDCC/HKK68gI0O7EaBSeJ7Hrl27sGXLFnAcN/YNKRrQ7Ny5U/W7qKZu2nnRgkUF6bBwQN9IEN1DAeRnqHe8yo6MlhVlwEJjAVTBKhXHB7KsUnHPnj2oqalRZW1n+kcRFETYeYs03ZioQy9Ju7LuvOzevRuiKE74YIELEEkyveWWW6R/cxyHHTt2oK2tDT6fD/v378eKFSvkXKau1NTUYM+ePSgtLR3zuCu3UNWLDyOKIu28aITTxksJmmp32qV8F3VpvVKRHRmV5aTQzCuV6WVAI5WC6FBNTQ2ampqwb98+bL//YRTecD8uuvcZ1QMXAOgeCmA4IIDjIJXqEvXEOu2qeyGizrrq0nql4uleuuHRCrbz0tw7gpFASOXVTI2CF53ieR7r1q3Dvd+4Fc7ylTjWMYzOQfXnlbBdlxJXChxW9ZMAzU5K2u1UN2k31uOFjn/VoPVKRerxoh256Q7kpUeOmLWQ7D8VCl50LjfdgRWlkbvZ1090q7yauEqjPLoIaYEWBjR2DfrRPeQHx8Xu6oiytFypKAgCDr5+AMMf7cdo0xHVm2yS+KMj7Xasp+DFAC5ZnA8AOPBxl8oroUojrWEVRx93DKlWccTyXapy05Bql7VGgExBq5WKtbW1qKysxF8euB3df/4x/v2bN2qi+sns9JC0S8GLAVRHg5fXT3YjHFa3JJbtvFRR8KIJC/LTYOGAgdEgugbVmXFEybrqY5WKwMQO5pxKlYqs+ml8Lo6W5rSZlR4GNFLwYgBrKrKRaufRPRTAUZW3+Zq62c4LHRtpgdPGS8Mx1Tq/jgUvdGSkpqkqFUtKShWvVNR69ZPZ6aHiiIIXA7BbLVi7IBcAcOBj9fJexkyTpkZTmrFY6rSrzoWIKo20I75S8dybv4vCG+7HT56tV7xSUevVT2a3uCADHBeZSq/Wju1MKHgxiEuWRI6O6k+ol/fSNxLEoC9SWldOLb41Q0raVWFAoy8o4JOuSEBLwYs2sErFm//+y3CWr8RrKtzwaL36yexS7Lw0pkGruy8UvBhE9eLIvKlDTX2q1OYLgoA9f3kZwx/tR2rPMdjolaUZsV4vyh8bnewcghAW4UqxodhFnVO15DPLI8NuD3zcDV9Q2eMZLVc/kYhY0q42K47oLcYgqvLS4M5OQUAI42+nehX92axi4J+u/wK6//xjHP3VXVQxoCFL4o6NlK44+igu32WqSheijrNLMlGU6cRoUMDBUz2K/mytVj+RmKVFkZ1S2nkhsuI4Tqo62q9gyTRVDGhfVV4aeAuHQV8IHV5lz6+p0ki7OI7DZ84qAAD89aMORX82q36aLJZWq/qJjMWSdtXuzj0VCl4M5NIlkaMjpfJeqGJAHxxWXhqSqfSFiAUvZ1Hwokns6OjVo52K78rV1NTg3H/4PviMvDGPu91uTcxpM7ulRbHjZrVbcEyGghcDWbswDxYO+KRrGK39o7L/PKoY0I8lKsw4EkUxbiwABS9adNGCXKTaebR7ffjwjLK5DSc7h9BXsBoVd/waf37xFTz99NPYt28fGhsbKXDRgMrcNDisFowGBTT3jqi9nAkoeDEQV4oN55ZlAQDqFTg6oooB/VBjQGPbgA8Do0FYLZxUrk20xWnjpQ7dexU+OvrLkTMAgEuWFuLzV16BG264AevWraOjIo3g4/5utdhpl4IXg2El0wcUODqiigH9iCXtKldxxI6MFuan05BODfvMWZGjo78eVS54EUURfzkSuan5/MoSxX4umZ2lhdpN2qXgxWCkUQEnuiHIfE5JFQP6wY6NTnYqN+OIOuvqw/ql+eA44MMzXpxR4LgZAI53DOJk5xDsvAVXnF2oyM8ksyd12u3QXrk0BS8Gs8rtQqbTCq8vhPda+mX9WfHzUsajigFtqcxNg43nMOQP4cyAT5GfSfku+pCb7sCa8mwAwKvHOhX5mc9Hd10uWZKPTKdNkZ9JZk/LAxopeDEYK2/B3y2KVh0p0DmzpqYGNf/yE6oY0Di71YKqPDbjSJkLEZVJ64d0dKRA3kv8kdGmVXSkrGUseGnqHla8keFMKHgxICVHBfiCAo6nnIXS2x/Hz377B6oY0DCWtHtCgeBlJBBCYw+NBdCLzyyP9Hs5+EkPhvzyduj+8IwXjd3DcFgtuHw5HRlpWUGGA1mpNoTFyJGzllDwYkBsVMBhTz+8vqCsP+ulD9ox6AvBnZOOr9/4BaoY0LAlBcqNCTjePghRBPLSHcjPcMj+88j8LMxPR2VuKgJCGK/LfNPDdl0uW1aAdIdV1p9F5ofjOCwt1ObREQUvBuTOTsWC/DQIYRFvnJS37ffv3vYAAK493w2Lhdq/axmrOFJi5yWW70LJunrAcZzUsG7vR/LlvUSOjCIl0lRlpA9a7bRLwYtBsd4NcpZMn+4ZxsFTPeA44Nrzy2T7OSQ5pGOjTvk7ZlJnXf1hRzj7jnfKVqn4XssAWvpGkWLjcdmyAll+BkkuNuOIdl6IIi6Jjgo48HGXbKWxe96JdNf99KI8lGalyPIzSPJU5qbCzlswEhBk78D8ESXr6s75ldlwpdjQOxzA4eY+WX7GX96L7LpcvrwAKXY6WtYDlrR7XGPTpSl4MagLq3Jh4zm09I2iqSf5rZ2FsCgFL1+iXRddsPIWLMiPVByd6JTvLiocFnGMghfdsfEWrF8a7bYrQ8O6cFjEC+9TYzq9YcFLh9eP/pGAyquJoeDFoNIcVqypiPRukKPq6MCJLrQN+JCVasMGajKlG7ExAfIl7Xr6RjAcEGCPC5aIPshZMn3Y04czAz6kO6xYFw2SiPalO6xwZ0d21rV0dETBi4FJowJkmHP0+0ORRN2rzy2l1u86sqSAjQmQ7yLE8l0WF6bDxtMlRk8uWZIPq4XDJ13DaOweTur3/vN7kV2XK84qhNNG1ww9kTrtUvBClMCSdg9+0oNAKJy079sz5JeGuNGRkb7Eer3It/PyEXXW1a1Mpw0XLcgFALyaxKMjYcyRETWm0xstdtql4MXAzirORG6aHcMBAe8mMQHvucOtCAoizil14awSeoPSE1YufVLGiiPqrKtvrGFdMqdMv93Ui85BPzKcVmn+GtEPVnGkpaRdCl4MzGLh8Olow7pk5b2Iooj/iR4ZfelTtOuiNxW5abBbLRgNCmjpk6fiiMqk9Y2VTB863Ze0BE3W2+WzZxfBbqW3Hb1ZGpcrp9Rg15nQq8jgpH4vSZpz9F7LAD7uGILDasEXVlHFgN7wFg4L8+XLe/H6glJQRMGLPpXlpGJZUQaEsIi64/O/6QkJYbz0QTsAOjLSqwX5scGuct30zBYFLwbHRgV8cGYAPUP+eX8/1lF344oiuFJoGqwesaOjj2Uolz4WzXcpcTnhSqXXh15J3XaTkPfyt8ZedA8FkJ1qk4bGEn2x8RZZb3rmgoIXgyvIdGJZUQZEEfjfT+Y3KmAkEMKfo02m6MhIv5bImLRL+S7GcHk072X/8a55J/uzI6MrVxRR9ZmOaS1pl15JJpCskukX32/HkD+E8pxUXFSVm4ylERUslrFcmoIXY1jlzkJeugND/hDeauyd8/cJCmG8KB0Z0TGzni3VWLk0BS8mwPJe6k/Mb1TA76KJuteuoSGMesZ2Xk52DiV9hg0FL8ZgsXBS1dFf53F09L8nu9E/EkRumh0XVuUka3lEBVrr9ULBiwmcX5kNh9WCDq9/zp1VG7uH8VZjLywcsOV8d5JXSJRUlpMKh9UCfyiM5t7kjY4QwiKOd9A0aaO4XJoy3THnm56/HIn0dtl4ThGsdGSka6xc+pOuoaT2DZsrejWZgNPG48Jo46m5lkyzjrqXLMlHsYuGMOoZb+GwSIajo8buYfiCYaTYeFTk0lgAvfv0ojw4rBa09o9KQelsBEJhvPwhHRkZRYnLiQynFaGwiFPd8jW5TBQFLyZxSbTqaP8c8l5CQlgawngdddQ1hFjSbvKCF3ZktLQoAzwdK+peip2XqhXnMuuo/kQXBn0hFGQ48KlKOjLSO47jpH4vWjg6ouDFJFjS7luNvfAFhVl97f6Pu9A56EdOml3aSib6tigvFb7mI/jLc3tQV1cHQZjda2IylO9iPLGS6c5Zfy07MrrqnGIKZg1CSxVHFLyYxOKCdBRlOuEPhWddPcA66l6zupS6YxpAbW0tdnx5HTqeuRcvP3wv1q9fj8rKStTW1s7r+34kddalfBejuGxZJGn3PU8/Or2+hL/OFxSk8QKbVlFjOqPQUtIuvROZBMdx0hbwbPJeugb9eDV610VDGPWvtrYWW7ZsQXdH25jHW1tbsWXLlnkFMLTzYjwFmU6sKssCALx2LPHdl7rjXRjyh1DscmJ1WbZMqyNKi804Mnjw8sMf/hAXX3wxUlNTkZWVldDX3HLLLeA4bszHRRddJOcyTSPW7yXxUQHPHW5BKCxiVVmWtGVI9EkQBGzdunXSyhH22LZt2+Z0hNQ7HECHN9LBeRkFL4ZyxRxKplljus+dU0xtFQyE5by09o/C6wuquhZZg5dAIIBrr70WX/va12b1dVdeeSXa2tqkjxdeeEGmFZrLpxflgeOA4x2D6EhgCzgyhJESdY2ivr4eLS0tU35eFEV4PB7U19fP+nuzXZfynFSkO6xzXiPRHpbnVn+iG6OBmQPb0YAg7dZ+nuafGYor1YbCdBt8zUfwH4/tTlq+3FzIepX53ve+BwDYvXv3rL7O4XCgqKhIhhWZW3aaHStLXXivZQD1J7qxZc30/Vrebe7Hyc4hOG0WOrc2gLa2tpmfBOC0p3XW3zt2ZES7c0azrCgDpVkpaO0fxf+e7MZnzpo+af+1Y50YDQooy0nBKrdLoVUSJdTW1uKDh76Gkb5O/H/PRB5zu93YtWsXampqFF2LJnNe6urqUFBQgCVLluC2225DZ+fUZ61+vx9er3fMB5la9eLERwX8T3QI41XnFCPDSUP29K64OLEAdMdfW3FP7ft4t7kv4eZkH1G+i2FxHIcrogFLIkdHsSOjEnAcHRkZBcuXG+kb+36cjHy5udBc8LJx40Y89dRTeO211/CTn/wEb7/9Ni677DL4/ZNPRH7ggQfgcrmkj7IyOt6YDst7ef1kN8LTtIYf9oekixAdGRlDdXU13G73NG8oHOyufIQLl+GZt5pR88gbuOKnB/Cf+z9B5+Dkx4yCIKCurg6v/uU5+JqPYEkBNaczosulvJfOaa8bQ/6QlNj7+ZW0W2sUcubLzdWsg5cdO3ZMSKgd/3Ho0KE5L+i6667D5z73OaxYsQKbNm3Ciy++iI8//hjPP//8pM+/5557MDAwIH14PJ45/2wzWF2ehTQ7j97hAD48M/Uu1fPvt2E4IKAqLw0X0EwSQ+B5Hrt27QKACQFM5G8XeOpXv8B//9PfoWZ1KZw2C052DuGBF49h7QOv4au738ZLH7RJrcFra2tRWVmJ9evX48hvv4+OZ+7FrRsvVPwOjMjvwqpcpDus6B7y40jrwJTPe/VoB/yhMCpzU3F2Ce3CGYWc+XJzNeuclzvvvBPXX3/9tM+prKyc63omKC4uRkVFBU6cODHp5x0OBxwOR9J+ntHZeAvWLszDX4924MCJLpwzxZk0OzK69vzp7tSJ3tTU1GDPnj3YunXrmIuR2+3Gzp07pXPrtQtz8b3NZ+P5I234/TsteOd0H1491olXj3UiJ82OxSMf4vf/vn3CnVh72xls2bIFe/bsUfwMnMjHbrXg0qX5eP5IG/76UQfOjZZPj/fn9yJ5VZ9fSUdGRpJovlyiz0uGWQcveXl5yMvLk2Mtk+rp6YHH40n4vJ7M7NIl0eDl4y58ff2iCZ//pGsIh073wcIBXzyPhjAaTU1NDTZv3oz6+nq0tbWhuLgY1dXV4Hl+zPMynDZcf0E5rr+gHCc7h7DnnRbUvtuCjoERHHn0h1NuIXMch23btmHz5s0TvifRr88sL4gEL0c7cNdnl074vNcXlHLpPk8J/oaS6Puvku/Tsua8NDc3o6GhAc3NzRAEAQ0NDWhoaMDQUGyo07Jly/Dcc88BAIaGhnDXXXfh4MGDaGpqQl1dHTZt2oS8vDxcc801ci7VVFjS7rvNfRjyhyZ8nnXUXb+0AIWZTkXXRpTB8zzWrVuHG264AevWrZsxyFhUkI5vb1yGN759GbatECAMTt0rSI0tZCK/9UsLwFs4HGsfhGeSaeR7P+xAQAhjUUG61A+EGMNM+XIcx6GsrAzV1dWKrUnW4OW73/0uVq9ejfvuuw9DQ0NYvXo1Vq9ePSYn5vjx4xgYiJyh8jyP999/H5s3b8aSJUtw8803Y8mSJTh48CAyMuiPIVkq89JQnpOKoCDizU96xnwuKITx7DuRUtkvfYoSdclYVt6CIltibeKV3EIm8stKteP8iki33FcnqTpiCf6fX1lMR0YGM1O+HADs3LlT0Z1WWYOX3bt3QxTFCR/r1q2TniOKIm655RYAQEpKCl5++WV0dnYiEAjg9OnT2L17N1UQyWCqUQF1x7vQPeRHXrpdmmtCSDwtbiETZcRKpseWy/aPBFB/IrIbR1VGxsTy5UpLS8c87na7Vclx01ypNFGGNCrgxNjt/99FE3VrznPDxtPLg0ykxS1kogzWbffNUz1j2sO//GE7QmERy4oysKiAdsmNqqamBk1NTdi3bx+efvpp7Nu3D42Njaok59O7k0mtXZgL3sKhsXtYOr/u9Pqw7zgbwkiJumRyWtxCJsqoykvDwvw0hMLimEaXfznCqoxo18XoZpsvJxcKXkwq02nDeeVZAIAD0aOj2sOtEMIizivPorsnMi2tbSET5bDxAH/9KJL30jPkxxvR3LnPr6RZRkQZFLyYGKs6qv+4OzKEMXpkdB0l6pIEaGkLmSjnM9Gjo33HuxASwnjpw3YIYRErSjNRmUcdlokyaPyriV2yJB8/efkoXn71Vfxrzzv46N0eZC9Yic/R3RNJENtCJuZxXnk2slNt6B3y4bHf/RlPvvYefCM2XLXhi2ovjZgIBS8mdvJvr6LtP29H0NuN+34TeWw4txCvnPcI3T0TQibFWziUDXyA93/5AL4e1+/nB/t/jqL/+BldO4gi6NjIpGpra/GlL12LoHdstdFQb6cqE0IJIfpQW1uLP//0rgmNCtloCLp2ECVwYqIz73XC6/XC5XJhYGAAmZk0GGwygiCgsrJyykFbHMfB7XajsbGRKkYIIRK6dhA5zeb9m3ZeTEiLE0IJIdpH1w6iFRS8mJAWJ4QSQrSPrh1EKyh4MSFq704ImQu6dhCtoODFhKi9OyFkLujaQbSCghcTovbuhJC5oGsH0QoKXkyK2rsTQuaCrh1EC6hU2uQEQUB9fT3a2tpQXFyM6upqumsihMyIrh0k2Wbz/k3BCyGEEEJUR31eCCGEEGJYFLwQQgghRFcoeCGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdIWCF0IIIYToCgUvhBBCCNEVCl4IIYQQoitWtReQbKxhsNfrVXklhBBCCEkUe99OpPG/4YKXwcFBAEBZWZnKKyGEEELIbA0ODsLlck37HMPNNgqHwzhz5gwyMjImjGyfL6/Xi7KyMng8HpqbJCP6PSuDfs/Kod+1Muj3rAy5fs+iKGJwcBAlJSWwWKbPajHczovFYoHb7Zb1Z2RmZtIfhgLo96wM+j0rh37XyqDfszLk+D3PtOPCUMIuIYQQQnSFghdCCCGE6AoFL7PgcDhw3333weFwqL0UQ6PfszLo96wc+l0rg37PytDC79lwCbuEEEIIMTbaeSGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdIWClwQ98sgjqKqqgtPpxJo1a1BfX6/2kgxnx44d4DhuzEdRUZHay9K9AwcOYNOmTSgpKQHHcfjDH/4w5vOiKGLHjh0oKSlBSkoK1q1bhw8//FCdxerYTL/nW265ZcLr+6KLLlJnsTr2wAMP4FOf+hQyMjJQUFCAq6++GsePHx/zHHpNz18iv2c1X9MUvCTgd7/7HbZt24bvfOc7OHz4MKqrq7Fx40Y0NzervTTDOfvss9HW1iZ9vP/++2ovSfeGh4exatUqPPzww5N+/kc/+hEeeughPPzww3j77bdRVFSEK664QpoTRhIz0+8ZAK688soxr+8XXnhBwRUaw/79+/H1r38db775Jvbu3YtQKIQNGzZgeHhYeg69pucvkd8zoOJrWiQzuuCCC8Tbb799zGPLli0Tv/3tb6u0ImO67777xFWrVqm9DEMDID733HPSv8PhsFhUVCT+27/9m/SYz+cTXS6X+Oijj6qwQmMY/3sWRVG8+eabxc2bN6uyHiPr7OwUAYj79+8XRZFe03IZ/3sWRXVf07TzMoNAIIB33nkHGzZsGPP4hg0b8MYbb6i0KuM6ceIESkpKUFVVheuvvx6nTp1Se0mG1tjYiPb29jGvb4fDgUsvvZRe3zKoq6tDQUEBlixZgttuuw2dnZ1qL0n3BgYGAAA5OTkA6DUtl/G/Z0at1zQFLzPo7u6GIAgoLCwc83hhYSHa29tVWpUxXXjhhXjiiSfw8ssv45e//CXa29tx8cUXo6enR+2lGRZ7DdPrW34bN27EU089hddeew0/+clP8Pbbb+Oyyy6D3+9Xe2m6JYoitm/fjk9/+tNYsWIFAHpNy2Gy3zOg7mvacFOl5cJx3Jh/i6I44TEyPxs3bpT+9znnnIO1a9di4cKF+M1vfoPt27eruDLjo9e3/K677jrpf69YsQLnn38+Kioq8Pzzz6OmpkbFlenXnXfeiSNHjuD111+f8Dl6TSfPVL9nNV/TtPMyg7y8PPA8PyFi7+zsnBDZk+RKS0vDOeecgxMnTqi9FMNi1Vz0+lZecXExKioq6PU9R9/4xjfwpz/9Cfv27YPb7ZYep9d0ck31e56Mkq9pCl5mYLfbsWbNGuzdu3fM43v37sXFF1+s0qrMwe/34+jRoyguLlZ7KYZVVVWFoqKiMa/vQCCA/fv30+tbZj09PfB4PPT6niVRFHHnnXeitrYWr732GqqqqsZ8nl7TyTHT73kySr6m6dgoAdu3b8dNN92E888/H2vXrsVjjz2G5uZm3H777WovzVDuuusubNq0CeXl5ejs7MQPfvADeL1e3HzzzWovTdeGhoZw8uRJ6d+NjY1oaGhATk4OysvLsW3bNtx///1YvHgxFi9ejPvvvx+pqam48cYbVVy1/kz3e87JycGOHTvwxS9+EcXFxWhqasK9996LvLw8XHPNNSquWn++/vWv4+mnn8Yf//hHZGRkSDssLpcLKSkp4DiOXtNJMNPveWhoSN3XtCo1Tjr085//XKyoqBDtdrt43nnnjSkXI8lx3XXXicXFxaLNZhNLSkrEmpoa8cMPP1R7Wbq3b98+EcCEj5tvvlkUxUhp6X333ScWFRWJDodDvOSSS8T3339f3UXr0HS/55GREXHDhg1ifn6+aLPZxPLycvHmm28Wm5ub1V627kz2OwYg/td//Zf0HHpNz99Mv2e1X9NcdJGEEEIIIbpAOS+EEEII0RUKXgghhBCiKxS8EEIIIURXKHghhBBCiK5Q8EIIIYQQXaHghRBCCCG6QsELIYQQQnSFghdCCCGE6AoFL4QQQgjRFQpeCCGEEKIrFLwQQgghRFcoeCGEEEKIrvz/JWsDNP8xzmoAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from pylab import plot,show\n", "\n", "# Constants\n", "N = 26\n", "C = 1.0\n", "m = 1.0\n", "k = 6.0\n", "omega = 2.0\n", "alpha = 2*k-m*omega*omega\n", "\n", "# Set up the initial values of the arrays\n", "A = np.empty([N,N],float)\n", "for i in range(N):\n", " if i>0:\n", " A[i,i-1] = -k\n", " A[i,i] = alpha\n", " if i r2 and c - r2 - 1 < mupper):\n", " dup[c - r2 - 1,r2] -= dlow[r2 - r - 1, r2] * dup[c-r-1,r]\n", " \n", " # Backsubstitution\n", " x = np.empty(N,float)\n", " for r in range(N-1,-1,-1):\n", " x[r] = v[r]\n", " for c in range(r+1,min(r + mupper + 1,N)):\n", " x[r] -= dup[c - r - 1,r] * x[c]\n", " \n", " return x\n", "\n", "# Returns band-diagonal elements d,l,u of matrix A\n", "# See also scipy.sparse.diags at https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.diags.html\n", "def get_banded(A, mlower, mupper):\n", " n = len(A[0])\n", " d = np.zeros(n, float)\n", " l = np.zeros((mlower,n), float)\n", " u = np.zeros((mupper,n), float)\n", " for r in range(n):\n", " d[r] = A[r][r]\n", " for c in range(max(0,r - mlower),r):\n", " l[r-c-1,r] = A[r][c]\n", " for c in range(r+1,min(r+mupper+1,n)):\n", " u[c-r-1,r] = A[r][c]\n", " return d, l, u" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [-2.04 0.08 -1.76 2.92]\n", "Ax = [-4. 3. 9. 7.]\n", "v = [-4. 3. 9. 7.]\n" ] } ], "source": [ "A = np.array([[ 2, 1, 0, 0 ],\n", " [ 3, 4, -5, 0 ],\n", " [ 0, -4, 3, 5 ],\n", " [ 0, 0, 1, 3 ]],float)\n", "\n", "v = np.array([ -4, 3, 9, 7 ],float)\n", "x = linsolve_banded(*(get_banded(A,1,1)), v)\n", "print('x =', x)\n", "print('Ax =', A.dot(x))\n", "print('v =', v)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "def random_banded(n, mlower, mupper):\n", " A = np.random.rand(n, n)\n", " for r in range(n):\n", " for c in range(0,r-mlower):\n", " A[r][c] = 0.\n", " for c in range(r + mupper + 1, n):\n", " A[r][c] = 0.\n", " return A" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A =\n", " -------- --------- --------- --------- --------- -------- --------- -------- -------- ---------\n", "0.339089 0.425374 0.608148 0.99282 0.559892 0 0 0 0 0\n", "0.682256 0.196732 0.047846 0.0809591 0.475877 0.871224 0 0 0 0\n", "0.51193 0.0589133 0.857931 0.817493 0.107546 0.686704 0.05595 0 0 0\n", "0.803584 0.620484 0.0260938 0.354983 0.0235425 0.666172 0.949641 0.696003 0 0\n", "0 0.923449 0.522509 0.961997 0.04002 0.964269 0.0749992 0.877955 0.242533 0\n", "0 0 0.908196 0.578138 0.756524 0.749277 0.957537 0.276867 0.108327 0.0749192\n", "0 0 0 0.739343 0.433746 0.118236 0.506705 0.897381 0.166318 0.717684\n", "0 0 0 0 0.989404 0.520105 0.465202 0.995558 0.505746 0.854393\n", "0 0 0 0 0 0.93188 0.762252 0.682169 0.775073 0.0679517\n", "0 0 0 0 0 0 0.509245 0.614805 0.142971 0.19313\n", "-------- --------- --------- --------- --------- -------- --------- -------- -------- ---------\n", " x = [ 0.75621087 -4.1587422 -1.41753296 1.64238818 2.25837649 0.08050577\n", " -1.20402935 4.17532493 -1.35788382 -5.93396431]\n", "Ax = [0.52037706 0.90776196 0.49941561 0.44271983 0.41303592 0.84240939\n", " 0.85558548 0.11629715 0.54984427 0.61370111]\n", " v = [0.52037706 0.90776196 0.49941561 0.44271983 0.41303592 0.84240939\n", " 0.85558548 0.11629715 0.54984427 0.61370111]\n" ] } ], "source": [ "n = 10\n", "mupper = 4\n", "mlower = 3\n", "A = random_banded(n, mlower, mupper)\n", "print(\"A =\\n\", tabulate(A))\n", "v = np.random.rand(n)\n", "x = linsolve_banded(*(get_banded(A,mlower,mupper)),v)\n", "print(\" x = \", x)\n", "print(\"Ax = \", A.dot(x))\n", "print(\" v = \", v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us apply it to the springs example from Section 6.1 of M. Newman *Computational Physics*" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbqdJREFUeJzt3Xd4XOWZN/7vmTNNddTrqLkbjI0xAcxGYENwMIljUEwoWQK7gXcJIbHjZXkD+b3BySawmxBiZwlhScg6hLLZGJFGdcCyxWICBgtTbGNjyRrJ6m3Upp05vz9mnjOjPpLm9PtzXbqueDSSHpTRmfs8z104URRFEEIIIYTohEXtBRBCCCGEzAYFL4QQQgjRFQpeCCGEEKIrFLwQQgghRFcoeCGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdMWq9gKSLRwO48yZM8jIyADHcWovhxBCCCEJEEURg4ODKCkpgcUy/d6K4YKXM2fOoKysTO1lEEIIIWQOPB4P3G73tM8xXPCSkZEBIPIfn5mZqfJqCCGEEJIIr9eLsrIy6X18OoYLXthRUWZmJgUvhBBCiM4kkvJBCbuEEEII0RUKXgghhBCiKxS8EEIIIURXKHghhBBCiK5Q8EIIIYQQXaHghRBCCCG6QsELIYQQQnSFghdCCCGE6IrhmtQRkkyCIKC+vh5tbW0oLi5GdXU1eJ5Xe1mEEGJqFLwQMoXa2lps3boVLS0t0mNutxu7du1CTU2NiisjhBBzo2MjQiZRW1uLLVu2jAlcAKC1tRVbtmxBbW2tSisjhOiBIAioq6vDM888g7q6OgiCoPaSDIWCF0LGEQQBW7duhSiKEz7HHtu2bRtdjAghk6qtrUVlZSXWr1+PG2+8EevXr0dlZSXd9CQRBS+EjFNfXz9hxyWeKIrweDyor69XcFVEaXTnTOaCdm2VQcELIeO0tbUl9XlEf+jOmcwF7doqh4IXQsYpLi5O6vOIvtCdM5kr2rVVDgUvhIxTXV0Nt9sNjuMm/TzHcSgrK0N1dbXCKyNyoztnMh+0a6scCl4IGYfneezatQsT374gBTQ7d+6kfi8GRHfOZD5o11Y5FLwQMomrr74G53zle+Az8sY8nldYjD179lCfF4OiO2cyH9XV1SgtdU/5edq1TR4KXgiZxEsftsNbdB6Wb/sNnn95L67e/m8ovOF+fO0XL1LgYmB050zmg+d5XP4Pd0/6Odq1TS4KXggZRxRF/MdrJwEA/1i9CFdt+Az+4Ss3wVm+EoeaB1ReHZET5TuR+WjqHsbB8CLkX30v8grHBrhut5t2bZOIghdCxnn1aCeOtnmRZufxD39XCQC4sCoHAPBRmxdeX1DF1RE5TZfvBACiCPz0pz+lO2cyqX/9y0cICGF89vNfQFtLM17+619RdPXdKLzhfux7+wMKXJKIghdC4oiiiP/YF9l1+crFlchKtQMACjOdqMxNhSgC7zT1qblEIrOamhrceM/OCflOfEYe8q++B0PFa1RaGdGy14514NVjnbBaONy36WxYrVZsuPxyXLxhM5zlK9HQ4lV7iYZCgxkJiVN/ohvvefrhtFnw1U9XjfncBVU5aOoZwZuNPVi/rEClFRIlhMo/hdLbH8ctC31YlimguLgYJzg3fvjicfzg+Y9wblkWVpVlqb1MohH+kIDv//kjAMBXP12FRQXp0ufWVGTjndN9eKe5D19cM3UyL5kdCl4IifNwNNflxgsqkJfuGPO5C6py8T+HWvBWY68aSyMKEcIi3m8dAGfh8eWrN2JJYQYA4FJRxDvNA3jpw3bc8dS7eP6bn5Z25oi5/aq+EU09IyjIcOAbly8e87nzyrMBAO+eph3bZJL12OjAgQPYtGkTSkpKwHEc/vCHP0z7/Lq6OnAcN+Hj2LFjci6TEADAm6d68FZTL+y8Bf906YIJn2d5L++3DGAkEFJ6eUQhJzuHMBIQkGbnsTA/dgfNcRx+dO1KlOekorV/FP/8P+8hHJ4qO4aYxZn+Uemm596rliPdMXZP4LyKLADA8Y5BDFK+XNLIGrwMDw9j1apVePjhh2f1dcePH0dbW5v0sXjx4pm/iJB5YhegL33KjcJM54TPu7NTUOxyIhQWcbi5X+HVEaW85+kHAJzjdoG3jK06ynTa8MiXz4PdasGrxzrxWP0pFVZItOT+F45iNCjgU5XZ2HxuyYTPF2Q4UZaTAlEEGqKvLTJ/sgYvGzduxA9+8INZZ1gXFBSgqKhI+qDMfiK3d5v78PrJblgtHG6/dOGkz+E4Ttp9+RsdHRlWQ0s/AEyZ07Ki1IUdm84GAPz45eN0jGhib3zSjb8caYOFA3Z84ewpS+zXRI+O3qGjo6TRZLXR6tWrUVxcjMsvvxz79u2b9rl+vx9er3fMByGzxXZdas4rhTs7dcrnXVCVCwB4q7FHkXUR5bGdl3PdWVM+54YLynD1uSUQwiK+8cy76B7yK7M4ohlBIYzv/SmSpPvlCytwdolryueuqYjmvdCObdJoKngpLi7GY489hmeffRa1tbVYunQpLr/8chw4cGDKr3nggQfgcrmkj7KyMgVXTIzgg9YBvHasExYOuGPdommfe0F05+Vwcz/8IRrOZzS+oIBj7YMApt55ASK7cD+85hwsKkhHh9ePbf/dAIHyX0zltwdP43jHILJTbfjnDUumfe7q6M7L4dN9lCeVJJoKXpYuXYrbbrsN5513HtauXYtHHnkEn/vc5/Dggw9O+TX33HMPBgYGpA+Px6PgiokRsF2XL6wqQWVe2rTPXZifhtw0O/yhMI60ULddo/nwzACEsIj8DAeKXRPznuKlOaz4xZfPQ4qNx+snu/GzV08otEqitq5BP36692MAwN1XLpux6mxZUQZS7TwG/SGc6BxSYomGp6ngZTIXXXQRTpyY+qLgcDiQmZk55oOQRH3cMYiXPmwHxwFfXz/9rgsQueNmuy+U62A8DZ5IQLrKnTVl/kK8xYUZuL9mBQDgZ6+dQP2JLlnXR7ThRy8dw6A/hHNKXfjS+TPv9lt5C86N7uS920x5L8mg+eDl8OHDmhiCJggC6urq8Mwzz6Curg6CQEcGRsB2XTauKMLiaD+PmVxASbuGJeW7lE2dvzDeNavduOGCMogisO2/G9A+4JNpdUQL3m3uw+/faQEAfG/z2RMq0qbC8l4oaTc5ZG1SNzQ0hJMnT0r/bmxsRENDA3JyclBeXo577rkHra2teOKJJwBEpm1WVlbi7LPPRiAQwJNPPolnn30Wzz77rJzLnFFtbS22bt2KlpYW6TG3241du3bRrAodO9U1hL8cOQMgsV0X5sJo0u47Tb0ICWFYec3fA5AEvTdDpdFU7tt0Nt7zDOCjNi++8cy7ePq2i2Cj14XhCGER9/3xQwDAljVuqQFdIqhZXXLJ+td16NAhrF69GqtXrwYAbN++HatXr8Z3v/tdAEBbWxuam5ul5wcCAdx1111YuXIlqqur8frrr+P5559XNUCora3Fli1bxgQuANDa2ootW7agtrZWpZWR+Xqk7hOEReAzywumrRQYb2lRBjKdVgwHBHzURtVtRtE3HMDpnhEAwMrSrFl9rdPG45Evn4d0hxVvN/XhwVeOy7BCorb/OeTB+60DyHBY8X+vXDarr11dngUAONU9jN7hgAyrMxdZg5d169ZBFMUJH7t37wYA7N69G3V1ddLz7777bpw8eRKjo6Po7e1FfX09rrrqKjmXOC1BELB161aI4sTscPbYtm3b6AhJhzy9I3jucCsA4M7LZtcEkbdw+FQl5b0YDdt1WZCXBleqbdZfX5mXhh9vWQkA+M/9p7D3o45kLo+orH8kgB+9FOn2vu2KJcjPcMzwFWNlpdqlmUeHKe9l3mhfcxr19fUTdlziiaIIj8eD+vp6BVdFkuHR/Z9ACIuoXpwnJdLNBuW9GM97LFl3HgMXN55TjH/4u0oAwD//TwM8vSNJWBnRgof2foy+kSCWFKbjK2sr5vQ9qFld8lDwMo22trakPo9oQ/uAD78/FAlKvzHLXReGBS9vN/VS3waDkPJd3IkfIU7mno3LsaosC15fCF9/+l2M+AOU7K9zH53x4sk3TwOIdNKdaz4Tm3NEwcv8UfAyjUSrnLRQDUUS958HPkFACOPCqhwpCJmtFaUupNp59I8E8XHnYJJXSJQmiqJUaTSfnRcAsFst+PmNq+FKseHNV19EcWk51q9fjxtvvBHr169HZWUl5crpiCiKuO9PHyAsAp9bWYyLF+bN+XuxiqP3WvoRFMLJWqIpUfAyjerqarjd7in7PXAch7KyMlRXVyu8MjJXXYN+PP23SJL4XHddAMDGW6QLEeW96F9L3yh6hgOw8RyWF8+/V5Q7OxU12S3o+sP98PaMzX2hZH99+WPDGbzd1IcUG4/vXLV8Xt9rQV46XCk2+IJhHGujm575oOBlGjzPY9euXQAwIYBh/965cycNjtSRX71+Cv5QGOeWZeHvFuXO63tdUEl5L0bBjoyWF2fCaZv/37MgCPjVj++b9HOU7K8fQ/4Q7n/hKADgzssWoSQrZV7fz2LhcF606uid03TdmA8KXmZQU1ODPXv2oLS0dMzjbrcbe/bsoT4vOtI3HMBvD0bOrb95+aKEOqhOJ77T7mQVaUQ/pCOjaYYxzgYl++tXfEPSb/30KXQMjKAyNxW3Vlcl5fuzfi/v0JDGeZG1SZ1R1NTUYPPmzXj6jy9j++46pGfn4cTj22G10q9PT/7rfxsxEhBwdkkm1i8tmPf3W1WWBTtvQdegH009I6iaYS4S0a5kVBrFo2R/fZqsISmfkYd//MG/w2Fdn5SfIU2YpqTdeaGdlwTxPI/rvnAlMlesg1h8NrqHQ2ovicyC1xfEf73RBAD4xmXz33UBIo3JWJn1W4098/5+RB0hIYz3WyPBy2zGAkyHkv31Z6qGpMJgN7637dak5SitKsuChQNa+0dplMQ8UPAyC3arBRU5qQCAkzQZVFeeeKMJg74QlhSmY8NZRUn7vlK/l1N0fq1XJzqHMBoUkO6wYkFeelK+JyX768t0DUmZZOUopTmsUlI4DWmcOwpeZmlBfuTi9kkXBS96MewP4fHXGwFEZhhZEhyklogLF1DSrt6xfJeVblfSXhuU7K8vSuconUfN6uaNgpdZYu2daedFP57622n0jQRRlZeGz68sSer3Pq88G7yFQ2v/KFr6qJuqHjUkqb/LeJTsrx9K5yjRhOn5o4zTWaLgRR8EQUB9fT1Oe1rxk7p2iLmL8bV1CxMeX5+oNIcVK0pdeM/Tj7ebeuHOTk3q9yfya0hypVE8luz/jw88gRfeOoqrLz4bj/7L39OOi8YonaPEgpcPzwzAFxSSUp5vNhS8zJIUvNCxkWZNVjFgd+UDFzwCnF+W9J93YVUO3vP0463GXlyz2p3070/kMxII4eOOSLOwucy4SgTP87j88vXYP1QAsbiQAhcNYjlKra2tk+a9cBwHt9udtBwld3YK8tId6B7y44PWAZxfObdO32ZGx0aztCA/Ug7bNejHwGhQ5dWQ8aaqGAgMdOP6L31Jlq6m1KxOvz5o9SIsAoWZDhS5nLL9HJYIfKqbbnq0KD5HaTw5cpQ4jsMamnM0LxS8zFKm04bCzMgodEra1ZbpKwbk62r6qcoccBxwqmsYXYP+pH5vIq9kN6ebysKCyE1Pc88IzbTRKJajZM8cO7tIrhwlynuZHwpe5oDyXrRJra6mrlQblhVFSh9pzpG+NLBJ0jIdGTFFmU6k2nmEwiKaeymxW6s++7kvoOifHkfhDffjl//1BPbt24fGxkZZkqulZnXN/dShew4oeJmDRaxcmoIXTVGzq+mF0qgAalanJ2znRa58F4bjOKkDM103tKuxexichUfJ8vNx6y03Yd26dbLlKJ1d4oKN59A95Iend1SWn2FkFLzMwcIC6vWiRWp2NZWa1dHOi250D/nR0jcKjgPOcSens+50FuazvJdh2X8WmZvG6P83Soz6cNp4rCiNvO7eaabrxmxR8DIHbOeFjo20Rc2upp+KJu0e7xhE/0gg6d+fJN+R6JHRwvx0ZDptsv88lux/im56NEvJ4AUA1lCzujmj4GUOWM5Lc+8IfEEaaa8VYysGlO1qmp/hwIL8NIgicKiJLkR60MCGMcqcrMsslLpz086LVrHghQWacosNaexX5OcZCQUvc5Cf4UCGw4qwCDT10IVIS1jFQGbe2KnRSnQ1lfJemmgLWA9i+S7yHxkBtPOiB+z/mwUK7bycFw1ejrV7MeSnYb+zQU3q5oDjOCwsSEeDpx+fdA5LlSZEG2pqavDiYCn+/PJruHpJKq5ftwrV1dWyNwe7oCoHz7zlobwXHRBFEe8pVGnEsKOIvpEgeocDyEmzK/JzSWJEUZTykaqSNKBzJoWZTpRmpaC1fxTvefrxd4vyZv4iAoB2XuaMyqW1rbnXD2f5Stxw442yVgzEu7AqFwDwQesA3UVpXHPvCPpHgrDzFsVuPlLtVpRmpQCg3Rct6hkOYNAXAscBFbnKjfmgfi9zQ8HLHNGYAO0Kh0XpOK8qV5ntXwAoyUqBOzsFQljEu3Qh0jQ2z+iskkzYrcpdBtnREVUqas+paC5SaVaKorOGYv1e6JoxGxS8zNFCqjjSrDavD/5QGDaeQ0mWfC3fJ3OB1O+Fjo607L1osq7c/V3Gk8qlKWlXcxqjoxuUqjRizitnSbt9CIepWV2iKHiZI7bzcqpriF5wGtMUPbcuy0mFlVf2JX4hBS+6EMt3USZZl4ntvFDwojUs30WpZF1mWXEGUmw8vL4Q7cjNAgUvc1SWnQI7b4E/FEZrP3VH1BK1LkIAcEE076XB009l9BoVFML4oFXZMmkmtvNCb1Ja09jFyqSVSdZlbLxFCqIp7yVxFLzMkZW3SNuLdHSkLWznpVLBfBemMjcV+RkOBISwVIpLtOV4+yD8oTAynVbFXyNs56W5lwY0as0phRvUxZOOjijvJWEUvMwDmxRLW33aIgUvKlyEOI6joyONiy+Rtlgm78Ysl/gBjad7aECjVghhEad71AteqOJo9ih4mQcaE6BNjSoeGwGxvBfq96JNbEdM6SMjIBLcUrM67WntG0VQEGG3WlASLWdX0urozssnXcPoG6bxIomg4GUeFlKvF80JCWE090buaNXYeQFieS/vnO6jowENYpVGSjWnG29BHo0J0JpTrNIoNw28wrtxAJCTZpeC2sMe2n1JBAUv8xDf60UUqeJIC1r7RxEKi3BYLSjKVLZMmllckI6sVBtGg4KUGEq0YcgfwsedgwCAVQpMkp4MJe1qDytdV+PIiImVTPertgY9oeBlHhbkpYPjgP5ou2+ivvikO6XzGRiLhZOmTFPei7Z80DoAUQRKXE4UqBTcSsdG3bTzohXSNGmFBjJOhvJeZoeCl3lIsfNSu286OtIGNSuN4lHSrjZJ+S4qHRkB8dOl6ZqhFWrnyQGx4KXB048QHTfPiIKXeaIxAdqiZqVRPDbn6K2mXgjUxFAzlB7GOBl2NEE7ttohTZNWcedlUX46MpxWjAYFHGsfVG0dekHByzxRxZG2qNmgLt7y4gykO6wY9IVwnC5EmtHQ3A9AnUojJn7HlnZf1DcaEHBmwAdAuWnSk7FYOKnqiPq9zIyCl3liFUdUOaANbCCj2jsvVt4ibQP/rbFH1bWQiE6vD2cGfOA44ByVknUZKpfWDnbNcKXYkJ1qU3Uta8op7yVRFLzMEzs2+oR2XlQXCIXR2hcZ1VCZp9xI+6nQkEZtea8lUvm1uCAd6Q6rqmuhAY3aIeW75KeB49RJ8mcoaTdxFLzMEzs2au0fxbA/pPJqzK25dwRhEUiz88hPd6i9nDFJu1RKrz41m9ONtzCfunNrBdv9UrNMmllV5oKFA1r6RtHp9am9HE2j4GWestPsyEmzA6C7KLXFlzuqfQcFRI4mHFYLeoYDdKyoAVpI1mUW0M6LZmglTw4AMpw2LCnMAEB5LzORNXg5cOAANm3ahJKSEnAchz/84Q8zfs3+/fuxZs0aOJ1OLFiwAI8++qicS0yKRVT6qAlaKZNmHFZeajxFR0fqCodFaeflXA0EL+zY6HTvCAIhKotVU+zYSL1k3Xh0dJQYWYOX4eFhrFq1Cg8//HBCz29sbMRVV12F6upqHD58GPfeey+++c1v4tlnn5VzmfNGYwK0oVHFwWpTieW9UNKumpp6huH1hWC3WrC0KEPt5aAw04E0Ow8hLErjLIjyRFHURHfdeBS8JEbWrLWNGzdi48aNCT//0UcfRXl5OXbu3AkAWL58OQ4dOoQHH3wQX/ziF2Va5fwtouBFExo1dhECxg5pFEVRE8dZZsSOjFaUZMLGq39aznEcqvLT8EGrF590DUnXEKKsvpEgBkaDALSzY8t2az9o9cIfEuCw8iqvSJvU/yuOc/DgQWzYsGHMY5/97Gdx6NAhBIPBSb/G7/fD6/WO+VAaJd9pg1bKpOOtLs+G1cKhbcCHlmglFFGe2sMYJ0MVR+prjA5kLHE5kWLXRpBQkZuK3DQ7AkIYH7Qq/36mF5oKXtrb21FYWDjmscLCQoRCIXR3d0/6NQ888ABcLpf0UVZWpsRSx2B3TU09w9TWWSWjAQFtrNGURu6ggEhDspXRniJ/o7wX1TRoKN+FYdOlqdeLeljgqJV8FyCyK3deBRvSSEdHU9FU8AJgwrY6KzGdarv9nnvuwcDAgPTh8XhkX+N4Ja4UpNh4BAURp+n8WhVs1yUr1YbsaPWXVlzARgVQ3osqAqEwPjoTuYPVQpk0s7CAdmzVFj/IVUso72VmmgpeioqK0N7ePuaxzs5OWK1W5ObmTvo1DocDmZmZYz6UZrFw0oWI8l7UobVKo3gXLqBmdWo61u5FQAjDlWJDRa76zQsZtvPySdcw9QFSiRbz5IBY3ss7zX302piCpoKXtWvXYu/evWMee+WVV3D++efDZlO3bfNMaFKsurRYacSsqciGhQOaekbQQY2nFBc/SVpLCdNVeWngOGBglAY0qiW+u66WrHS7YLVw6Br0U67cFGQNXoaGhtDQ0ICGhgYAkVLohoYGNDc3A4gc+XzlK1+Rnn/77bfj9OnT2L59O44ePYpf//rXePzxx3HXXXfJucykoAGN6mJ3UFrcecl02rC8KA2+5iP48SO/Rl1dHQRBUHtZptEQTdY9V+V5RuOl2HmUuCIDGtnxBVGOEBalm54FKg5knIzTxuPs0sjrlZrVTU7W4OXQoUNYvXo1Vq9eDQDYvn07Vq9eje9+97sAgLa2NimQAYCqqiq88MILqKurw7nnnot//dd/xc9+9jNNl0kzNONIXSznpUpjd1AAUFtbi9f/9Xp0PHMvfnLP17F+/XpUVlaitrZW7aWZgpY6645HAxrVc6Z/FIFQGDaeQ2l2itrLmYCGNE5P1j4v69atm/a8bvfu3RMeu/TSS/Huu+/KuCp5LCoYe36tpe1pM2jsjiRKa6nSCIgELlu2bJnwd9Da2ootW7Zgz549qKmpUWl1xuf1BaWj3JUaStZlFuano/5EN42PUAE7MqrITQNv0d71+ryKLPz6f2nnZSqaynnRM/YHMOQPocPrV3s5pjLoC6J7KPI718I0aUYQBGzdunXSAJ49tm3bNjpCktEHLQMQRaA0KwX5GeoP6xxvIe28qKZRQzONJsMqjo62DdLQ30lQ8JIkdqsFFTmRN07Ke1FWU3TXJS/dgQyndhK76+vr0dLSMuXnRVGEx+NBfX29gqsyl4bokZGW+rvEiyX6086L0qRp0ho8agaAYlcKSlxOCGFROvokMRS8JFFsxtGgyisxl1ilkXZ2XYBITlcyn0dmL1ZppK1kXYY1R2umAY2K09I06alQs7qpUfCSRNKMI9oCVpRWe7wUFxcn9Xlk9qSxABrMdwHGD2ik3RclaW2a9GRYv5d3m/vVXYgGUfCSRNIWcCddhJTELkJa2/6trq6G2+2eMnmb4ziUlZWhurpa4ZWZQ/uAD+1eHywcsKJUmzsvHMdJb550dKQcX1BAa3+kf4oWe0MxayqyIYYFHNhfh6eeepraLMSh4CWJaOdFHVLworGdF57nsWvXLgATx1uwf+/cuRM8r42BcEbD8gSWFGYgzSFrYeW8xMqlKXhRyumeEYgikOG0Ildj40TiHf/bX3Hm0a/ik9134+///svUZiEOBS9JxCoHugb90ph1Ij8tTpNmampqsGfPHpSWlo55vNTtpjJpmb2nwWGMk6Hu3Mpj06QX5KVptq1FbW0trv/SlxAaHDuUmLVZMHsAQ8FLEmU4bSjMjJRj0oVIGX3DAfSPRAJFreW8MDU1NWhqasJfX30NhZv/BYU33I//PfwRBS4y03JzunjUqE55pzSe70JtFmZGwUuSSUdHVC6tCFZpVOxyIsWu3eMXnudx+WXrsfzTV8FZvhJnBqgXkJzCYRFHNJ6sy8SXS9MQPmWc0uhARobaLMyMgpckW5RPYwKUpNVKo6mURXsBNfeOqLwSYzvVPYRBfwhOmwVLCrV5d83QgEblSXlyGg1eqM3CzCh4STLaeVEWuwhpMd9lMuU5kRkqHgpeZMWGMZ5T6oKV1/ZlzmnjUZoVeV1QxZEytDpNmqE2CzPT9l+1DlHynbK03uJ7vHLaeZGVIAioq6vDU08/DV/zEZxTnKH2khLCci8o70V+/SMBaYdLqzu21GZhZhS8JBnbeWnuHYEvaN5kKqVoudJoMhS8yKe2thaVlZVYv349/vDQ/0XHM/fiZ//0WV1UZbDgmyWSEvmw33FRplOzJfTUZmFmFLwkWX6GAxlOK8Ji7I2VyEMURTR2aXM0wFRYzgsdGyUXm949Psmxt7NdF2WlbLQI5crJr1HjybrMVG0W3NRmAQAFL0nHcRzlvSika8iP4YAACxcLCrSOrbN7KECTYpPECGWlC2nnRTFaz3eJx9osfOunTyJv07/gi//vMTQ2Npo+cAEoeJEFjQlQBpsmXZqdAodVH9unmU4bslIjk689fbT7kgxGKCtdWEADGpVyKtqgTus7LwzP8/jM5Zch7axLES46y9RHRfEoeJEBjQlQht7KpBkp76WHgpdkMEJZaUEGDWhUCuvxooedF6YsO3rc3Deq8kq0g4IXGbBeL3RsJC89jLSfDPV6SS4jlJVyHBfLe6FyadmEw6KUi7ggT9v9f+KxG56uQT9GA9o9/lQSBS8yYDsvp7qGIISpY6ZcmnTW44Upp6TdpDJKWSkLwqnNgnzavD74gmFYLRzc2SlqLydhrlQbMp2Ryig6bo6g4EUG7uwU2HkL/KEwzvTTNp9c9FYmzUjBC20BJ4VRykpjvV5o50UurNKoPDdV880Lx6NKxbH09f+eTlh5i5QMRkdH8giHxViLb73mvNBFKGlYWWlRccmYx/VUVkoNLuUXP01ab+i6MZY2O/QYwKKCdBzvGMTJziGsX1ag9nIMp93rgz+kv+1fYOyxUTgswmKZ/LiDzE5NTQ3yV/wdrt3xa+RwI3jktstRXV2t+R0XJjZdOjKgcapjMDJ3Wp8mPR0KXsai4EUmC/Np50VOLN+lPEd/27/FLid4Cwd/KIyuIT8KM51qL8kwzgz44SxfifOX5GPdugvUXs6sxA9o7BkOIC/dofaSDEfr06Sn46ZjozH0ddXXkVjlAAUvcjil02RdIHKsyAbx0V1Ucnl6I3lEZTrbjQPGDmikvBd5aH2a9HRiO7aUKwdQ8CKb+F4vk3X+JPPTpOOLEEC9XuTCgsFynXRcHm8hDWiUjT8koCVaqaOnHi9M/LERvadQ8CKbBXnp4DigfySyBUySS6+VRgz1epEH+33qZVzEeOxNlXZsk6+5ZwRhEUh3WJGvwyO5kiwnOA4YDQroHqL3FApeZJJij20B07C15Dul00ojpiwn8tqg8+vkYnfWet15oXJp+ZyK263VYzK0w8qjOJofR71eKHiRFY0JkEdICEtv+lU63P4FqHJADsP+kHRHytqp681C2nmRjZ7zXRjq9RJDwYuMaEyAPM70+xAURDisFulORG8oeEm+lmjTv0ynFa7o8Eu9YTkvnr5RGtCYZI06nGk0XhnlykkoeJGRtPNCwUtSNUbzXSpyU3XbI4UFL52DfviCNKskGaRk3Vx97roAkQGN6Q4rDWiUgd6mSU8m1p2bghcKXmS0sIDOr+XQ2KX/i5ArxYaM6KySFroQJQXbStfrkREQGWewQOoRRdeNZGrs1t9AxvFoxzaGghcZsWOj1v5RDPtDKq/GOJqiW6Z6rTQCIm9SdCFKLr2XSTNSuXQ37dgmy8BoUMqH0mueHBCf6E+9Xih4kVF2mh25aXYAtPuSTHqdaTQe9XpJLraD5dZ58CJNl6adl6Rh1wx2LKdXLOelbYByoih4kdlCqeJoUOWVGEejjrvrxovtvNBdVDIYZedlAe28JF2jAfJdACA/3QGnzYKwCJzpN/d1g4IXmUmTYukuKikCoXCsS6bOL0TUqC55RFHU9WiAeAsL2M4LdedOFiNUGgGR42aW02X26wYFLzKjiqPk8vRFumSm2XnkZ+ivS2a8curZkDTdQwGMBgVwHFCq8+ClMjcyoNHrC1F37iQ5ZYBkXYYqjiIoeJEZNapLLnYHVZGrzy6Z8WhWSfKwC3lRphMOK6/yaubHaePhzqYBjcmk52nS49GObQQFLzJjwUtT9zCCgrkTrJKBzTTSc8UAU5KVQrNKksSj85lG47EdAuq0O3+iKMaS/A1w3aAd2wgKXmRWnOlEio1HKCyaPlJOBqNUGgGA3WpBiStyh02vjfkxQo+XeDRdOnk6vH6MBgXwFk73ydxA/IgAStiV3SOPPIKqqio4nU6sWbMG9fX1Uz63rq4OHMdN+Dh27JgSS006i4WTEvAo72X+jFJpxNCAxuQwSqURE5suTcdG88UCwPKcVNh4/d+vU3+oCNn/n/zd736Hbdu24Tvf+Q4OHz6M6upqbNy4Ec3NzdN+3fHjx9HW1iZ9LF68WO6lyoZmHCVPkwGGq8WjC1FySJVGOfpO1mVY8EI7L/N3ymDXDPYaHxgNYmA0qPJq1CN78PLQQw/hq1/9Km699VYsX74cO3fuRFlZGX7xi19M+3UFBQUoKiqSPnhev0l4Urk0XYjmxRcUcGbAB8A4FyIKXpLDaDsv7IanuXcE/hDNvpoPI0yTjpdqtyIvPdL81Mw7trIGL4FAAO+88w42bNgw5vENGzbgjTfemPZrV69ejeLiYlx++eXYt2+fnMuUHUva/YR2XuaFJetmOq3I1unU4PFoxP38BYUw2gbYzosxgpf8aCfYsEgdmOdLmmlkgGRdhq4bMgcv3d3dEAQBhYWFYx4vLCxEe3v7pF9TXFyMxx57DM8++yxqa2uxdOlSXH755Thw4MCkz/f7/fB6vWM+tEYKXrqGqSR2HuKPjPReJs1Q5cD8nekfRVgEHFYL8tP13fuH4TgOCynvJSlOGWCQ63i0YwsoMuRh/BuNKIpTvvksXboUS5culf69du1aeDwePPjgg7jkkksmPP+BBx7A9773veQuOMkqctPAWzgM+UNo9/pQ7DLGubzSGrsjf6hGvAi1eX3whwTd9yhRQ3NcmbTFYoygFoiMCXivZYCOm+chEArD0xfZlTNCgzqGuuzKvPOSl5cHnucn7LJ0dnZO2I2ZzkUXXYQTJ05M+rl77rkHAwMD0ofH45nXmuVgt1pQEX2TojEBc9dksEojAMhJsyPNzkMUgdY+c5c+zpVRxgKMt1BK2qVrxlx5+kYghEWk2nkUZhpjVw6I77Jr3muGrMGL3W7HmjVrsHfv3jGP7927FxdffHHC3+fw4cMoLi6e9HMOhwOZmZljPrRIGtDYSQMa58poiXdAdFYJbQHPi9GSdZkFlOg/b/GddY1y1AxQzgugwLHR9u3bcdNNN+H888/H2rVr8dhjj6G5uRm33347gMjOSWtrK5544gkAwM6dO1FZWYmzzz4bgUAATz75JJ599lk8++yzci9VVosK0rH3ow4aEzAPjT3GC16AyJvusfZBU1+I5oONBjBKsi4TXy493VE7mZpRpkmPx8qlW6I7S7yBjksTJXvwct1116Gnpwff//730dbWhhUrVuCFF15ARUUFAKCtrW1Mz5dAIIC77roLra2tSElJwdlnn43nn38eV111ldxLlRX1epmfIX8IXYN+AMY6NgJoVsl8GW00ADN+QGOeQZKRlSRVGhnsmlHsSoHVwiEoiOjw+lCSZawj00QokrB7xx134I477pj0c7t37x7z77vvvht33323AqtS1sK4iiMyeyzfJTfNjkynMcqkGaocmB+jjQZg2IBGT+8oPukcouBlDtixETuCMwrewsGdnYKmnhE0946YMnjRf69knWDJd12DflN3RZwrI+a7MLHgxbzJd3M16AuibyTy92SU7rrxpBlH3XTTMxdG664bz+w7thS8KCTDaUNRphMAHR3NhRErjZj45DvqAzQ7rNIoO9WGDIPtyAFx06XpmjFrg76gdNRshGnS47HrRgsFL0Ru1Gl37oy88+KOlvgO+UPSLgJJjFErjRg21JV2XmavKdoXKi/dYbijZoCOmyl4UVCsYyYFL7Nl1EojIJLbwHblqOJodlqilUZugwYv0s4LXTNm7VS00shoybqM2RvVUfCioEUFVHE0V9KxUa4xL0Rmv4uaK8PvvERveDw0oHHW4nu8GJHZG9VR8KIgqVEd3UXNSv9IQDpOqcwz5puU2ZPv5sqolUZMfoYDGTSgcU6ko2YD5rsAseCla9CP0YD5AlsKXhTEdl48vSPwBc33YpsrdhEqzHQg1a5Idb/iaEDj3Bh954XjOKlZHR0dzY5Re7wwrlQbMpyR6yFr1GgmFLwoKD/dgQxn5C6qqYcS8BLVZOB8F6Y8N5K0SzsviQuHRbREt8yNGrwAsXJp6hGVOFEUpWnSCwy68wKY+6aHghcFcRyHhXkp8DUfweO/eRJ1dXUQBNqBmYkRp0mPRzkvs9c15Ic/FIaFA4qznGovRza08zJ7XYN+DAcEWDjjdV6OZ+brhjH34DWqtrYWr+64A0M9HXjwGeBBAG63G7t27UJNTY3ay9OsRoMn6wKxnI0z/aMICmHYeLqvmAm7YJdkpRj69yU1qqOdl4Sx0vKynFQ4rLzKq5GPmXPljPsXrzG1tbXYsmULhno6xjze2tqKLVu2oLa2VqWVaV+TgXu8MPkZDjisFoTFSABDZmb0ZF1mgRS8DFETwwQZvdKIiTW4NN81g4IXBQiCgK1bt0564WGPbdu2jY6QJiGKoimCF47jTL0FPBdGT9ZlKnJTpQGN3UMBtZejC0adJj0e5bwQWdXX16OlpWXKz4uiCI/Hg/r6egVXpQ/dQwEM+kPgDH52DZj7/Hou2N2mEWcaxXPaeGl36RTlvSTE6JVGTFl2LNHfbLtyFLwooK2tLanPMxNWaVSalQKnzbhn14C5z6/nQjo2MnhQC8Qn7VLeSyJYzovRpkmPV5qdAo4DRoOC6XblKHhRQHFxcVKfZyZGnmk0Xrk0aM1859dzwXpbmCF4WRiX90KmFxTCUkM/o183HFYexWy0iMl6vVDwooDq6mq43W5wHDfp5zmOQ1lZGaqrqxVemfaZodKIoWOjxPlDAtq9PgDGz3kBqFx6Nlr6RhEKi3DaLNLMMCNzmzTvhYIXBfA8j127dgHAhACG/Xvnzp3geWMfi8yFNNPI4HdQAFCeS8FLolr7RiGKQIqNR26aXe3lyE7aeaHp0jOKJeumw2KZ/IbRSKSbHpONj6DgRSE1NTXYs2cPSktLxzzudruxZ88e6vMyBbMk3gGxkt+B0SAGorOcyOTiK42m2tE0kgX5aRDDAj4+/CaeePIpanA5DVYmbYZrBhA/oNFcwQs1qVNQTU0NNm/ejIee+CN+9NybWFRZhtd/egftuEwhHBalhF0z7Lyk2HnkZzjQNeiHp28ErlSX2kvSLDZJ1+iVRkz9K8/jzH/ejpC3Gzc/E3mMGlxO7pSJ8uSA2N+A2XZsaedFYTzPY9PGzyDtrEsxkruUApdpdAz64AuGwVs4uLPN8SZFeS+JMVOlUW1tLa699lqEvN1jHqcGl5NrNEmDOqbcpI3qKHhRAXux0fHA9NiRUXlOqqHbv8eL79tApuYxSYM6anA5e9JRs4EHMsZjAXzbwCgCobDKq1GOOd4RNCbVbkVeugMAvUlNJ1ZpZOw3qHi085KYZpOMBqAGl7Mz7A9JVWhm2XnJT3fAaTPfaBEKXlRSQZUlMzJTpRFTZtKyx9mSdl4MHthSg8vZYTc8OWl2ZKUavwoNiLbayDbf+wkFLyphd9ine6n0cSqN3ZE/RLNUDQC085KIgZEgvL4QABg+F4oaXM6OmaoT45mx4oiCF5XQHfbMzFRpxLCdhNa+UQhhc80qSRS7QOelO5BqN3bBJDW4nB2zTJMez4yjRSh4UUkF23kxWWOhRAlhUWq6ZIbuukxhhhN23oJQWETbgHnOr2dDyncxQZk0NbhMnCAIqD9Qh+GP9iPY+oGpkpjNeDNMwYtKqJvq9M70jyIghGG3WlCSZfw3KcZi4eA2ad+GRDWbpNKIoQaXM6utrUVlZSWe/df/g+4//xi7vvX3qKysNE0ZuRnLpSl4UQnbeTnTP4qgYJ7ytkSxs+uKnFTwJmjxHa/chHdRs+ExSaVRvJqaGjQ1NeHq7/wn8jb9C+7a9RQaGxspcEEkcNmyZcuEqiwz9cExY64cBS8qyc+Ilbe19pknWk5UowkrjRgzXohmw2w7LwzP87i4+hKknXUpnOUr6agI1AeHYYnrA6NBDIyao3cYBS8q4TiO3qSm0WiyFt/xYq8LCmon0xIN9t0myHkZryIn8vdwmq4ZAKgPDpPmsCIvPVIabpYdWwpeVBQrlzbHi202WKWRGYMXM1YOJEoIi2jpM+fOCxCXK9dDLRYA6oMTz2xJuxS8qMhsL7bZiHXXNWHwkk2vi6l0eH0ICiKsFg7FLhPuvESDl5a+UYQoV4764MQxW6M6Cl5UxJJ2m6lceoygEJaOBsy58xJ5U+4dDmDQZ47z60SxC3NpdorpErmBaCm9lZXS+9RejuqoD06M2RrVUfCiIrYFTMdGMYIg4Pd/fhneD+oQbv0QeWnGbkI2mQynDTlp7Pya8l7imbHSKJ7FwknDO6lH1Ng+OOOZrQ+O2XLlKHhRUXk0+c7TOzJptrzZsF4NX77mKnT/+cfwPPl/UVVVZYpSx/Eo72VyUvBiwnwXpiKXJe1S3gsQKSP//e9/D2tm3pjHzdYHhyWwm+W42Xy3tRrCytuG/CH0DgeQG500bUasV8P4II71ajDTRQiI3EW95+k3zYUoUZ7ocaIZk3UZqlKc6LIrN6Hkn2zwt3yIf99YjoqyUlRXV5tix4Vhr4uWvhEIYdHwx6q086Iip41HUaYTgLkvRNSrYaJy6rI7KTONBpiKNJGejo0kp3tHwFl4VKz4FG6+6ctYt26dqQIXACh2pcBq4RAURHR4jZ8PRcGLymhMAPVqmIzZku8S5TFpg7p4LHihnJeY0z2sI7f5EvwZ3sKhNNs8Nz0UvKisnCqOqFfDJCjnZSJfUEDnoB+AeRN2gViuXDPlyklYIMcCO7My05EiBS8qqzDRi20q1KthIun8uncU4TC9QQGQmtNlOKzISrWpvBr1uLNTwHGxXDlCwQtTJl03jP9+okjw8sgjj6CqqgpOpxNr1qyZcft///79WLNmDZxOJxYsWIBHH31UiWWqgsqlqVfDZNj5dUAIo2PQ+OfXiWABvjsndcrXihnE58qZ+boRTzo2MmFTy3hmalQne/Dyu9/9Dtu2bcN3vvMdHD58GNXV1di4cSOam5snfX5jYyOuuuoqVFdX4/Dhw7j33nvxzW9+E88++6zcS1UFTRAe16th3JuS2Xo1MLyFk6rRzHykGI/1vCk3cbIuQ8fNY7Egzuw7L3RslEQPPfQQvvrVr+LWW2/F8uXLsXPnTpSVleEXv/jFpM9/9NFHUV5ejp07d2L58uW49dZb8Y//+I948MEH5V6qKtiLrd3rgy9onmqa8WpqarBnzx5k5xWOedxsvRriUd7LWM0mb1AXj5J2Y4b9IXRFc6HMnLALxCf6G79Rnax9XgKBAN555x18+9vfHvP4hg0b8MYbb0z6NQcPHsSGDRvGPPbZz34Wjz/+OILBIGw2Y51156TZkWbnMRwQ0NI3ikUF6WovSTU1NTV4R1yIx599AZeV2/F/rlxjul4N8Wj21VgseCk3+d01EDseocA29jvISrXBZeJcKCDWQqBr0I/RgIAUu3GvnbLuvHR3d0MQBBQWjr2bLiwsRHt7+6Rf097ePunzQ6EQuru7Jzzf7/fD6/WO+dATjuNQnhvrtGt2zf0+OMtX4uovfsmUvRrimWkLOBFmHw0QL/baoC67UrKuicvnGVeKDRnOyJ6E0dssKJKwOz65ThTFaRPuJnv+ZI8DwAMPPACXyyV9lJWVJWHFymJn+KdpzD1VDcSh4CVGFEUaDRCHjo1iKFk3huM40+RRyhq85OXlgef5CbssnZ2dE3ZXmKKiokmfb7VakZubO+H599xzDwYGBqQPj8eTvP8AhcS2gI1/TjmdQCgslcNWmnCa9HhmG7Q2nb6RIIYDkZwwlshsZiy3ozN6PGBmlKw7llkqjmQNXux2O9asWYO9e/eOeXzv3r24+OKLJ/2atWvXTnj+K6+8gvPPP3/SfBeHw4HMzMwxH3pTRlvAAIDW/lGERcBps6Agw7xznhj2uuge8mMkEFJ5NepiF+LCTAecNvMeJTKuVBtcKZHrodHfpGZCOy9jmaVru+zHRtu3b8evfvUr/PrXv8bRo0fxrW99C83Nzbj99tsBRHZOvvKVr0jPv/3223H69Gls374dR48exa9//Ws8/vjjuOuuu+ReqmroeCCiKXoRqsxNM3UfD8aVEnuDajFB9cB0aCzAROx3YfbjZjpqHiuW6G/sa4bsU6Wvu+469PT04Pvf/z7a2tqwYsUKvPDCC6ioqAAQafke3/OlqqoKL7zwAr71rW/h5z//OUpKSvCzn/0MX/ziF+Veqmriu+zOlA9kZKe72R0UXYSY8pxUvN86gOaeESwpzFB7OaqhMumJynOjrw0T3/QEQmGc6Y+8SdN1I6Iseqxq9JwX2YMXALjjjjtwxx13TPq53bt3T3js0ksvxbvvvivzqrSjJCsFFg7wBcPoGvSjINo902yaeijfZTwpeDH4hWgmLBeKknVjKnIoabelbwRhEUi188hPp6NmYOxOvpFvhmm2kQbYrRaUZJlnGuhUTscdG5EIalQX0UyVRhNUmCS3YTqn444TjfomPVul0dlXo0EBPQaefUXBi0aU010UnV1PwixljzOJjQag1wYTP13arOioeSKHlUdxdPfeyK8NCl40wuxJuyEhLP23085LjNlfF0DktdEazWsoo7lGEvaG3dI3AsGkk8djZdJ0zYjnNsFNDwUvGmGW8rapnOn3IRQWYbdapIm5JPZmzc6vzahtwAchLMLOW1CYQa8NpijTCTtvQVAQpaRVs6Hd2smZYceWgheNMPsdNiuTrshJhcVCZ9cMS+b2h8LS8DmzYRdgd3YKvTbiWCwc3DnmzpWTeryYfCDjeGZ4P6HgRSMqTH5+TY2mJmfjKZnbQ5VGUzJzxZEQFqVcKNp5GavMBEEtBS8awSLlrkFzdlOVyqTpIjSBGe6iphOrNKJ8l/FYsH/ahN25270+BIQwbDwnBfgkotwEjeooeNGI+HbfRn7BTUXaeaEeLxOYPXihSqOpSa8NE+68sGuGOzsVPB0njsF2KdsGRhEIhVVejTwoeNEQM7f7pp2XqZml3fdUqLvu1Mzc64WSdaeWn+6A02ZBWIRhk7kpeNEQs95hC2FRunOkMumJzFA5MB3qrjs1KXjpMV81mhS80OtiAo7jDD9dmoIXDWHl0mZ7k6Kz6+mZNagFgGF/CN1DkS6hFLxM5M5OBccBg/4Q+kaCai9HUZTkPz3ppqfPmNcNCl40RDo2MtmbFOuSWZZDZ9eTYa+Ldq8PvqCg8mqUxS688RO2SYzTxkt9kcx23EzHRtMz+mgRCl40pMLgL7apNNJMo2llpdqQ4YjMUG3pM+b59VQoWXdmRn+TmowoirTzMoMygx83U/CiIezF1tI7aqp233QHNT2O4wx/IZoKlUnPzIy9XnqGAxgOCOA4em1Mxejl0hS8aEixywmrhUNACKPD61N7OYpp6qadl5mYoenUZDy9lKw7Exb0myl4Yf+txZlOOKy8yqvRJqNfMyh40RArb4E7O/KCM+OFiHZepmbWpF0PlUnPqDyXdec2T84LHRnNjP3NDIwGMTBqvGRuCl40xmzHA+GwKHUHpZ2XqZk2eIkm7FLOy9TMmCtHNzwzS3NYkZduB2DM9xMKXjTGbE2nOgf98AXD4C0cSrPp7HoqZgtqgUhSJjuvp2OjqbFrRofXb5pqNNp5SYw727jXDQpeNMZs5dJNUovvFNh4ejlOJX7nxSzNyLqHAhgNRpIyS6n/z5SyUu3IdEaq0cxy08Ouj7TzMj0j79jSu4XGlJtsujTdQSWmNDsFHAeMBAT0DAfUXo4i2N9AcaYTditdqqYjDWg0Sa4cHRslxsiN6uiKoDGxQWvmSL5jM42q6CI0LYeVR3G0GZkRt4AnQ2MBEmemuWiDviB6owE83fRML1ZxZLxyaQpeNIaNCOgbCcLrM16G+Hi085I4szUjY/OuKHiZWbmJcuXYrktumh3p0eaNZHJGzpWj4EVj0h1W5KYZN0N8vKbu6EDGPHqDmonZBjRSpVHizNSojo6MEsf+dlr7jNf4lIIXDZLusA1+IRJFUUrYpZ2XmRk5+W4y1F03cWYa6spaK9A1Y2bFrhTDNj6l4EWDzFIu3TXkx0hAgIWD1JyPTM1MRwMAzTWaDfZG7ukbMdwd9ninu2nnJVHxLSiMdt2g4EWDzFIuzbZ/S7JSqMV3AmI9G4yXfDdeUAijbYB6vCSqKNMJO29BUBCl35tRxXZe6HWRCKPu2FLwokFGTrKKRzONZoddhM4MjCIQCqu8Gnmd6R9FWAScNgvy0x1qL0fzeAsHN6ssMfhxcyznha4biWA3PS0Gez+h4EWDzNLumxLvZicv3Y4UGw9RBFr7jX133Rw304jjOJVXow8VJtix9QUFtEdzNypoRy4htPNCFMNyG1r7RhESjHuHzZJ1aeclMRzHGfZCNF4sWZfeoBJVboKKo5a+EYhipCozJ1qVSaZn1GsGBS8aVJgR6SgaCotoGzBWhng82nmZPbP0eqFk3dkzw3TpprhkXdqRSwyr1vP0GWu3loIXDbJYOJRFM8SNehcVXyZdmUc7L4kyS68X9t9HVWiJM0OvF5ppNHvsmtE16MdowDiDOyl40aiKXGPPOOobCWLQFwLH0d31bJSzuyiDvi4YalA3e1KLhR7jDu+kjtyz50qxISM6uNNIM44oeNGoWLm0MbeA2a5LcaYTThuVSSeqNMsBX/MRvLn3T6irq4MgGOdOKh7lvMwe+10N+kPoHzHmaBHpqJleFwmLz5Uz0k0PBS8aZfRyaVYmTXdQiautrcVXrvgUOp65Fw1PfB/r169HZWUlamtr1V5aUnl9QenNl4KXxDltPIqiwzuNWnFEOy9zU5ZtvFw5Cl40yujl0myaNM00SkxtbS22bNmC9rYzYx5vbW3Fli1bDBXAsIA9hwbvzRqrVDTidOmQEEZLNOmUcl5mx4jduSl40ajYRciY59d0B5U4QRCwdevWSV8H7LFt27YZ5giJVRrRrsvsVRh4LtqZfh9CYRF2q0XaYSKJie3kG6fiiIIXjWLbfIO+EAZGjXd+Le280B3UjOrr69HS0jLl50VRhMfjQX19vYKrko9HalBHlUazZeTRIiz/rzwnFRYLlUnPBvtbMlIaAgUvGpVi51GQEWmLbsTSR9p5SVxbW1tSn6d1VGk0d+W5xt15oWTduYtvVGeUnXwKXjTMqJ0R+0cCUkImnV3PrLi4OKnP0zr2eqfgZfbYzYARqxTphmfuSrNTwHHAaFBAz3BA7eUkBQUvGmbEJCsgdgdVkOFAqp0SMmdSXV0Nt9s9ZUdRjuNQVlaG6upqhVcmDw+VSc8Z25Xo8PrhCxojB4qhjtxz57DGKtGM8n4ia/DS19eHm266CS6XCy6XCzfddBP6+/un/ZpbbrkFHMeN+bjooovkXKZmlRs0+Y5mGs0Oz/PYtWsXAEwIYNi/d+7cCZ7Xf7+ccFiU2pjTzsvsZaXGNSQzyJsUQ8HL/LijPaKeeuppQ/SIkjV4ufHGG9HQ0ICXXnoJL730EhoaGnDTTTfN+HVXXnkl2trapI8XXnhBzmVqVoXBd17oIpS4mpoa7NmzB6WlpWMeLywuwZ49e1BTU6PSypKra8iPQCgM3sKh2EUVJbPFcZz0d2WkXDlRFKWjMDo2mr3a2lq89P+2oOOZe/HQvXcaokeUbHv2R48exUsvvYQ333wTF154IQDgl7/8JdauXYvjx49j6dKlU36tw+FAUVGRXEvTDaPmvNBMo7mpqanB5s2bUV9fj395og4enwOP/svfY/PqMrWXljTstV6S5YSVp1PtuajIScMHrV5DVRx1DvrhC4Zh4YDSLKpCmw3WI2p8oi7rEaXXmx/Zrg4HDx6Ey+WSAhcAuOiii+ByufDGG29M+7V1dXUoKCjAkiVLcNttt6Gzs1OuZWpaeU7kzf3MwCgCobDKq0me01KZNAUvs8XzPNatW4fLrqqBs3wlPuk2Tt8GIL5Mmnbl5ipWcWScpF12zSjNToHdSkFtoozcI0q2V0F7ezsKCgomPF5QUID29vYpv27jxo146qmn8Nprr+EnP/kJ3n77bVx22WXw+/2TPt/v98Pr9Y75MIq8dDtSbDxEEWjtN86bVGw0AL1BzdXiwnQAwImOQZVXklxUaTR/Ruz1IlUa5dANz2wYuUfUrIOXHTt2TEioHf9x6NAhABOTC4HIL2uqqgkAuO666/C5z30OK1aswKZNm/Diiy/i448/xvPPPz/p8x944AEpIdjlcqGszDhb6PEDtYzS7tvrC0qlehS8zN2SwgwAwMcGDV6o0mjujNhll/Lk5sbIPaJmnfNy55134vrrr5/2OZWVlThy5Ag6OjomfK6rqwuFhYUJ/7zi4mJUVFTgxIkTk37+nnvuwfbt26V/e71eQwUw5bmpON4xaJjKAXZBzUu3I8NpU3k1+sWCl6aeEfhDAhxW/VcaAUALjQaYN3Zs5OkbgRAWwRugGy3bRaLgZXaM3CNq1sFLXl4e8vLyZnze2rVrMTAwgLfeegsXXHABAOBvf/sbBgYGcPHFFyf883p6euDxeKb85TocDjgcjoS/n97Edl6MEbw0UaOppCjMdCDDacWgL4TG7mEsK8pUe0lJ0UyjAeat2JUCG88hKIho9/oMkeBKDermhvWIam1tnTTvheM4uN1uXfaIki3nZfny5bjyyitx22234c0338Sbb76J2267DZ///OfHVBotW7YMzz33HABgaGgId911Fw4ePIimpibU1dVh06ZNyMvLwzXXXCPXUjXNaOXStP2bHBzHxR0dDam8mvkTBAGv/PVVfPK3l+FrPoJSl3FvSOTGWzgp4dkox8103ZgbI/eIkjVt+6mnnsI555yDDRs2YMOGDVi5ciV++9vfjnnO8ePHMTAwACDyi37//fexefNmLFmyBDfffDOWLFmCgwcPIiMjQ86lalaZwcqlWbIuVRrN3xKDJO3W1taisrISn73iM+j+04/R8cy9WLNiqa57UKjNSDOO+kcC0nBaSuSeval6RLndbt2WSQMy9nkBgJycHDz55JPTPid+KyslJQUvv/yynEvSnfEDtaZLdtYDuoNKnsUFkYD+eLt+gxej9qBQW4WBKo5onMj8sR5Ru2tfxLd/ewDZeQU4+stv6XLHhaGCeY1zRwdqjQSMMVCLRgMkDzs2OtGpz2MjI/egUFuZgSqOKFk3OXiex42bNyL97EsRKFiO3pGQ2kuaFwpeNM5h5VEcHail96TdkUAInYORfj0UvMwfOzY63TOsyyF8Ru5BoTYjTZc+3U3JusmSYuelXTk979gCFLzoglT6qPMtYBZ8Zafa4EqlMun5ys9wwJViQ1gEPunS3+6LkXtQqC1+vtFkO1t6Iu28UL5LUiwtiuzYHmvXd0NXCl50wCjl0k10B5VUkYojlrSrv+DFyD0o1MauGYO+EPpHgiqvZn6kMmmahZYUS6NtFWjnhcjOKAMam6SZRnQHlSx67rTLelBMlYTOcRzKysp02YNCbU4bj8LMSLm53q8bUpI/7bwkxbLozstxHV4z4lHwogPl0Z0K/R8b0c5Lsum510t8DwoYrAeFFrA5QHquOIrPk6OE3eRgx0YfdwxCCOv3SJGCFx2IDVrTd/KdVGmURxehZJEGNHbq8y6K9aDIzhs7MkTvPSi0wAjTpdmukSvFhqxUu8qrMYbK3DQ4rBb4gmFd78pR0bwOsO3SDq8fvqAAp02fd6KxHi+085IsbOeluXcEowEBKXb9vTZqampwwF+B//7zK/hslRP/cMVqVFdX047LPFUYIFeO+kIlH2/hsLgwHR+0enG83YsqneYS0c6LDmSl2pDhiMSZej068gUFtA34AFCZdDLlpTuQk2aHqNOKI+ZYxzCc5Svx91++EevWraPAJQnYzouej43oqFkeSwsjSbvHdJy0S8GLDnAcF9sC1umFiK07w2lFNpVJJ9XigsjRkR6TdgHAHxJwMtpob3mxOceAyKHcAI3qKFlXHsuK9Jvoz1DwohN6rziKn2mk9xEHWqPnpF0AONk5hFBYRKbTaogJyFrBdivavT5dNjEE6NhILrFeLxS8EJnpvdcLXYTko/cBjUfbIuteXpxJgW0SZRvguJkVKdCxUXKxnZembn125wYoeNENvXfZZZVGek0O07LFbOdFpxVHR9sinT6XF2eqvBJj0ftxcyAURmvfKAC66Um2/AwHslIj3blP6nQ2GgUvOlGu8ymxVGkkH3Zs5OkdxUhAf8PWWPByFgUvSRc/JkBvWvtHERYBp82CggyH2ssxFI7jsLRQ30dHFLzoBGs45ekdQViHjYUapZwXuoNKtpw0O/LSIz0w9HYXJYoi7bzIqDx63dDjzotUaZRDeXJykDrt6nTGEQUvOlGc5QRv4eAPhaWOk3rhDwk4M8C2f2nnRQ6LC/SZtNvh9aNvJCj1niDJFdt50V+jOsqTkxebcUQ7L0RWNt6CkiwnAP3dRXl6RyGKQJqdl3YISHLpNWmX7bosyEvTbfNFLdPzcTMFL/JaKu286OuawVDwoiMVOt0Cjm80Rdu/8lis0wGNH9GRkaxY8NLSO6q7OTbUoE5eLHjpHPSjbzig8mpmj4IXHSnL0eesEmmaNM00ko1ee71Qvou8SrJSYOM5BIQw2r0+tZczK2y3iHZe5JHusMKdHemrpMejIwpedKRCp2WPdAclP3Zs1No/imG/fiqOYsELddaVA2/h4M7WX95LOCxK1zm240yST89JuxS86Ihez6+lnRe6g5JNVqod+dFy0hM6qTjyBQWpCo3KpOXDrht66hHV7vUhEArDauGkXD+SfFLei852bAEKXnRFjxchgHZelMJ2Xz7WyRbw8fZBhEUgNy0WeJHk02OvF7ZWd3YKrDy9TcmFVRzRzguRFeuW2T0U0M3RQFAIoyXaJZOmScsrVi6tj+AlPt+FErnlo8cdW7rhUUZsQOMQRFFfCd0UvOhIptOGrOhEZr3kvbT2RaocnDYLCjPp7lpOUtKuTo6NKN9FGSwA0NN0aUrWVUZVXhpsPIchf0i6ydQLCl50pkJn06Ube2iatFL01uslfiAjkU9sqKt+EnZp50UZNt6ChfmR64be+r1Q8KIzsXJpfQQvp7vZRYjuoOTGer20Dfjg9QVVXs30aCyAcljw4vWF0D+ij34eUoO6HLpuyC2WtEvBC5GR3sqlY5VGdAclN1eKDUWZkcqMExqvHmjpG8WgPwQbz0l3fkQeKXZeGmyoh6RdURSpu66CWPCit14vFLzoTLnOjo1o+1dZi3VydMQ66y4uyIDdSpchuUkVRzq4bvQOBzDkD4HjYjvNRD567fVCVw2dKdNd8EI9XpSkl067dGSkrPK4qfRaxwKsokwnzbtSACuXPtU1jEAorPJqEkfBi86wHYyWvhHNzyoJCWF4+qLbv3m086IEKWm3U9s7L1RppCw9TZeO7dbSDY8SSlxOZDitCIVFfNKl7ZueeBS86ExRphM2nkNQENE2oO3StrYBH4KCCLvVguJM6pKpBL0MaGSVRtRZVxl6alQXS9alGx4lcByHpTq5bsSj4EVneAuHsmx9HB01Re+gynNSYbFQmbQSFhdEdl46vH4MjGqz4mjQF5Reu3RspAw95cpJwQsNclWMHpN2KXjRIb2US9NMI+VlOG0ocbGKI21eiFg/iaJMJ7LT7CqvxhxY8NLu9cEXFFRezfSkYyPaeVFMLGlXm9eMyVDwokN6KZdmPV6oTFpZizWetEv5LsrLSbMj3WGFKEby5bSMyqSVF5txRMELkZFetoCbeihZVw3SgEaN7rx8RJ11FcdxXFynXe1eN4b8IfQMRxrplVPwohiW89LaP6r5BpcMBS86pJdy6SZpNABdhJTEdl60WnFEZdLq0EPSLjsyykmzI9NpU3k15uFKjTW41MtUegpedEgPx0ZCWJRycujYSFla7vUihEVpa5qCF2WV6+C6QUdG6tFb0i4FLzrEqo36R4KarShp9/oQEMKw8RyKXVQmrSRWcdQ16NfcLJvTPcMYDQpw2iyoouNERbEEWF0EL9RZV3F6S9ql4EWH0hxW5KVHZpVotWMmS9Yty06FlaeXmZLSHFaUZqUA0N7uC+vvsrQwAzyVzytKD43qaJyIepZS8EKUUJ4TeXPS6l1UE23/qkqrSbuU76IelrDr6RtFWKPduenYSD2xYyMvRFGbr494sgYvP/zhD3HxxRcjNTUVWVlZCX2NKIrYsWMHSkpKkJKSgnXr1uHDDz+Uc5m6pPXKAbqDUhfLe9FarxcKXtRT7HLCauEQCIXR7vWpvZxJ0WgA9SwqSAdv4eD1hdDh9au9nBnJGrwEAgFce+21+NrXvpbw1/zoRz/CQw89hIcffhhvv/02ioqKcMUVV2BwUFsXYbW5s5zwNR/BS3/cg7q6OgiCthpPUaWRurTa6+UjCl5UY+UtcGdHdmy1eNPjDwloiwZVdNOjPIeVl/LQjulgwrSswcv3vvc9fOtb38I555yT0PNFUcTOnTvxne98BzU1NVixYgV+85vfYGRkBE8//bScS9WV2tpa/NvNl6PjmXvxp53fxvr161FZWYna2lq1lyY5TT1eVKXFY6P+kQDaBiJvTsuoQZ0qynNZ0q728l48vaMQRSDNziOXOi+rQk95L5rKeWlsbER7ezs2bNggPeZwOHDppZfijTfemPRr/H4/vF7vmA8jq62txZYtW9DT2Tbm8dbWVmzZskUTAYwoinE7LxS8qGFRtOKoZziAniFtbAGzXRd3dgr18FBJhYaPm+OPmjmOkrnVsKyQgpc5aW9vBwAUFhaOebywsFD63HgPPPAAXC6X9FFWVib7OtUiCAK2bt06aTIVe2zbtm2qHyF1DvrhC4bBWzhpm5ooK9VuRVmOtiqOjlJnXdVpuUcUJeuqT0+9XmYdvOzYsQMcx037cejQoXktanzULYrilJH4Pffcg4GBAenD4/HM62drWX19PVpaWqb8vCiK8Hg8qK+vV3BVEzVGy6Td2SmwUZm0apYUaKvTLkvWPYuCF9VoebQIJfmrjwUvJ7uGEBLCKq9metbZfsGdd96J66+/ftrnVFZWzmkxRUVFACI7MMXFxdLjnZ2dE3ZjGIfDAYfDMaefpzdtbW0zP2kWz5MLXYS0YUlRBl491qmZvBeqNFIf+5vU5LFRL+28qK0sOxWpdh4jAQFNPcNYVKDd3LRZBy95eXnIy8uTYy2oqqpCUVER9u7di9WrVwOIVCzt378f//7v/y7Lz9ST+IAuGc+TS5M0FoAuQmqKJe2qf2wUFMI4EV0H7byopywnBWJYQMfxI3h8dwcWVpahuroaPM+rvTTqrqsBFguHxYUZeM/Tj2Ptg5oOXmTd029ubkZDQwOam5shCAIaGhrQ0NCAoaHYxXTZsmV47rnnAESOi7Zt24b7778fzz33HD744APccsstSE1NxY033ijnUnWhuroabrd7yiM0juNQVha5GKmJdl60YXFBrNeL2k2nTnUNIyCEke6wUh6Uil76y5/Q9p9fRccz9+LWf/iKZioVhbCIlj6qUNQCvSTtznrnZTa++93v4je/+Y30b7absm/fPqxbtw4AcPz4cQwMDEjPufvuuzE6Ooo77rgDfX19uPDCC/HKK68gI0O7EaBSeJ7Hrl27sGXLFnAcN/YNKRrQ7Ny5U/W7qKZu2nnRgkUF6bBwQN9IEN1DAeRnqHe8yo6MlhVlwEJjAVTBKhXHB7KsUnHPnj2oqalRZW1n+kcRFETYeYs03ZioQy9Ju7LuvOzevRuiKE74YIELEEkyveWWW6R/cxyHHTt2oK2tDT6fD/v378eKFSvkXKau1NTUYM+ePSgtLR3zuCu3UNWLDyOKIu28aITTxksJmmp32qV8F3VpvVKRHRmV5aTQzCuV6WVAI5WC6FBNTQ2ampqwb98+bL//YRTecD8uuvcZ1QMXAOgeCmA4IIDjIJXqEvXEOu2qeyGizrrq0nql4uleuuHRCrbz0tw7gpFASOXVTI2CF53ieR7r1q3Dvd+4Fc7ylTjWMYzOQfXnlbBdlxJXChxW9ZMAzU5K2u1UN2k31uOFjn/VoPVKRerxoh256Q7kpUeOmLWQ7D8VCl50LjfdgRWlkbvZ1090q7yauEqjPLoIaYEWBjR2DfrRPeQHx8Xu6oiytFypKAgCDr5+AMMf7cdo0xHVm2yS+KMj7Xasp+DFAC5ZnA8AOPBxl8oroUojrWEVRx93DKlWccTyXapy05Bql7VGgExBq5WKtbW1qKysxF8euB3df/4x/v2bN2qi+sns9JC0S8GLAVRHg5fXT3YjHFa3JJbtvFRR8KIJC/LTYOGAgdEgugbVmXFEybrqY5WKwMQO5pxKlYqs+ml8Lo6W5rSZlR4GNFLwYgBrKrKRaufRPRTAUZW3+Zq62c4LHRtpgdPGS8Mx1Tq/jgUvdGSkpqkqFUtKShWvVNR69ZPZ6aHiiIIXA7BbLVi7IBcAcOBj9fJexkyTpkZTmrFY6rSrzoWIKo20I75S8dybv4vCG+7HT56tV7xSUevVT2a3uCADHBeZSq/Wju1MKHgxiEuWRI6O6k+ol/fSNxLEoC9SWldOLb41Q0raVWFAoy8o4JOuSEBLwYs2sErFm//+y3CWr8RrKtzwaL36yexS7Lw0pkGruy8UvBhE9eLIvKlDTX2q1OYLgoA9f3kZwx/tR2rPMdjolaUZsV4vyh8bnewcghAW4UqxodhFnVO15DPLI8NuD3zcDV9Q2eMZLVc/kYhY0q42K47oLcYgqvLS4M5OQUAI42+nehX92axi4J+u/wK6//xjHP3VXVQxoCFL4o6NlK44+igu32WqSheijrNLMlGU6cRoUMDBUz2K/mytVj+RmKVFkZ1S2nkhsuI4Tqo62q9gyTRVDGhfVV4aeAuHQV8IHV5lz6+p0ki7OI7DZ84qAAD89aMORX82q36aLJZWq/qJjMWSdtXuzj0VCl4M5NIlkaMjpfJeqGJAHxxWXhqSqfSFiAUvZ1Hwokns6OjVo52K78rV1NTg3H/4PviMvDGPu91uTcxpM7ulRbHjZrVbcEyGghcDWbswDxYO+KRrGK39o7L/PKoY0I8lKsw4EkUxbiwABS9adNGCXKTaebR7ffjwjLK5DSc7h9BXsBoVd/waf37xFTz99NPYt28fGhsbKXDRgMrcNDisFowGBTT3jqi9nAkoeDEQV4oN55ZlAQDqFTg6oooB/VBjQGPbgA8Do0FYLZxUrk20xWnjpQ7dexU+OvrLkTMAgEuWFuLzV16BG264AevWraOjIo3g4/5utdhpl4IXg2El0wcUODqiigH9iCXtKldxxI6MFuan05BODfvMWZGjo78eVS54EUURfzkSuan5/MoSxX4umZ2lhdpN2qXgxWCkUQEnuiHIfE5JFQP6wY6NTnYqN+OIOuvqw/ql+eA44MMzXpxR4LgZAI53DOJk5xDsvAVXnF2oyM8ksyd12u3QXrk0BS8Gs8rtQqbTCq8vhPda+mX9WfHzUsajigFtqcxNg43nMOQP4cyAT5GfSfku+pCb7sCa8mwAwKvHOhX5mc9Hd10uWZKPTKdNkZ9JZk/LAxopeDEYK2/B3y2KVh0p0DmzpqYGNf/yE6oY0Di71YKqPDbjSJkLEZVJ64d0dKRA3kv8kdGmVXSkrGUseGnqHla8keFMKHgxICVHBfiCAo6nnIXS2x/Hz377B6oY0DCWtHtCgeBlJBBCYw+NBdCLzyyP9Hs5+EkPhvzyduj+8IwXjd3DcFgtuHw5HRlpWUGGA1mpNoTFyJGzllDwYkBsVMBhTz+8vqCsP+ulD9ox6AvBnZOOr9/4BaoY0LAlBcqNCTjePghRBPLSHcjPcMj+88j8LMxPR2VuKgJCGK/LfNPDdl0uW1aAdIdV1p9F5ofjOCwt1ObREQUvBuTOTsWC/DQIYRFvnJS37ffv3vYAAK493w2Lhdq/axmrOFJi5yWW70LJunrAcZzUsG7vR/LlvUSOjCIl0lRlpA9a7bRLwYtBsd4NcpZMn+4ZxsFTPeA44Nrzy2T7OSQ5pGOjTvk7ZlJnXf1hRzj7jnfKVqn4XssAWvpGkWLjcdmyAll+BkkuNuOIdl6IIi6Jjgo48HGXbKWxe96JdNf99KI8lGalyPIzSPJU5qbCzlswEhBk78D8ESXr6s75ldlwpdjQOxzA4eY+WX7GX96L7LpcvrwAKXY6WtYDlrR7XGPTpSl4MagLq3Jh4zm09I2iqSf5rZ2FsCgFL1+iXRddsPIWLMiPVByd6JTvLiocFnGMghfdsfEWrF8a7bYrQ8O6cFjEC+9TYzq9YcFLh9eP/pGAyquJoeDFoNIcVqypiPRukKPq6MCJLrQN+JCVasMGajKlG7ExAfIl7Xr6RjAcEGCPC5aIPshZMn3Y04czAz6kO6xYFw2SiPalO6xwZ0d21rV0dETBi4FJowJkmHP0+0ORRN2rzy2l1u86sqSAjQmQ7yLE8l0WF6bDxtMlRk8uWZIPq4XDJ13DaOweTur3/vN7kV2XK84qhNNG1ww9kTrtUvBClMCSdg9+0oNAKJy079sz5JeGuNGRkb7Eer3It/PyEXXW1a1Mpw0XLcgFALyaxKMjYcyRETWm0xstdtql4MXAzirORG6aHcMBAe8mMQHvucOtCAoizil14awSeoPSE1YufVLGiiPqrKtvrGFdMqdMv93Ui85BPzKcVmn+GtEPVnGkpaRdCl4MzGLh8Olow7pk5b2Iooj/iR4ZfelTtOuiNxW5abBbLRgNCmjpk6fiiMqk9Y2VTB863Ze0BE3W2+WzZxfBbqW3Hb1ZGpcrp9Rg15nQq8jgpH4vSZpz9F7LAD7uGILDasEXVlHFgN7wFg4L8+XLe/H6glJQRMGLPpXlpGJZUQaEsIi64/O/6QkJYbz0QTsAOjLSqwX5scGuct30zBYFLwbHRgV8cGYAPUP+eX8/1lF344oiuFJoGqwesaOjj2Uolz4WzXcpcTnhSqXXh15J3XaTkPfyt8ZedA8FkJ1qk4bGEn2x8RZZb3rmgoIXgyvIdGJZUQZEEfjfT+Y3KmAkEMKfo02m6MhIv5bImLRL+S7GcHk072X/8a55J/uzI6MrVxRR9ZmOaS1pl15JJpCskukX32/HkD+E8pxUXFSVm4ylERUslrFcmoIXY1jlzkJeugND/hDeauyd8/cJCmG8KB0Z0TGzni3VWLk0BS8mwPJe6k/Mb1TA76KJuteuoSGMesZ2Xk52DiV9hg0FL8ZgsXBS1dFf53F09L8nu9E/EkRumh0XVuUka3lEBVrr9ULBiwmcX5kNh9WCDq9/zp1VG7uH8VZjLywcsOV8d5JXSJRUlpMKh9UCfyiM5t7kjY4QwiKOd9A0aaO4XJoy3THnm56/HIn0dtl4ThGsdGSka6xc+pOuoaT2DZsrejWZgNPG48Jo46m5lkyzjrqXLMlHsYuGMOoZb+GwSIajo8buYfiCYaTYeFTk0lgAvfv0ojw4rBa09o9KQelsBEJhvPwhHRkZRYnLiQynFaGwiFPd8jW5TBQFLyZxSbTqaP8c8l5CQlgawngdddQ1hFjSbvKCF3ZktLQoAzwdK+peip2XqhXnMuuo/kQXBn0hFGQ48KlKOjLSO47jpH4vWjg6ouDFJFjS7luNvfAFhVl97f6Pu9A56EdOml3aSib6tigvFb7mI/jLc3tQV1cHQZjda2IylO9iPLGS6c5Zfy07MrrqnGIKZg1CSxVHFLyYxOKCdBRlOuEPhWddPcA66l6zupS6YxpAbW0tdnx5HTqeuRcvP3wv1q9fj8rKStTW1s7r+34kddalfBejuGxZJGn3PU8/Or2+hL/OFxSk8QKbVlFjOqPQUtIuvROZBMdx0hbwbPJeugb9eDV610VDGPWvtrYWW7ZsQXdH25jHW1tbsWXLlnkFMLTzYjwFmU6sKssCALx2LPHdl7rjXRjyh1DscmJ1WbZMqyNKi804Mnjw8sMf/hAXX3wxUlNTkZWVldDX3HLLLeA4bszHRRddJOcyTSPW7yXxUQHPHW5BKCxiVVmWtGVI9EkQBGzdunXSyhH22LZt2+Z0hNQ7HECHN9LBeRkFL4ZyxRxKplljus+dU0xtFQyE5by09o/C6wuquhZZg5dAIIBrr70WX/va12b1dVdeeSXa2tqkjxdeeEGmFZrLpxflgeOA4x2D6EhgCzgyhJESdY2ivr4eLS0tU35eFEV4PB7U19fP+nuzXZfynFSkO6xzXiPRHpbnVn+iG6OBmQPb0YAg7dZ+nuafGYor1YbCdBt8zUfwH4/tTlq+3FzIepX53ve+BwDYvXv3rL7O4XCgqKhIhhWZW3aaHStLXXivZQD1J7qxZc30/Vrebe7Hyc4hOG0WOrc2gLa2tpmfBOC0p3XW3zt2ZES7c0azrCgDpVkpaO0fxf+e7MZnzpo+af+1Y50YDQooy0nBKrdLoVUSJdTW1uKDh76Gkb5O/H/PRB5zu93YtWsXampqFF2LJnNe6urqUFBQgCVLluC2225DZ+fUZ61+vx9er3fMB5la9eLERwX8T3QI41XnFCPDSUP29K64OLEAdMdfW3FP7ft4t7kv4eZkH1G+i2FxHIcrogFLIkdHsSOjEnAcHRkZBcuXG+kb+36cjHy5udBc8LJx40Y89dRTeO211/CTn/wEb7/9Ni677DL4/ZNPRH7ggQfgcrmkj7IyOt6YDst7ef1kN8LTtIYf9oekixAdGRlDdXU13G73NG8oHOyufIQLl+GZt5pR88gbuOKnB/Cf+z9B5+Dkx4yCIKCurg6v/uU5+JqPYEkBNaczosulvJfOaa8bQ/6QlNj7+ZW0W2sUcubLzdWsg5cdO3ZMSKgd/3Ho0KE5L+i6667D5z73OaxYsQKbNm3Ciy++iI8//hjPP//8pM+/5557MDAwIH14PJ45/2wzWF2ehTQ7j97hAD48M/Uu1fPvt2E4IKAqLw0X0EwSQ+B5Hrt27QKACQFM5G8XeOpXv8B//9PfoWZ1KZw2C052DuGBF49h7QOv4au738ZLH7RJrcFra2tRWVmJ9evX48hvv4+OZ+7FrRsvVPwOjMjvwqpcpDus6B7y40jrwJTPe/VoB/yhMCpzU3F2Ce3CGYWc+XJzNeuclzvvvBPXX3/9tM+prKyc63omKC4uRkVFBU6cODHp5x0OBxwOR9J+ntHZeAvWLszDX4924MCJLpwzxZk0OzK69vzp7tSJ3tTU1GDPnj3YunXrmIuR2+3Gzp07pXPrtQtz8b3NZ+P5I234/TsteOd0H1491olXj3UiJ82OxSMf4vf/vn3CnVh72xls2bIFe/bsUfwMnMjHbrXg0qX5eP5IG/76UQfOjZZPj/fn9yJ5VZ9fSUdGRpJovlyiz0uGWQcveXl5yMvLk2Mtk+rp6YHH40n4vJ7M7NIl0eDl4y58ff2iCZ//pGsIh073wcIBXzyPhjAaTU1NDTZv3oz6+nq0tbWhuLgY1dXV4Hl+zPMynDZcf0E5rr+gHCc7h7DnnRbUvtuCjoERHHn0h1NuIXMch23btmHz5s0TvifRr88sL4gEL0c7cNdnl074vNcXlHLpPk8J/oaS6Puvku/Tsua8NDc3o6GhAc3NzRAEAQ0NDWhoaMDQUGyo07Jly/Dcc88BAIaGhnDXXXfh4MGDaGpqQl1dHTZt2oS8vDxcc801ci7VVFjS7rvNfRjyhyZ8nnXUXb+0AIWZTkXXRpTB8zzWrVuHG264AevWrZsxyFhUkI5vb1yGN759GbatECAMTt0rSI0tZCK/9UsLwFs4HGsfhGeSaeR7P+xAQAhjUUG61A+EGMNM+XIcx6GsrAzV1dWKrUnW4OW73/0uVq9ejfvuuw9DQ0NYvXo1Vq9ePSYn5vjx4xgYiJyh8jyP999/H5s3b8aSJUtw8803Y8mSJTh48CAyMuiPIVkq89JQnpOKoCDizU96xnwuKITx7DuRUtkvfYoSdclYVt6CIltibeKV3EIm8stKteP8iki33FcnqTpiCf6fX1lMR0YGM1O+HADs3LlT0Z1WWYOX3bt3QxTFCR/r1q2TniOKIm655RYAQEpKCl5++WV0dnYiEAjg9OnT2L17N1UQyWCqUQF1x7vQPeRHXrpdmmtCSDwtbiETZcRKpseWy/aPBFB/IrIbR1VGxsTy5UpLS8c87na7Vclx01ypNFGGNCrgxNjt/99FE3VrznPDxtPLg0ykxS1kogzWbffNUz1j2sO//GE7QmERy4oysKiAdsmNqqamBk1NTdi3bx+efvpp7Nu3D42Njaok59O7k0mtXZgL3sKhsXtYOr/u9Pqw7zgbwkiJumRyWtxCJsqoykvDwvw0hMLimEaXfznCqoxo18XoZpsvJxcKXkwq02nDeeVZAIAD0aOj2sOtEMIizivPorsnMi2tbSET5bDxAH/9KJL30jPkxxvR3LnPr6RZRkQZFLyYGKs6qv+4OzKEMXpkdB0l6pIEaGkLmSjnM9Gjo33HuxASwnjpw3YIYRErSjNRmUcdlokyaPyriV2yJB8/efkoXn71Vfxrzzv46N0eZC9Yic/R3RNJENtCJuZxXnk2slNt6B3y4bHf/RlPvvYefCM2XLXhi2ovjZgIBS8mdvJvr6LtP29H0NuN+34TeWw4txCvnPcI3T0TQibFWziUDXyA93/5AL4e1+/nB/t/jqL/+BldO4gi6NjIpGpra/GlL12LoHdstdFQb6cqE0IJIfpQW1uLP//0rgmNCtloCLp2ECVwYqIz73XC6/XC5XJhYGAAmZk0GGwygiCgsrJyykFbHMfB7XajsbGRKkYIIRK6dhA5zeb9m3ZeTEiLE0IJIdpH1w6iFRS8mJAWJ4QSQrSPrh1EKyh4MSFq704ImQu6dhCtoODFhKi9OyFkLujaQbSCghcTovbuhJC5oGsH0QoKXkyK2rsTQuaCrh1EC6hU2uQEQUB9fT3a2tpQXFyM6upqumsihMyIrh0k2Wbz/k3BCyGEEEJUR31eCCGEEGJYFLwQQgghRFcoeCGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdIWCF0IIIYToCgUvhBBCCNEVCl4IIYQQoitWtReQbKxhsNfrVXklhBBCCEkUe99OpPG/4YKXwcFBAEBZWZnKKyGEEELIbA0ODsLlck37HMPNNgqHwzhz5gwyMjImjGyfL6/Xi7KyMng8HpqbJCP6PSuDfs/Kod+1Muj3rAy5fs+iKGJwcBAlJSWwWKbPajHczovFYoHb7Zb1Z2RmZtIfhgLo96wM+j0rh37XyqDfszLk+D3PtOPCUMIuIYQQQnSFghdCCCGE6AoFL7PgcDhw3333weFwqL0UQ6PfszLo96wc+l0rg37PytDC79lwCbuEEEIIMTbaeSGEEEKIrlDwQgghhBBdoeCFEEIIIbpCwQshhBBCdIWClwQ98sgjqKqqgtPpxJo1a1BfX6/2kgxnx44d4DhuzEdRUZHay9K9AwcOYNOmTSgpKQHHcfjDH/4w5vOiKGLHjh0oKSlBSkoK1q1bhw8//FCdxerYTL/nW265ZcLr+6KLLlJnsTr2wAMP4FOf+hQyMjJQUFCAq6++GsePHx/zHHpNz18iv2c1X9MUvCTgd7/7HbZt24bvfOc7OHz4MKqrq7Fx40Y0NzervTTDOfvss9HW1iZ9vP/++2ovSfeGh4exatUqPPzww5N+/kc/+hEeeughPPzww3j77bdRVFSEK664QpoTRhIz0+8ZAK688soxr+8XXnhBwRUaw/79+/H1r38db775Jvbu3YtQKIQNGzZgeHhYeg69pucvkd8zoOJrWiQzuuCCC8Tbb799zGPLli0Tv/3tb6u0ImO67777xFWrVqm9DEMDID733HPSv8PhsFhUVCT+27/9m/SYz+cTXS6X+Oijj6qwQmMY/3sWRVG8+eabxc2bN6uyHiPr7OwUAYj79+8XRZFe03IZ/3sWRXVf07TzMoNAIIB33nkHGzZsGPP4hg0b8MYbb6i0KuM6ceIESkpKUFVVheuvvx6nTp1Se0mG1tjYiPb29jGvb4fDgUsvvZRe3zKoq6tDQUEBlixZgttuuw2dnZ1qL0n3BgYGAAA5OTkA6DUtl/G/Z0at1zQFLzPo7u6GIAgoLCwc83hhYSHa29tVWpUxXXjhhXjiiSfw8ssv45e//CXa29tx8cUXo6enR+2lGRZ7DdPrW34bN27EU089hddeew0/+clP8Pbbb+Oyyy6D3+9Xe2m6JYoitm/fjk9/+tNYsWIFAHpNy2Gy3zOg7mvacFOl5cJx3Jh/i6I44TEyPxs3bpT+9znnnIO1a9di4cKF+M1vfoPt27eruDLjo9e3/K677jrpf69YsQLnn38+Kioq8Pzzz6OmpkbFlenXnXfeiSNHjuD111+f8Dl6TSfPVL9nNV/TtPMyg7y8PPA8PyFi7+zsnBDZk+RKS0vDOeecgxMnTqi9FMNi1Vz0+lZecXExKioq6PU9R9/4xjfwpz/9Cfv27YPb7ZYep9d0ck31e56Mkq9pCl5mYLfbsWbNGuzdu3fM43v37sXFF1+s0qrMwe/34+jRoyguLlZ7KYZVVVWFoqKiMa/vQCCA/fv30+tbZj09PfB4PPT6niVRFHHnnXeitrYWr732GqqqqsZ8nl7TyTHT73kySr6m6dgoAdu3b8dNN92E888/H2vXrsVjjz2G5uZm3H777WovzVDuuusubNq0CeXl5ejs7MQPfvADeL1e3HzzzWovTdeGhoZw8uRJ6d+NjY1oaGhATk4OysvLsW3bNtx///1YvHgxFi9ejPvvvx+pqam48cYbVVy1/kz3e87JycGOHTvwxS9+EcXFxWhqasK9996LvLw8XHPNNSquWn++/vWv4+mnn8Yf//hHZGRkSDssLpcLKSkp4DiOXtNJMNPveWhoSN3XtCo1Tjr085//XKyoqBDtdrt43nnnjSkXI8lx3XXXicXFxaLNZhNLSkrEmpoa8cMPP1R7Wbq3b98+EcCEj5tvvlkUxUhp6X333ScWFRWJDodDvOSSS8T3339f3UXr0HS/55GREXHDhg1ifn6+aLPZxPLycvHmm28Wm5ub1V627kz2OwYg/td//Zf0HHpNz99Mv2e1X9NcdJGEEEIIIbpAOS+EEEII0RUKXgghhBCiKxS8EEIIIURXKHghhBBCiK5Q8EIIIYQQXaHghRBCCCG6QsELIYQQQnSFghdCCCGE6AoFL4QQQgjRFQpeCCGEEKIrFLwQQgghRFcoeCGEEEKIrvz/JWsDNP8xzmoAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from pylab import plot,show\n", "\n", "# Constants\n", "N = 26\n", "C = 1.0\n", "m = 1.0\n", "k = 6.0\n", "omega = 2.0\n", "alpha = 2*k-m*omega*omega\n", "\n", "# Set up the initial values of the arrays\n", "A = np.empty([N,N],float)\n", "for i in range(N):\n", " if i>0:\n", " A[i,i-1] = -k\n", " A[i,i] = alpha\n", " if i