"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.linalg as sl\n",
"import scipy.optimize as sop\n",
"from pprint import pprint"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Checking properties of transposes\n",
"\n",
"Check the four properties of the transpose operator from the lecture for some example matrices.\n",
"\n",
"Hint: \n",
"\n",
"```Python\n",
"A = np.random.random((m,n))\n",
"```\n",
"\n",
"is a convenient way to generate an arbitrary m by n matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- I made a deliberate mistake here - what was it? Can you fix this test?\n",
"\n",
"- Why did I make the above choices in the use of `allcose` vs `array_equal`?"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(A+B)^T = A^T + B^T: True\n",
"(alpha * A)^T = alpha * A^T: True\n",
"(A^T)^T = A: True\n"
]
},
{
"ename": "ValueError",
"evalue": "matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 10 is different from 20)",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'(alpha * A)^T = alpha * A^T: '\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0malpha\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'(A^T)^T = A: '\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray_equal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'(AB)^T = B^T A^T: '\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m@\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mB\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m \u001b[1;33m@\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;31mValueError\u001b[0m: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 10 is different from 20)"
]
}
],
"source": [
"m = 10\n",
"n = 20\n",
"A = np.random.random((m,n))\n",
"B = np.random.random((m,n))\n",
"print('(A+B)^T = A^T + B^T: ',np.allclose((A+B).T, A.T + B.T))\n",
"alpha = np.pi\n",
"print('(alpha * A)^T = alpha * A^T: ',np.allclose((alpha * A).T, alpha * A.T))\n",
"print('(A^T)^T = A: ',np.array_equal((A.T).T, A))\n",
"print('(AB)^T = B^T A^T: ',np.allclose((A@B).T, B.T @ A.T))\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(A+B)^T = A^T + B^T: True\n",
"(alpha * A)^T = alpha * A^T: True\n",
"(A^T)^T = A: True\n",
"(AB)^T = B^T A^T: True\n"
]
}
],
"source": [
"m = 20\n",
"n = 20\n",
"A = np.random.random((m,n))\n",
"B = np.random.random((m,n))\n",
"print('(A+B)^T = A^T + B^T: ',np.allclose((A+B).T, A.T + B.T))\n",
"alpha = np.pi\n",
"print('(alpha * A)^T = alpha * A^T: ',np.allclose((alpha * A).T, alpha * A.T))\n",
"print('(A^T)^T = A: ',np.array_equal((A.T).T, A))\n",
"print('(AB)^T = B^T A^T: ',np.allclose((A@B).T, B.T @ A.T))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Prove that $(AB)^T = B^T A^T$\n",
"\n",
"As per the section title!\n",
"\n",
"Hint: Writing in index notation $A = [a_{ij}]_{m\\times n}$ and $B = [b_{ij}]_{n\\times p}$\n",
"the product can be defined as: the $(i,j)$-th entry of the product $AB$ is given by\n",
"$$\\sum_{k=1}^n a_{ik} b_{kj}$$\n",
"with $A$ being an $m\\times p$ matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution\n",
"\n",
"Something like this:\n",
"\n",
"Suppose $A = [a_{ij}]_{m\\times n}$ and $B = [b_{ij}]_{n\\times p}$.\n",
"\n",
"Define $C=(AB)^T$ and $D= B^T A^T$. We need to show that $c_{ij}=d_{ij}$ $\\forall i, j$.\n",
"\n",
"The $(i,j)$-th entry of the product $AB$ is given by \n",
"\n",
"$$\\sum_{k=1}^n a_{ik} b_{kj}$$\n",
"\n",
"$c_{ij}$, the $(i,j)$-th entry of the transpose of this product ($C:=(AB)^T)$, \n",
"is given by the above but with the $i$ and $j$ indices swapped:\n",
"\n",
"$$c_{ij} = \\sum_{k=1}^n a_{jk} b_{ki}$$\n",
"\n",
"$d_{ij}$, the $(i,j)$-th entry of $B^T A^T$, is given by the product of $B^T$ with $A^T$ - we take rows of $B^T$ and multiply them against columns of $A^T$ (that's the same as columns of $B$ against rows of $A$):\n",
"\n",
"$$d_{ij} = \\sum_{k=1}^n b_{ki} a_{jk}$$\n",
"\n",
"Thus $c_{ij} = d_{ij}$ $\\forall i,j$, and so $C=D$, i.e. $(AB)^T = B^T A^T$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Symmetric matrices\n",
"\n",
"1. Given an $n\\times n$ symmetric matrix $A$ and an arbitrary $m\\times n$ matrix B, show that the matrix $BAB^T$ is symmetric.\n",
"\n",
"Verify through an example using NumPy.\n",
"\n",
"\n",
"\n",
"2. If $A$ and $B$ are both symmetric square matrices, are $AB+BA$ and $AB-BA$ symmetric or skew-symmetric?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution\n",
"\n",
"Using the properties of transposes\n",
"\n",
"$$(BAB^T)^T = ((BA)B^T)^T = (B^T)^T (BA)^T = B A^T B^T = B A B^T $$\n",
"\n",
"and so $BAB^T$ is symmetric.\n",
"\n",
"\n",
"Similarly $AB+BA$ is symmetric but $AB-BA$ is skew symmetric."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 21 51 81]\n",
" [ 51 121 191]\n",
" [ 81 191 301]]\n",
"Is symmetric: True\n"
]
}
],
"source": [
"# symmetric 2x2 matrix\n",
"A = np.array( [[1, 4], [4,1] ] )\n",
"# arbitrary 3x2 matrix\n",
"B = np.array( [[1, 2], [3,4], [5,6] ] )\n",
"\n",
"print((B@A)@B.T)\n",
"\n",
"print('Is symmetric: ',np.array_equal((B@A)@B.T, ((B@A)@B.T).T))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Matrix-vector multiplication as a weighted sum of columns\n",
"\n",
"In the lecture we pointed out \"another useful interpretation\" of matrix vector multiplication.\n",
"\n",
"Code up this approach and perform the same testing that we did in the lecture for `def mat_vec_product(A, x)`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 5. -3.]\n",
"mat_vec_product2(A, x) == A@x: True\n",
"[ 5. -3. 2.]\n",
"mat_vec_product2(A, x) == A@x: True\n",
"mat_vec_product2(A, x) == A@x: True\n"
]
}
],
"source": [
"def mat_vec_product2(A, x):\n",
" \"\"\"Function to compute the matrix vector product Ax - a different approach.\n",
"\n",
" Parameters\n",
" ----------\n",
" A : ndarray\n",
" Some matrix A with shape (m,n)\n",
" x : array_like\n",
" Some vector x with shape (n,)\n",
"\n",
" Returns\n",
" -------\n",
" b : ndarray\n",
" RHS vector b with shape (m,)\n",
" \"\"\"\n",
" m, n = np.shape(A)\n",
" assert x.ndim == 1 # restrict to the case where x is 1D\n",
" assert n == len(x) # as 1D we can check the length of x is consistent with A\n",
" b = np.zeros(m) # and then initialise the appropriate length array for b\n",
" for j in range(n):\n",
" b[:] += x[j] * A[:, j]\n",
" return b\n",
"\n",
"\n",
"# Let's test against the Numpy product for a simple 2x2\n",
"A = np.array([[2, 3], [1, -4]])\n",
"x = np.array([1, 1])\n",
"\n",
"print(mat_vec_product2(A,x))\n",
"print('mat_vec_product2(A, x) == A@x: ', np.allclose(mat_vec_product2(A, x), A@x))\n",
"\n",
"# Check a non-square matrix\n",
"A = np.array([[2, 3], [1, -4], [1,1]])\n",
"# based on the above weighted sum of columns, what should Ax be for the following x?\n",
"x = np.array([1, 1])\n",
"\n",
"print(mat_vec_product2(A,x))\n",
"print('mat_vec_product2(A, x) == A@x: ', np.allclose(mat_vec_product2(A, x), A@x))\n",
"\n",
"\n",
"# and we can use random matrices to test a larger case\n",
"A = np.random.random((100,100))\n",
"x = np.random.random(100)\n",
"print('mat_vec_product2(A, x) == A@x: ', np.allclose(mat_vec_product2(A, x), A@x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Fitting two data points\n",
"\n",
"We stated in class that:\n",
"\n",
"\n",
"The polynomial that fits the two data points $\\{(x_0,y_0),(x_1,y_1)\\}$ is clearly the linear function given by\n",
"\n",
"$$ y = f(x) \\equiv a_0 + a_1\\,x \\;\\;\\;\\;\\; \\text{i.e. the degree one polynomial:} \\;\\;\\;\\;\\; y = P_1(x) \\equiv a_0 + a_1\\,x$$\n",
"\n",
"where through substitution we arrive at two simultaneous equations (or a $2\\times 2$ matrix system) which can fairly easily be solved by substituting one equation into the other to conclude that\n",
"\n",
"$$ a_0 = y_0 - \\frac{y_1-y_0}{x_1-x_0}x_0, \\;\\;\\;\\;\\;\\;\\;\\; a_1 = \\frac{y_1-y_0}{x_1-x_0}. $$\n",
"\n",
"Form the set of two simultaneous equations and solve by hand to derive this solution for the coefficients.\n",
"\n",
" \n",
"\n",
"Confirm that this result is exactly the same as the Lagrange polynomial (we saw in the lecture) that you can just write down without needing to invert any system for coefficients."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution\n",
"\n",
"We have two pieces of information - let's write these out:\n",
"\n",
"\\begin{align*}\n",
"(1) & \\;\\;\\;\\; y_0 = a_0 + a_1\\,x_0, \\\\[5pt]\n",
"(2) & \\;\\;\\;\\; y_1 = a_0 + a_1\\,x_1. \n",
"\\end{align*}\n",
"\n",
"We are assuming we know the $x$'s and the $y$'s, and we want to find the $a$'s.\n",
"\n",
"There are multiple ways we could solve this, one way is to rearrange the second equation to give an expression for $a_1$ in terms of $a_0$:\n",
"\n",
"$$ a_0 = y_1 - a_1x_1,$$\n",
"\n",
"and substitute this into the first equation:\n",
"\n",
"$$ y_0 = a_0 + a_1\\,x_0 = y_1 - a_1x_1 + a_1\\,x_0 = y_1 - (x_1-x_0)a_1,$$\n",
"\n",
"which we can solve for $a_1$:\n",
"\n",
"$$ a_1 = \\frac{y_1 - y_0}{x_1 - x_0}.$$\n",
"\n",
"We can now substitute this into either one of the two original equations to find $a_0$ - consider (1) for example:\n",
"\n",
"$$y_0 = a_0 + a_1\\,x_0 = a_0 + \\frac{y_1 - y_0}{x_1 - x_0} x_0.$$\n",
"\n",
"Rearranging we get the solution we want:\n",
"\n",
"$$a_0 = y_0 - \\frac{y_1 - y_0}{x_1 - x_0} x_0.$$\n",
"\n",
"Note that this is equivalent to forming and solving the linear system\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 & x_0 \\\\\n",
"1 & x_1 \n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"a_0\\\\\n",
"a_1\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"y_0\\\\\n",
"y_1\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We therefore have\n",
"\n",
"$$ y = f(x) = P_1(x) \\equiv a_0 + a_1\\,x$$\n",
"\n",
"where \n",
"\n",
"$$ a_0 = y_0 - \\frac{y_1-y_0}{x_1-x_0}x_0, \\;\\;\\;\\;\\;\\;\\;\\; a_1 = \\frac{y_1-y_0}{x_1-x_0}. $$\n",
"\n",
"Collecting terms we have\n",
"\n",
"\\begin{align}\n",
"y = a_0 + a_1\\,x \n",
"& = y_0 - \\frac{y_1-y_0}{x_1-x_0}x_0 + \\frac{y_1-y_0}{x_1-x_0} x\\\\\n",
"& = y_0 - \\frac{y_1\\, x_0}{x_1-x_0} + \\frac{y_0\\, x_0}{x_1-x_0} + \\frac{y_1\\, x}{x_1-x_0} - \\frac{y_0\\, x}{x_1-x_0} \\\\\n",
"& = \\frac{1}{x_1-x_0}\n",
"\\left[ \n",
"y_0(x_1-x_0) - y_1\\, x_0 + y_0\\, x_0 + y_1\\, x- y_0\\, x\n",
"\\right] \\\\\n",
"& = \\frac{1}{x_1-x_0}\n",
"\\left[ \n",
"y_0 (x_1 - x) + y_1 (x - x_0)\n",
"\\right] \\\\\n",
"& = \\frac{x - x_1}{x_0-x_1}\\,y_0 + \\frac{x - x_0}{x_1-x_0}\\,y_1\n",
"\\end{align}\n",
"\n",
"\n",
"and this does indeed agree with the expression we can just write down based on our knowldge of the Lagrange polynomial!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Squared error calculation (from Numerical Methods)\n",
"\n",
"As described in the docs ([numpy.polyfit](http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html)), least squares fitting minimises the square of the difference between the data provided and the polynomial,\n",
"\n",
"$$E = \\sum_{i=0}^{N} (p(x_i) - y_i)^2,$$\n",
"\n",
"where $p(x_i)$ is the value of the polynomial function that has been fit to the data evaluated at point $x_i$, and $y_i$ is the $i^{th}$ data value.\n",
"\n",
"Write a Python function that evaluates the squared error, $E$, and use this function to evaluate the error for each of the polynomials calculated below. \n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# create the polynomials using code from the lecture\n",
"\n",
"# consider the above example data again\n",
"xi=np.array([0.5,2.0,4.0,5.0,7.0,9.0])\n",
"yi=np.array([0.5,0.4,0.3,0.1,0.9,0.8])\n",
"\n",
"# Calculate coefficients of polynomial degree 0 - ie a constant value.\n",
"poly_coeffs=np.polyfit(xi, yi, 0)\n",
"\n",
"# Construct a polynomial function which we can use to evaluate for arbitrary x values.\n",
"p0 = np.poly1d(poly_coeffs)\n",
"\n",
"# Fit a polynomial degree 1 - ie a straight line.\n",
"poly_coeffs=np.polyfit(xi, yi, 1)\n",
"p1 = np.poly1d(poly_coeffs)\n",
"\n",
"# Quadratic\n",
"poly_coeffs=np.polyfit(xi, yi, 2)\n",
"p2 = np.poly1d(poly_coeffs)\n",
"\n",
"# Cubic\n",
"poly_coeffs=np.polyfit(xi, yi, 3)\n",
"p3 = np.poly1d(poly_coeffs)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGKCAYAAAARyQg4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAByeElEQVR4nO3dd3gUVRfA4d9NoYTQe03ogkoXQVRAFAhFeg1dCQqoiPqJoKJIEBtiQTEgIBBAem8CBrAg0hSk9xJqaAkhgST3++MuIQkBQrK7k2zO+zx5YGdnZ87cnd2zc+cWpbVGCCGEEK7DzeoAhBBCCGFfktyFEEIIFyPJXQghhHAxktyFEEIIFyPJXQghhHAxktyFEEIIFyPJ3QmUUr2UUjrBX7hS6h+l1ECllMcDbitEKRXioFCdQinlayuHXql4bYY//rRSSk1RSh21YL+9lFJ97rJcK6V8nR2TPSil3JRSY5VSp5VScUqphQ/4+ga242+QYNkd56lSqq5S6i+l1DXb+tVsy4cqpY4rpWKUUjvSejy2baar90Qp9YFSSvpdO9EDJRaRZh2Ak0Au2/+/AQoB71sZlAVOA3WBQ1YHkkF9BHxlwX57Yb4zJiVZvgzzfp52dkB20h54DXgD+BMIs8M2+yez7EfgOtASiAT2K6VqA4HAZ8BCINwO+4aM/56INJLk7lw7tNYHbf9frZQqBwwikyV3rXU0sMnqOJxFKZXVdsx2obVOVz+KtNbngfNWx5EGlWz/jtVax9ljg1rr3QkfK6XcgIpAoNZ6XYLlt/Y9Xmt9OK37VUp5AjEu8J6kmVJKAZ5a6xtWx2IFqZa31t9ATqVUIQClVFOl1J9KqetKqStKqYVKqYp3e7FSqohS6oZS6rVknvtAKRWplMprexyilPpNKfWsUmqb7bldSqnWybz2vnEk2F5TpdQO27rblVKPK6U8lFKjbNWcF23VyDkSvPaOanml1GNKqblKqZO2be2zbSN7KsoVpdSHtuO8opS6oJRap5Sqk2SdW9Wp7WwxXlJKXVVKBSul8idZVyulApVSwxLEuOFW1Woy5dLSVh7R2K7ilFK1lVJrlFIRtqrZtbYrt1uvHWl7Px9LsCyHrSz+VLZbOCpJtXyC8nxJKfWxUuqMMrd+piulvJRS5ZRSq2z7PaiU6pkk5nJKqWlKqSO24zqslPr+1rlz67iA+kA9dfv2UojtuTuqgJVSnrbjOWo7pqO2x57JxN1PKTXCdr5cVkotUUqVSBJjV1t5Rtje051KqX73Pgvufy7byvED28PYpOdlMtsrqJSaYTtPLiulpgJ5klkvvlretr1YzPfte7Z9HLU9P8X2kkO25R/YXuOhlHpHKbVXKRWtlApVSn2hlMqWTPn1V0p9qpQKBaKBPHd5T47azonOSqk9tnNwi1LqyWTif822fpRSarNS6gnb4ylJ103mtdWVUhttrz2llHoPUMmsd99jtK1XRim1XJnvrHO2dQLucXx9lFJ7gRtAc9tzVZVSi5X5jF9XSv2ulHoqmZjqK/O5DLeVzyql1CP3O+Z0SWstfw7+w1RnaqBckuVzgBjAC2iK+QL4BXge6AocxPz6Lp7gNSFASILHs4HdSbbrDpwApiR53WngP6CbbX+/2PZfLsF6DxLHGWAn0BloAewGzgITgMlAE0zNxE3g0wSv9bWVR68Ey9oB79q2Ux+TEM8As5IcW6Ljv0eZTwS6Aw1t25yF+bBXSbBOA1scJ2zxNgVewVSN/ppke7fW+x1oDXQC9mGqcPMlie8ccAToY9tHFdvfdWArphq4HebH3XWgqu21HsAfwAHA27ZsCnAFKJ1gH1OAo8mU5zHgJ1u5v24r96m29+hV4DlgPhAHPJzg9U8DHwOtbP/vBewH/kywTmVgG/APUMf2VznJ+e2bYP0ZmHNrBNAYGG6LZ0YycR+1re8H9AQuAOsTrPekLeaxwLO27b0KvH2fc+C+5zJQ3fbe6wTHVfAe29wIXAUG2sp5ku280ECD5M5ToCBQz7bORNs+qtvKdJRteRvb8hK218wCrmFq9Z7FnJeXgXnJlN8pTJV+C9t7mP0u78lRzDnyN+YcbAFst203T4L1XkwQaxPMZ/Gobb0p9ynzAsAlYA/mM9Ia85k5Aegk66bkGLNgbt+dsh1TM8w5fOwux3cK2AV0ARoBZYEatv38ZjvuZsBizA+hmgle3xxzzi6ylWMrzOfxElDSqvyR2j/LA8gMfwk+aBUxX+B5gX6YL56FtnW2YL7UPRK8rjTmC3FMgmUhJE7uDWzbfirBsudty+oked1NoHyCZYVsMQxNsOxB4rgJlElmv2uSHP984EiCx74kSe5J1le2cuqG+VLPf7fjT2H5u9u2tw/4KpmyW5lkfX/b8kYJlmlM0smR5DhuAh8liS8OqJZkm3O580s0F3ARmJ9km5cxSbqLbb9dk2xrCskn93XJlLsGuiVYlhfzBTb8HuXlgUmoGqie5Nh+u8f57Wt7/Ijt8QdJ1nvXtrxKkrjXJ1nvTdvyYgkeX0zF5y6l5/JIkiSeu2zvOVtcnZMsX8E9knuCMk2uTG4lUt8Ey56yLetxl/OyWpLy2waoe70ntmVHMYkqb4JltRKeY5jahRPA8iTba2tbb8p9yigQ8yO6VIJlOTCfHZ2KYwywPa6dYB2F+ZGZ3PFFAkWSbHMt5sdGliTfCXuwff/alh0E1iZ5bS5b7GMf9Pyz+k+q5Z1rL+aL5SLwHRAM9FGmyroG8LPWOubWylrrI5hfvfXvtkGtdQjmijlhFWU/4F+tddL72ge01gcSvPYc5iqzFJgq4AeMY79OfJ9wr+3fVckcdwml1B1Vc7copXIppT5RSh3C/KK+CUzDfJDL3+1199jes0qpX5VSYZhkdhOogPmBldTsJI/nYBJ03STLl2utr916oLU+imk7kHS9o1rrHUmWPQ0s1VpfTvD6q5griPoJlh0FXgJ6YK4op2qtZ9ztOJNYkeTxHe+H1voS5j0veWuZUiqLMi229yqlrmPKaqPt6bveFrqHp23/Tk+y/NbjpOfRsiSPd9r+LWX7928gr63KtYVSKs/9AkjLZ+oe6mJ+DM9LsnxWKrZ1L00xCXKereraQ5lbMqttzz+dZP2F2paJUuBP2zlwS9KyLmH7m5PkdYswn6P7qQts0lofv7XA9plZkmS9lB5jHeC41npzgu1p7nwPbtmktT5z64Eyt/Xq244nLsF+FLDm1n6UUuUxV/nBSeKJxDSyTFrm6Z4kd+dqAzwGPIS5Auyhtb6IuZpSJN+y9QyQ7z7b/R5or5TKr5TywXxwxiez3sVklkUDt+5xPWgcl5I8vnGP5R6YX8t3MxmT1L7GXCE9BgywPZftbi9KjlKqBrAciABewHxBPIb5tZ/cts4mfKBNA5xLQPF7rZdgWdL1kiu/fHdZfgZT7gktw1T3ZwW+TOY1d/Mg70fCcvgYc995OqZqsjbmSg0esOxtbp0nSY/3TJLnb0l6Xt5qfJgNQGu9HtO7pCSwADivTNuFKveIIa2fqeQUBS5prW8mWZ7ceZEWhTDV0RGYH1q3/s7Zns+fZP0HaRGfqKz17Yaet97norZ/zyVZLxZzBXs/Rbn75yShlB5j0aSx3GV7tyQti3yY7533kuznJubWSl5lGjsWsq3/YzLrteDOMk/3pLW8c+3St1vLJ3QJU8VUJJnninD/rjlTMV/QvTBfatcxtQIPKq1xpIqtAU0rTJXlVwmWP5rKTbbDXGW0TfhFrEwDscvJrF84STxZMOV46l7rJViWdL3krqIucvdyTZrcxmG+kA4BQUqpeskkFHvqjKkhGHlrgVLKOw3bu3U8RUjc3fHW8T/weaS1ngvMtcXVAPgEWKmUKqGTb+HuiHP5NCYZeCZ5P5I7L9IiDIjCVF0nJzTJ45RetafEreRYKOFCpZQ75n56Sl5/t89JQik9xtOYtgn3294tScviMqYWbhzme/LOF2gdZ6vhA3gHc0WfVIZrcS9X7umArdpqK9DB9iECwHYV/gSw/j6vv4pJ5v0wjbhm2JY5NY40yIpJZkkTWK9Ubs8LU30a/0FXSj3D7arHpDomedwB89n4M8nyZipJq39MrUDS9ZKzHmiulMqZ4PU5MX2e1ydY1hXTEDAA0yCpGqZRmiN5cWfZ905mvWhMY637uXU8nZMs97f9uyHloSWmtY7QWi8FfsBc1SV7ReWgc/lPzHnaLsnypMeZVisxV9K5tdZbkvlLmtzt6aTtr0OS5a1J2cXgn0AdpVTC2z45MOd5Qik9xk1AKZW4V4nizvcgWbbzYCNQFdiW3L5sq+7D3LN/+C7x/JuS/aUncuWefryHqY5dqpT6DvAGPsS0lP4iBa//jtv33ZOrkndWHA9Ma31FKbUJeEMpdRpT/deHO6u7U2olppX+FKXUZMy99ve48wr7lodt682yrRuIaeS1Nsl61zHjE3yG+UHyIabldEqqzj/CVO+tVUp9gvnh8TYmsY4AUEqVxtxi+VFrPce2bBgwWim1Wmv9awr2kxorgZ5KqZ2YRkVtMQkwqd1Af6VUJ8wVebjWel/SlbTW/ymlZgIf2O5b/oG5F/seMPNBvyiVUiMwV2q/Yq7oSmBay+/Qpj/33dj1XNZa/6KU+g34QSlVANNYrxOmAaHdaK1DbOU3Vyk1BtiMufr0xbT0fltrvd+e+0yw7zil1IfABKXURMy96jLAEEy53W8cgC8xretXK9OtLxp4C/PZSbiflB7jFMznZL7ts3Ae0wjx1q2slIxLMBjzg3KVUupHTG1AAUybDHet9RCttVZKDQAW2WruZmO+hwpjPgvHtdZjUrCvdEOu3NMJrfVKzP3OPJgTazymNeeTKfmlbvvC3A9s0VpvsyqONOiCudIah/lAn8GMGvbAtNarMF/+9YClmB8KPTCJKzmvYe7P/ozpmrQU02UmqamYZPEtpjX7eUyL+uTaMiSN6V9MdfJV22unYe431tda/2NLgjO487g/x7T2naaS9L23o1cwDfsCMWWQE/N+JPWJLZaJmEZuP9xjmz1t6/fBtH94wfa4Zyri+wvzpf8lplvbJ9hqQu71Igedy20xx/Mxpqw8MPdu7a0bph1Ee0xjtrm2/RzA/vf4E9FaT8R0pXzOtu8XuN2K/cp9XnsB0wXtAuY8H4f58Zh0VENIwTHa2r80Bv7FvH8/YVrzj7Nt457x2LaxDdPmJgzTpmc1ZoTHR0lQi6S1Xo5pOJcDc46vAj7F3MZJSe1cuqJS3shSpGdKqQqY1tF9tdY/Wh1PRqDMWOC/As9prZO7z5ZwXY0ZXexdJ4QmRLqizMBKmzFd16alg3iWApW01mWtjiW9kmr5DE6ZkbzKYaobT2Ou/oQQIlVst4cGcHvAnkrAUMzATHfrgubIeAZjarkOYGqVOmBqZF52diwZiST3jO9FzAhP+zEDUVy/z/pCCHEv1zHtCHpg7m1fwrQgH6K1jrQgnmjMbYJSmAaN+4AXpYby3qRaXgghhHAx0qBOCCGEcDGS3IUQQggX4zL33AsUKKB9fX3ttr1r166RI0eO+68oHoiUq/1JmTqGlKv9SZna39atWy9orQsmXe4yyd3X15ctW7bcf8UUCgkJoUGDBnbbnjCkXO1PytQxpFztT8rU/pRSx5JbLtXyQgghhIuR5C6EEEK4GEnuQgghhIuR5C6EEEK4GJdpUHc3cXFxXLhwgcuXLxMbG5vi1+XOnZs9e/Y4MLLMKXfu3Ozfv588efJQoEAB3Nzk96UQQtib05O7UmoSZurLc1rrO6ZKtM3V+xVm2r9IoFdaZjk7efIkSil8fX3x9PTEbP7+wsPDyZkz5/1XFA/k6tWrZMuWjbNnz3Ly5ElKlbrbFOtCCCFSy4rLpilA03s87weUt/0FYOa3TrVr165RvHhxsmTJkuLELhxHKUWWLFkoXrw4165dszocIYRwSU5P7lrrDcC95r9uBUzVxiYgj1KqaFr2KVW/6Y+8J0II4Tjp8Z57ceBEgscnbctOJ11RKRWAubqncOHChISE3LGx3LlzEx4e/sBBxMbGpup14t4SlmtUVFSy75l4MBEREVKODiDlan9Sps6THpN7cnXnyU5dp7UOAoIAatWqpZMb+WjPnj2puncu99wdI2G5ZsuWjerVq1scUcYno345hpSr/UmZOk96rBs9CZRM8LgEEGpRLE4xY8YMatWqhbe3N0WLFsXPz4/ffvvNIfsKCQmhRIkSdtvelClTePLJJ+22PSGEEGmXHpP7YqCHMuoAV7TWd1TJu4oxY8YwaNAghg4dytmzZzl+/Dj9+/dn0aJFVocmhBAig3J6cldKzQT+BCoqpU4qpV5QSr2klHrJtspy4DBwEJgA9Hd2jM5y5coV3n//fcaNG0fbtm3JkSMHnp6etGzZks8++4zo6GgGDRpEsWLFKFasGIMGDSI6Ohq4fQX+xRdfUKhQIYoWLcrkyZPjt718+XIqV65Mzpw5KV68OJ9//jnXrl3Dz8+P0NBQvL298fb2JjQ0lM2bN1O3bl3y5MlD0aJFGThwIDdu3IjfllKK8ePHU758efLmzcuAAQPQWrNnzx5eeukl/vzzT7y9vcmTJ4+zi1AIIUQynH7PXWvd5T7Pa2CAo/Y/aNAgduzYcd/1YmNjcXd3T9U+qlWrxtixY++73p9//klUVBRt2rRJ9vnAwEA2bdrEjh07UErRqlUrRo4cyUcffQTAmTNnuHLlCqdOneKXX36hffv2tG7dmrx58/LCCy8we/ZsnnrqKS5dusSRI0fIkSMHK1asoFu3bpw8eTJ+P6dPn+bLL7+kVq1anDx5Ej8/P7777jsGDRoUv87SpUv5+++/uXr1KjVr1qRly5Y0bdqU8ePHM3HiRIfdRhBCCPHg0mO1fKYRFhZGgQIF8PBI/jdWcHAw77//PoUKFaJgwYIMHz6cadOmxT/v6enJ+++/j6enJ82aNcPb25t9+/bFP7d7926uXr1K3rx5qVGjxl3jqFmzJnXq1MHDwwNfX1/69evH+vXrE60zZMgQ8uTJQ6lSpWjYsGGKfiAJIYSwRnpsLe9QKbmiBue0ls+fPz8XLlwgJiYm2QQfGhqKj49P/GMfHx9CQ0MTvT7h67y8vIiIiABg3rx5jBw5kiFDhlClShVGjx5N3bp1k41j//79DB48mC1bthAZGUlMTAw1a9ZMtE6RIkWS3Y8QQoj0R67cLVS3bl2yZcvGwoULk32+WLFiHDt2LP7x8ePHKVasWIq2/dhjj7Fo0SLOnTtH69at6dixI0Cyo/S9/PLLPPTQQxw4cICrV68yatQozN2R+5NR/4QQIv2R5G6h3LlzM2LECAYMGMDChQuJjIzk5s2brFixgv/973906dKFkSNHcv78eS5cuMCIESPo1q3bfbd748YNgoODuXLlCp6enuTKlSu+/UDhwoUJCwvjypUr8euHh4eTK1cuvL292bt3L99/n/IRfwsXLszJkycTNcATQghhrUxXLZ/eDB48mMKFCzNy5Ej8/f3JmTMnNWvWZNiwYdSoUYOrV69SpUoVADp06MC7776bou1OmzaNgQMHEhsbS8WKFZk+fToADz30EF26dKFMmTLExsaye/duPv/8cwICAvj000+pXr06nTp1Yt26dSnazzPPPMPDDz9MkSJFcHNz48KFC6krCCGEEHajUlr9mt7VqlVLb9my5Y7le/bsoVKlSg+8PRmhzjESlmtq3xuRmIz65RhSrvYnZWp/SqmtWutaSZdLtbwQQgjhYiS5CyGEEC5GkrsQQgjhYiS5CyGEEPeTwaYAl+QuhBBC3Mvy5VC6NPzzj9WRpJgkdyGEEOJuoqPhtdegQAHIQL17pJ+7EEIIcTdjxsDBg7BqFWTJYnU0KSZX7kIIIURyTp6EkSOhTRto3NjqaB6IJPd0aOPGjVSsWNHqMIQQInN7802IizNX7xmMJHeL+fr6smbNmkTLnnrqqfipW4UQQljg11/h55/hnXfA19fqaB6YJHcRLyYmxuoQhBDCejdvwiuvmBbyb71ldTSpIsk9HQoJCaFEiRLxj319ffn888+pUqUKuXPnplOnTkRFRcU/v3TpUqpVq0aePHl44okn+Pfff+OfGz16NGXLliVnzpxUrlyZBQsWxD83ZcoU6tWrx+uvv06+fPn44IMPnHJ8QgiRrn37Lfz3H3z5JWTPbnU0qZLpWssPGgQ7dtx/vdjY7NhmSX1g1arB2LGpe+3dzJ49m5UrV5ItWzbq1avHlClTeOmll9i2bRt9+vRhyZIl1KpVi+nTp/P888+zb98+smbNStmyZdm4cSNFihRhzpw5dOvWjYMHD1K0aFEA/vrrLzp37sy5c+e4efOmfYMWQoiMJjQUhg+H5s3h+eetjibV5Mo9g3j11VcpVqwY+fLlo2XLluyw/UKZMGEC/fr14/HHH8fd3Z2ePXuSNWtWNm3aBJhpYosVK4abmxudOnWifPnybN68OX67xYoV45VXXsHDw4PsGfQXqhBC2M2bb8KNG/DVV6CU1dGkWqa7ck/pFXV4+PV0NeVrkSJF4v/v5eVFaGgoAMeOHeOnn37im2++iX/+xo0b8c9PnTqVMWPGcPToUQAiIiISzblesmRJJ0QvhBAZwK+/wsyZ5sq9bFmro0mTTJfcXU3JkiUZNmwYw4YNu+O5Y8eO0bdvX9auXUvdunVxd3enWrVqaK3j11EZ+JepEELYzY0bMGCAaUT39ttWR5NmUi2fDty8eZOoqKj4vwdptd63b1/Gjx/PX3/9hdaaa9eusWzZMsLDw7l27RpKKQoWLAjA5MmT2bVrl6MOQwghMq6vvoI9e+DrrzNsI7qEJLmnA82aNSN79uzxfw/Sar1WrVpMmDCBgQMHkjdvXsqVK8eUKVMAqFy5Mm+88QZ169alcOHC7Ny5k3r16jnmIIQQThccHIyvry9ubm74+voSHBxsdUgZ04kT8OGH0KKF+XMFWmuX+KtZs6ZOzu7du5Ndfj9Xr15N1evEvSUs19S+NyKxX3/91eoQXFJ6L9fp06drLy8vDcT/eXl56enTp1sd2l2l2zJt21br7Nm1PnLE6kgeGLBFJ5MT5cpdCCEyoGHDhhEZGZloWWRkZLLtb8Q9LF8O8+fDe+9lyJHo7kaSuxBCZEDHjx9/oOUiGZGRMHAgPPQQvPGG1dHYlSR3IYTIgEqVKvVAy0UyRo2CI0fgu+8y1HSuKSHJXQghMqDAwEC8vLwSLfPy8iIwMNCiiDKYvXvh00+hWzdo2NDqaOxOkrsQQmRA/v7+BAUF4ePjg1IKHx8fgoKC8Pf3tzq09E9r6N8fvLzg88+tjsYhZBAbIYTIoPz9/SWZp8a0aWY0uu+/h8KFrY7GIeTKXQghROYRFmYaz9WpAwEBVkfjMJLchRBCZB5vvw2XLsEPP4Cb66ZA1z0ycYejR4+ilHqg4W3vZ9SoUbz44ot2254QQjjMxo3w448weDBUqWJ1NA4lyT0dmDJlCo8++iheXl4UKVKE/v37c+XKFavDukNISAglSpRItGzo0KFMnDjRooiEECKFbtyAfv3Ax8fM+ubiJLlb7IsvvuDtt9/ms88+48qVK2zatImjR4/SuHFjbt686bQ4tNbExcU5bX9CCOFUn39uJoYZNw5y5LA6GoeT5G6hq1evMnz4cL755huaNm2Kp6cnvr6+zJ49myNHjjBjxgx69erFu+++G/+apFfPo0ePpmzZsuTMmZPKlSuzYMGC+OdiY2N58803KVCgAGXKlGHZsmWJ9t+gQQOGDRtGvXr18PLy4vDhw0yePJlKlSqRM2dOypQpww8//ADAtWvX8PPzIzQ0FG9vb7y9vQkNDeWDDz6gW7du8dv87bffeOKJJ8iTJw8lS5aMn8RGCCEsc+AAjBgBHTpA8+ZWR+MUma8r3KBBsGPHfVfLHhsL7u6p20e1ajB27H1X++OPP4iKiqJt27aJlnt7e+Pn58fq1avx9PS85zbKli3Lxo0bKVKkCHPmzKFbt24cPHiQokWLMmHCBJYuXcr27dvJkSMH7dq1u+P106ZNY8WKFVSsWBGtNfv27WPp0qWUKVOGDRs24Ofnx2OPPUaNGjVYsWIF3bp14+TJk8nGcvz4cfz8/AgKCqJ9+/ZcvXqVEydO3LcchBDCYbQ21fHZsplpXTMJuXK30IULFyhQoAAeHnf+xipatCjnz5+/7zY6dOhAsWLFcHNzo1OnTpQvX57NmzcDMHv2bAYNGkTJkiXJly8f77zzzh2v79WrFw8//DAeHh54enrSvHlzypYti1KK+vXr07hxYzZu3Jii4wkODubZZ5+lS5cueHp6kj9/fqpVq5ai1wohhENMmWL6tH/6KRQtanU0TpP5rtxTcEUNcD08nJw5czo0lAIFCnDhwgViYmLuSPCnT5+mYMGC993G1KlTGTNmDEePHgUgIiKCCxcuABAaGkrJkiXj1/Xx8bnj9QmfB1ixYgUffvgh+/fvJy4ujsjISB599NEUHc+JEycoW7ZsitYVQgiHO3fO9Gl/8klI0qtHa01oeCiHLh3i0MVDHLp0iLMRZwm7HsbF6xe5FHWJm7E30WjidBxaa1Z3X41vHl9rjuUBZb7kno7UrVuXrFmzMn/+fDp27Bi//Nq1a6xYsYKRI0fyzz//JJrW8cyZM/H/P3bsGH379mXt2rXUrVsXd3d3qlWrhpni11z9J6wWT262KKVU/P+jo6Np164dU6dOpVWrVnh6etK6dev47SVcNzklS5aMrzUQQgjLvf46RERAUBCh186w4dgGtoZuZevprWw7vY0r0bd7JbkrdwrlKER+r/zky56P0nlK4+nuiZtyQ6FQSpHVPauFB/NgJLlbKHfu3AwfPpxXXnmFXLly0ahRI06dOkX//v0pUKAA/v7+KKX44osvePfdd7lx4wZjE9Q8XLt2DaVU/BX+5MmT2bVrV/zzHTt25Ouvv6ZFixbkyJGD0aNH3zOeGzduEB0dTcGCBfHw8GDFihWsXr2aRx55BIDChQsTFhbGlStXyJ079x2v9/f3Z9SoUcyePZu2bdty5coVTpw4IVXzQgini1i2gF+2zWDN0FqsW9eWvRf2ApDFPQtVCleh8yOdqVK4CuXylaNs3rKUyl0KT/d7t3HKSCS5W+x///sf+fPn58033+TgwYNER0dTv3591qxZQ44cOejevTtr1qzB19cXX19fevfuzRdffAFA5cqVeeONN6hbty5ubm706NGDevXqxW+7b9++7N+/n6pVq5IrVy7efPNN1q1bd9dYcubMyddff03Hjh2Jjo6mZcuWPP/88/HPP/TQQ3Tp0oUyZcoQGxvL7t27E72+VKlSLF++nDfffJMXX3yR3LlzM3LkSEnuQginuHj9Ikv2LWH+rtms3reCqM6Qw3MPT+d5mheqv0BD34Y8WvhRsri71vSuyVG3qlwzulq1auktW7bcsXzPnj1UqlTpgbcX7oR77smZNGkSw4cP5/fff3fJeZkTlmtq3xuRWEhICA0aNLA6DJcj5Wp/jijTm7E3WXVoFVN2TGHxvsXcjLtJiVhv2m6JoE3/b3jCL8Clk7lSaqvWulbS5XLlns706dMHT09P/vjjD5dM7kIIYQ/Hrxznu7+/Y8qOKZy9dpaCXgUZWHsgXfTD1Gr2ImrAQGgx0OowLSPJPR3q3r271SEIIUS6o7Vm4/GNfP3X1yzYawbsalGhBX2q9aFZ+WZ4xsRB9epQoiSMGmVxtNaS5C6EECJd01qzdP9SAjcG8tepv8ibLS9v1n2T/o/1xydPgi6+H31ghphdtgwsuK2ankhyF0IIkS7F6Tjm7Z5H4MZA/jn7D755fPmu2Xf0rNYTL0+vxCvv3Gmu1rt2hWbNrAk4HZHkLoQQIl3RWrP60GreWfsO289sp2L+ivzU+ie6PNIl+e5qMTHQpw/kyZOphpi9F0nuQggh0o3Npzbz9pq3CTkaQuk8pZnWZhpdHumCu9s95vr48kvYsgV+/hkKFHBesOmYJHchhBCWCw0PZciaIUz7dxoFvQryddOv6Ver3/27se3bB++9B23amFnfBCDJXQghhIWiY6IZu2ksIzeO5EbsDYbUG8LQp4aSM2sKGsTFxZkx47NnN/O032eI7MxEZoXLwHx9fVmzZk2yz23cuJGKFSs6OSIhhEi59UfXU3V8VYasHUJD34b81/8/Pn7245QldoDvvoPffjMTgmWiGd9SQpJ7OjBjxgxq1aqFt7c3RYsWxc/Pj99++y1N23zqqafYt2+fnSIUQgj7uXj9Ii8seoEGPzXgRuwNlnddzuIuiymXr1zKN3L4MLz9NjRtCj16OCzWjEqq5S02ZswYRo8ezfjx42nSpAlZsmRh5cqVLFq0iCeffNLq8IQQwq7m/DeHgSsGEhYZxv+e+B/DGwy/s1vb/cTFmdbxHh4wYYJUxyfDkit3pVRTpdQ+pdRBpdSQZJ7PrZRaopT6Ryn1n1KqtxVxOtqVK1d4//33GTduHG3btiVHjhx4enrSsmVLPvvsM3r16sW7774bv35ISAglSpRItI2///6bypUrkzdvXnr37k1UVFSy6544cYK2bdtSsGBB8ufPz8CBmXdYRiGE84VFhjFi9wg6zu1IyVwl2RKwhU+e++TBEzuY6vj1600r+STfiWkVHByMr68vbm5u+Pr6EhwcbNftO4vTr9yVUu7AOOA54CTwt1JqsdY64RRjA4DdWuuWSqmCwD6lVLDW+kZa9z9o5SB2nNlx3/ViY2Nxd79H14t7qFakGmObjr3ven/++SdRUVG0adMmVfsBcyKuWrWKHDly0LJlS0aOHMnIkSMTrRMbG0uLFi145plnmDZtGu7u7iQ3yY4QQjjC4n2LCVgSQFhkGB81/IghTw7Bwy2V6SdhdXxv+173BQcHExAQQGRkJADHjh0jICAAMFNaZyRWXLnXBg5qrQ/bkvUsoFWSdTSQUymlAG/gIhDj3DAdLywsjAIFCuDhkfrfWAMHDqRkyZLky5ePYcOGMXPmzDvW2bx5M6GhoXz22WfkyJGDbNmySZW/EMLhrt24RsCSAFrNakUR7yJ8X+N73n363dQndgdXxw8bNiw+sd8SGRnJsGHD7LofZ7Dinntx4ESCxyeBx5Os8y2wGAgFcgKdtNZxSTeklAoAAgAKFy5MSEjIHTvLnTs34eHh8Y8/qvdRioJMy5U7kGifd5M9e3YuXLjApUuXkk3wN2/eJDo6On5bkZGRaK3jH2utyZ8/f/zjAgUKEBoaSnh4eKJ19+/fT8mSJbl+/Xqqj8deYmNj4+ONiopK9j0TDyYiIkLK0QGkXNPmYMRBPtrzESciT9C5ZGf6+PYhOjI6TWVafMECyq9fz9633uLMwYNw8KD9AgaOHz9+1+UZ7VywIrkn91Mr6aTyTYAdwDNAWeAXpdRGrfXVRC/SOggIAjOfe3LzBO/ZsydV87I7Yz73Ro0akS1bNtauXUv79u3veD5PnjzExsbGx3H16lWUUvGPlVKEhYXFPw4LC6NYsWLkzJkTLy+v+HUrVKjAyZMnyZ49e5pqCewhYblmy5aN6tWrWxqPK5B5xx1DyjV1tNZ89ddXvP3b2xTwKsCaHmt4pvQzQBrL9MABc7Xu58dDn3zCQw5oRFeqVCmOHTuW7PKMdi5YUS1/EiiZ4HEJzBV6Qr2B+do4CBwBHnJSfE6TO3duRowYwYABA1i4cCGRkZHcvHmTFStW8L///Y9q1aqxfPlyLl68yJkzZxg7duwd2xg3bhwnT57k4sWLjBo1ik6dOt2xTu3atSlatChDhgzh2rVrREVF8fvvvzvhCIUQmcnlqMu0+bkNr696Hb9yfvz70r/xiT1NYmOhVy/Ilg0mTnRY6/jAwEC8vBI38PPy8iIwMNAh+3MkK5L730B5pVRppVQWoDOmCj6h40AjAKVUYaAicNipUTrJ4MGDGTNmDCNHjqRgwYKULFmSb7/9ltatW9O9e3eqVq2Kr68vjRs3TjZxd+3alcaNG1OmTBnKlCmTqHX9Le7u7ixZsoSDBw9SqlQpSpQowc8//+yMwxNCZBJbQ7dS44caLDuwjLFNxrKg0wLye+W3z8bHjIE//oBvv4VixeyzzWT4+/sTFBSEj48PSil8fHwICgrKcI3pAJTWSWvEnbBTpZoBYwF3YJLWOlAp9RKA1nq8UqoYMAUoiqnGH621nn6vbdaqVUsn1wJ8z549VKpU6YFjdEa1fGaUsFxT+96IxKT62DGkXFNGa03Q1iBeXfkqhXMUZnaH2dQpUSfZdVNVpv/9BzVqQIsWMHeu9GlPQim1VWtdK+lyS27Aaq2XA8uTLBuf4P+hQGNnxyWEECLlomOiGbh8IBO3T6RpuaZMazONAl52nJXt5k0z+lzu3PD995LYH4CMUCeEEOKBhYaH0m52Ozad3MTQJ4cyouGIe0/LmhojR8K2bTBvHhQqZN9tuzhJ7kIIIR7IppObaPNzG8Kjw5nbYS7tKrez/042b4bAQHPl3rat/bfv4mTiGCGEECkW/G8wDaY0IIdnDja9uMkxiT0yErp3h+LF4euv7b/9TECu3IUQQtxXnI7j/V/fJ3BjIPV96jOv4zz7tYZP6u23Yf9+WLfO3G8XD0ySuxBCiHuKvBlJjwU9mLdnHi9Wf5FxzceRxT2LY3b2yy+my9ugQdCwoWP2kQlIchdCCHFX566do+XMlvx96m/GNB7DoDqDUI5qtX7xopkMplIlGDXKMfvIJCS5CyGESNa+C/toNqMZp8NPs6DTAlo9lHSOLzvSGl56Cc6dg8WLIXt2x+0rE5AGdZnUlClTZGY4IcRd/Xb8N56Y9ATh0eGE9ApxbGIHmDYN5syBjz4yg9aINJHkLu7rgw8+oFu3blaHIYRwkgV7FvDs1Gcp4FWATS9uonbx2o7d4ZEjMHAgPP00vPmmY/eVSUhyv4fg4GB8fX1xc3PD19eX4OBgh+4vJsblpqwXQmQwQVuDaD+nPdWKVOOPPn9QJm8Zx+4wNtZ0e1MKpk6FNEy1LW6T5H4Xs2fPJiAggGPHjqG15tixYwQEBNg9wfv6+vLJJ59QpUoVcuTIQUxMDKNHj6Zs2bLkzJmTypUrs2DBgvj1fXx82Lp1KwDTp09HKcXu3bsBmDhxIq1bt052P2FhYTz//PPkypWL2rVrc+jQoUTPv/baa5QsWZJcuXJRs2ZNNm7cCMDKlSsZNWoUP//8M97e3lStWhWAyZMnU6lSJXLmzEmZMmX44Ycf7FouQgjn0lrz0fqP6Le0H03LNWVtj7WO6+qW0CefwO+/w3ffgY+P4/eXSUhyv4sPP/yQyMjIRMsiIyMZNmyY3fc1c+ZMli1bxuXLl/Hw8KBs2bJs3LiRK1euMHz4cLp168bp06cBqF+/PiEhIQBs2LCBMmXKsH79+vjH9evXT3YfAwYMIFu2bJw+fZpJkyYxadKkRM8/9thj7Nixg4sXL9K1a1c6dOhAVFQUTZs2ZejQoXTq1ImIiAj++ecfAAoVKsTSpUu5evUqkydP5vXXX2fbtm12LxshhOPF6TheWfEK74e8T4+qPVjYaSE5suRw/I7/+gvefx86d4auXR2/v0xEkvtdnDx5Mtnlx48ft/u+Xn31VUqWLEl2W+vQDh06UKxYMdzc3OjUqRPly5dn8+bNgEnut5L5xo0beeedd+Ifr1+/PtnkHhsby7x58xgxYgQ5cuTgkUceoWfPnonW6datG/nz58fDw4M33niD6Oho9u3bd9eYmzdvTtmyZVFKUb9+fRo3bhx/tS+EyDhi4mLoubAn4/4exxt132BKqyl4uns6fsfh4Sahlyghk8I4gCT3uyhRokSyy0uVKmX3fZUsWTLR46lTp1KtWjXy5MlDnjx52LVrFxcuXABMct+4cSNnzpwhNjaWTp068fvvv3P06FGuXLlCtWrV7tj++fPniYmJSbQfnyTVX1988QWVKlUid+7c5MmThytXrsTvMzkrVqygTp065MuXjzx58rB8+fJ7ri+ESH+iYqJoP7s90/+dTuAzgXz23GeO68Oe1MCBcPQoTJ8OefI4Z5+ZiCT3uxg+fDheXl6Jlnl5eREYGGj3fSX8MB07doy+ffvy7bffEhYWxuXLl3nkkUfQWgNQrlw5vLy8+Prrr3n66afJmTMnRYoUISgoiCeffBI3tzvf0oIFC+Lh4cGJEyfilyWsgdi4cSOffPIJs2fP5tKlS1y+fJncuXPH7zPphz06Opp27drx5ptvcvbsWS5fvkyzZs3i1xdCpH8RNyJoMaMFi/Yt4lu/bxn61FDnJfaZM03juXffBemS6xCS3O+iY8eOBAUF4ePjg1IKHx8fgoKC8Pf3d+h+r127hlKKggULAqbh2q5duxKtU79+fb799tv4KvgGDRokepyUu7s7bdu25YMPPiAyMpLdu3fz008/xT8fHh6Oh4cHBQsWJCYmhhEjRnD16tX45wsXLszRo0eJi4sD4MaNG0RHR8f/aFixYgWrV6+2azkIIRznStQVmkxvQsjREKa2nsqA2gOct/OjR81gNXXrwnvvOW+/aRQRYXUED0aS+z34+/vHJ7WjR486PLEDVK5cmTfeeIO6detSuHBhdu7cSb169RKtU79+fcLDw3n66aeTfZycb7/9loiICIoUKUKvXr3o3bt3/HNNmjTBz8+PChUq4OPjQ7Zs2RJV4Xfo0AGA/PnzU6NGDXLmzMnXX39Nx44dyZs3LzNmzOD555+3ZzEIIRzk4vWLPDvtWf4+9TezO8yme9XuTtu3iomBLl3Mg+Bg8Ejfg6TGxsKKFdChg5lOPjTU6ohSTrlKVWqtWrX0li1b7li+Z88eKlWq9MDbCw8PJ2fOnPYITSSQsFxT+96IxEJCQmjQoIHVYbgcVyzX89fO03h6Y3af3828jvNoUaGFU/d/zN8fnxkzYNYs6NTJqft+EPv3w5Qp8NNPJqHnz2+64r/1FhQrZnV0iSmltmqtayVdnr5/NgkhhLCLsxFnaTS1EYcuHWJJlyU0LtvYuQGsXUupmTPhhRfSZWIPDzej306eDL/9Bm5u4OcH33wDLVpAFgdNgucoktyFEMLFnYk4Q8OfGnL8ynGWd11Ow9JOnkr13Dno1o3IkiXJ8dVXzt33PWhtEvmkSSaxX7sGFSvC6NHmSj29XaU/CEnuQgjhwm4l9hNXTrDCfwVP+9y9bY5DxMVBr15w6RK7x43jsRxOGBznPk6eNI31J0+GgwfB29s0Bejd27Tzc4Uu95LchRDCRZ0OP80zU5/hxJUTLPdf7vzEDjBmjGmV9u23XCtb1vn7t4mOhkWLTEJfvdr85nj6adMbr317SAe/OewqUyT3uLi4ZPt/C+vc6lYnhHCMMxFn4hP7Cv8VPOXzlPOD2LQJ3nkH2raF/v3BNpqmM23fbhJ6cDBcvAglS8LQoaYywcLfGg7n8sk9R44cnDp1isKFC+Pp6em8QRpEsrTW3Lhxg7Nnz5LD1X4qC5FOnL92nkZTG3H8ynFW+q+0JrFfvGgazpUsCT/+6NS67rAwk8wnT4YdOyBrVmjd2lS7P/ts5ph4zuWTe4kSJbhw4QLHjh17oClVo6KiyJYtmwMjy5yioqLw9vYmd+7cFChQwOpwhHA5YZFhPDvtWY5cOsJy/+XWJHatTSY9fdrM+OaE4WVjY011+6RJpvr95k2oUQO+/dbcT8+Xz+EhpCsun9zd3NwoVKgQhQoVeqDXhYSEUL16dQdFlXlJuQrhOJeuX+K5ac+x78I+lnRZQgPfBtYE8tVXsHgxfPklPPaYQ3d14IC5Qr/VJ71AARgwwPy2qFLFobtO11w+uQshRGZwNfoqTYOb8t/5/1jYaSHPlX3OmkD++gv+9z9o1Qpee80hu4iIMF3XJk1yjT7pjiDJXQghMrjIm5G0nNmSraFbmddxHn7l/awJJCzMjNVavLi5nLbjffZbfdInT4bZsxP3Se/RA4oWtduuXIIkdyGEyMCiY6Jp83MbNh7byIx2M2j1UCtrAomLMyO/nD1r7rPnzWuXzZ46ZfqkT5qUuE96nz5Qp45r9El3BEnuQgiRQd2MvUmnuZ1YfWg1Pz7/I50f6WxdMJ98YvqzjxsHte4Y6vyBREebW/aTJt3uk16/vplErl071+uT7giS3IUQIgOK03H0XtSbRfsW8Y3fN/Sp3se6YEJCzGgwnTvDyy+nejOZtU+6I0hyF0KIDEZrzSvLXyF4ZzCjnhnFwNoDrQsmNNQk9XLlICjogevJ79YnvU8faNQoc/RJdwRJ7kIIkcG89+t7fLflO9564i2GPDnEukBu3oSOHc2UamvXQgqnyU7YJ33xYrhxA2rWNDX6XbrY7XZ9pibJXQghMpDP//icwI2B9K3Rl0+e/cTaUTf/9z/TeG7GDHj44fuufvJkdoYONQ3kTp0yfdL795c+6Y4gyV0IITKISdsn8dYvb9Hp4U583/x7axP77Nkwdiy88oq53L6LxPOkPx7fJ/3rr6VPuiNJchdCiAxg4d6F9F3SlyZlmzC1zVTc3Sy8Gb1nj7kpXrcufP75HU/fbZ70gIBDDB9eNkPPk55RSHIXQoh0LuRoCJ3nduaxYo8xr+M8srhbeLl79Sq0aWP6o82Zk+jS+9Y86VOmmGFhvb1NW7veveGJJ2D9+hMUKybN3p1BkrsQQqRj209v5/mZz1MmbxmWdV1GjiwWdvKOizPDwR08aBrQFS+ebJ/0p582Xdg6dJA+6VaR5C6EEOnUwYsHaRrclLzZ87K6+2rye+W3NqCPPzZTro0dy47c9Zn06u0+6SVKSJ/09ESSuxBCpENnI87SZHoTYuNiWdVtFSVylbA2oBUr0O+9x76a/nSZ/Co7Bkmf9PRMkrsQQqQz4dHhNJvRjDMRZ1jXYx0PFXjIslhiY+G3nw5R86WuHKQqT2wNolINxTffQNeumW+e9IxCkrsQQqQjN2Jv0HZ2W/458w9Luizh8RKPWxLHrXnS500JZ97pVtxQbizpNZ9Ng7ykT3oGIMldCCHSiTgdR6+FvVhzeA1TWk1x+tStSedJd1dxbCzUg8pue4ldtor3mpZ2ajwi9dysDkAIIYTx9i9vM3PXTD5u9DE9q/V0yj61ho0bzX3zIkXMv+fPm3nSLw7+iLpnF+L2xed4Nm3klHiEfciVuxBCpANfbfqKz//8nIGPDeTtem87fH/JzZPeufPtsWnUwgUw5APo2RNee83h8Qj7kuQuhBAWm/PfHF5f9TptHmrD2KZjHTasbIrnSd+50/Rnr10bxo9/4JnehPUkuQshhIU2HNtAtwXdeKLkEwS3DXbIsLIPNE/6+fPw/PNmhrf58yFbNrvHIxxPkrsQQlhkz/k9tJrVijJ5y7C4y2Kye2a327ZvzZM+aRL8808K+6TfuGEu4c+cgfXroXhxu8UjnEuSuxBCWOB0+Gn8gv3I6p6VFf4ryJc97R3G0zRPutZm/tWNG80UrrVrpzkeYR1J7kII4WQRNyJoMbMF5yPPs6HXBnzz+KZpe/v3m8lafvoJQkMhf/5UzJP+9dfw44+mvv4eU7iKjEGSuxBCOFFMXAwd53Rkx5kdLO68mJrFaqZqO4nnSSfRPOktWz7gPOnLlsHgwabe/qOPUhWPSF8kuQshhJNorRmwbAArDq7ghxY/0LxC8wd8vUnkkyfD7Nm350kfPRq6dyd186Tv3Gn6wFWtCtOnm18JIsOT5C6EEE7y2R+fEbQtiCH1hhBQMyDFr0tunvQuXUy1e926aeipdvYstGgBuXLBkiUyP6sLkeQuhBBOMPu/2by95m06PdyJwEaB913/bn3Shw2D9u3tkIevX4dWreDCBdOITlrGuxRJ7kII4WC/H/+dHgt6UK9kPaa0noKbunvV944dJqE7dJ70uDizwc2bTV/2GjXstGGRXliS3JVSTYGvAHdgotZ6dDLrNADGAp7ABa11fSeGKIQQdnHw4kFazWpFqdylWNh5Idk87hwU5laf9MmTTXJ3+DzpQ4eam/affWZ2JFyO05O7UsodGAc8B5wE/lZKLdZa706wTh7gO6Cp1vq4UqqQs+O0h3PXzlEoR4YMXQhhB2GRYTQLbgbAcv/lFPAqEP/c3fqkf/utuZ/usHnSf/gBPvkEXn4Z3njDQTsRVrOiWWRt4KDW+rDW+gYwC2iVZJ2uwHyt9XEArfU5J8eYZj/t+IliXxTj6OWjVocihLBAdEw0bX5uw7Erx1jUeRHl8pUDTIO4oUPBxweaNYNffzV59p9/YMsWGDDAgYl95Uqzg2bNTJ85GTPeZSmttXN3qFR7zBX5i7bH3YHHtdYDE6wzFlMd/zCQE/hKaz01mW0FAAEAhQsXrjlr1iy7xRkREYG3t3eqX38m6gxd/+pKT5+e9PR1ztSNGUFay1XcScrUMdJSrlprAvcGsvbcWt6r9B51cz5HSEhBVqwows6deXBz09SufRE/v9PUrRuGp6fjv4e9Dx6k2quvcr14cXZ8/TWx2e031G1Kyblqfw0bNtyqta6VdLkV99yT+6mY9Mz2AGoCjYDswJ9KqU1a6/2JXqR1EBAEUKtWLd2gQQO7BRkSEkJatzfp/CR+DfuVH+v/eM8GNJmJPcpVJCZl6hhpKdfhvw5n7bm19C0TyMmVQ+lg65NeoQJ8/DH06KEoViw/kN+uMd/VsWPQtSsUKEDOkBCesqhlvJyrzmNFxjkJlEzwuAQQmsw6K7XW17TWF4ANQFUnxWc3far34diVY6w7ss7qUIQQTvJVyDRGbBhBrkO9mdDjHebMMWPE/P477N0LQ4akcrCZ1Lp40QxdFxkJK1ZIl7dMwork/jdQXilVWimVBegMLE6yziLgKaWUh1LKC3gc2OPkONOs9UOtyZstL5O2T7I6FCGEA0VHm6FgH++4gUHrXoAjDal6YjxTpijOnIGJE+GJJyy4xR0VZVrDHzoEixbBww87OQBhFadXy2utY5RSA4FVmK5wk7TW/ymlXrI9P15rvUcptRL4F4jDdJfb5exY0yqbRzb8H/VnwrYJXLp+ibzZ7zUlkxAio0nUJ52DuAW0Ib9bGVb/bx41Kj/I4O4OEBsLPXqYAWpmzTIj4IhMw5J+7lrr5cDyJMvGJ3n8GfCZM+NyhD7V+/Dt398yY+cMBtQeYHU4Qog0CgszM6JOmnS7T7pf24v8XbU5UUqx6cWllMtn8Q95reG110x1wuefQ6dO1sYjnE5aeTlY9aLVqVakGpN2SNW8EBlVbKy5Xd2xo7lf/uqrZmCZb7+FoyducLlJO87fPMrCzgvju7xZKjDQTOL+5pvSlz2TkuTuBC9Uf4Ftp7ex48wOq0MRQjyA/ftNn/RSpUzX8HXrzDzpt/qk9++vGfrHS4QcDeHH53/kyVJPWh0yTJgA771npon75BOroxEWkeTuBF0f7UpW96xM3j7Z6lCEEPcRHm6q3F99tRoVK5r8WL06zJ0LoaHw5ZdQpYpZ99PfP2Xyjsm89/R7dKvSzdrAARYuhJdeMq3jf/xRpm/NxOSdd4J82fPRplIbpv07jaiYKKvDEUIkobVpd9anDxQtCi+8AJcvZ2H0aDhxApYuhXbtIEuCNnLz98xnyNohdHq4Ex82+NC64G9Zt87cW3/sMXOv3dPT6oiEhSS5O0nfGn25FHWJubvnWh2KEMLm5EkYNcoMLvP00yTqk/7TT5t5++3k+6RvDd1Kt/ndqFOiDpNbTUZZPYzr33+b6VvLl4dly2RediHJ3Vka+jakfL7y/LD1B6tDESJTu9Un3c/PjO8+bJgZ12XKFFLUJ/3k1ZO0nNmSQjkKsbDTQrJ7On8Y10R27zYHU6CAmYkmv5NGvRPpmiR3J1FKEVAzgN+O/8Z/5/6zOhwhMp3t200r92LFTKv3XbvgnXfg4EEICYGePe9/wRtxI4KWM1sScSOCpV2XUti7sFNiv6ujR6FxY/DwgF9+cfLQdyI9k+TuRL2q9SKLexaCtgZZHYoQmUJYmJn8rHp1qFHDzHb63HOwapXJiyNHQtmyKdtWnI6j2/xu/Hv2X35u/zOPFHrEobHf1+nT8OyzZtD61auhXDrogifSDUnuTlTAqwDtKrVj6r9TuX7zutXhCOGSbvVJ79DBXMi+9trtPumnT5vB2ho3NssexDtr3mHRvkWMbTIWv/J+jgk+pS5cMIn9zBlzsLea7wthI8ndyfrV7MflqMvM/m+21aEI4VKSmye9f38zilxa50mftH0Sn/7xKf1r9Wdg7YH3f4EjXb5sfp0cPmya8depY208Il2yZPjZzOxpn6d5qMBD/LD1B3pWk3nehUiLiAjTOG7SJPjtN9Otu2lTUxXfokXirmupFXI0hH5L+9G4bGO+8vvK2pbxERHQvLlpMLBwIcj0qeIu5MrdyZRSBNQI4M+Tf7Lz7E6rwxEiw0nYJ71IEfPvuXPE90lftgzatrVPYj948SDtZrejfL7y/Nz+ZzzcLLweioyEli1h0yaYOdNUTwhxF5LcLdCzWk+yumdl/Jbx919ZCAHAqVPw8cdQseLtPuldutyeJ/1ufdJTK/xmOC1mtEChWNp1KXmy5bHfxh/U9eumH/uGDTBtmhlRR4h7kORugXzZ89HpkU5M/XcqV6OvWh2OEOlWwj7ppUqZe+rFisFPP5m2ZBMmOGae9JuxN/lg9wccvnSYBZ0WUCZvGfvu4EFER5uqiLVrYfJk6NrVulhEhiHJ3SIDHxtIxI0Ipv4z1epQhEh3duy4s0/60KG3+6T36OG4Qdi01ryy4hW2Xd7GhJYTeMrnKcfsKCVu3DDN/leuNL9kevSwLhaRodw3uSulOjsjkMzmseKPUbt4bb7d/C1aa6vDEcJyYWHwzTemT3r16hAUlLhP+kcfpbxPelp8/dfX/LD1B7qW7Gpto9foaGjfHpYsge+/NwPeC5FCKblyn6qUWqeUquzwaDKZgY8NZF/YPtYeWWt1KEJY4l7zpKelT3pqLT+wnMGrB9PmoTa8UNrCZJo0sb/0knWxiAwpJcm9JuAJbFdKfa6U8nZwTJlGx4c7UtCrIN9u/tbqUIRwqlt90u82T/qAAZA3r3Nj2nVuF53ndqZakWpMazMNN2XRXcvoaNNgbulSSewi1e579mqtd2qtnwICgG7APqVUF4dHlglk9chKQM0AluxfwtHLR60ORwiHujVP+lNPmVnYbs2TPm/enfOkO9u5a+doMaMF3lm8Wdx5MTmyWDSrWlQUtGlj+vONHy+JXaRain+aaq1/AioCC4FpSqlflVIPOyqwzKJfzX4olHSLEy4puXnSz50zXdpuzZNurz7pqRUVE0XrWa05d+0cS7osoXiu4tYEcu2a6ce+YoUZBL9fP2viEC7hgeqdtNZXtNYDgMeAApiq+i+UUjkdEl0mUDJ3SVo/1JqJ2ybKePPCZSQ3T3qnTrf7pA8Zkj4mMNNa88LiF/jz5J9MazONmsVqWhNIePjt+xNTpkBAgDVxCJeRouSulPJUStVWSr2qlJoBzAMexgxfOwDYq5R63oFxurSBtQcSdj2MmbtmWh2KEKl2v3nSf/zRMX3S0yJwYyAzds4g8JlA2lW2aGCYK1egSRPzy2f6dDP3rBBplJKucH8AV4E/gS+ACsASoBNQAigEzALmKqXkBlEq1PepT5XCVRi7aax0ixMZTnLzpCfsk56SedKtMPu/2bz363t0r9Kdd558x5ogzp+HZ54xrQhnzzZD7glhBym5co8APgYaA3m01rW01q9predorUO11le11m8A7wJDHRlsRhYcHIyvry9ubm74+voSHBwc/5xSitfrvM7OczulW5zIEJKbJ/3ZZ53fJz21Np/aTM+FPalXsh4TWk6wZjKYU6egfn3YvRsWLDCND4Swk5S0lm+stR6htV6rtb52j1U3YK7kRRLBwcEEBARw7NgxtNYcO3aMgICARAm+yyNdKJyjMGP+HGNhpELcXdI+6UnnSf/5Z+f2SU+t41eO8/zM5ynqXZQFnRaQ1SOr84M4fNh0Gzhxwow+17y582MQLk3ZqxpYKZUdeFZrvcQuG3xAtWrV0lu2bLHb9kJCQmiQhukU3+7fH99Vq5hTogR//PUX0dHRd6yTNWtW6iSYi/mYzzGOljlKrb9qkSMyHdZj2sHly5fJkyeP1WG4FEeXaWRkcc6c8ePs2SbcuFEQD48rFC78C0WKrMDb+5DD9usIMe4x7Kixg6hsUVTfWv2enzNHlavPtWt8/u+/ZImL4+1HH2Vvrlx230d6lZE//9WqVWPs2LFWh3EHpdRWrXWtpMvtNn+h1vo65l68AKrv30/nw4epGB5Oo2QSO3BHwi8aWpTjPsc5VeIUFfZXcEaYQiQrNjY758/X5/TpZly9WgWIJV++zRQp8i358/+Om1uM1SE+MI1mb+W9XPO6xqM7H7XkB3Tlq1cZvXMnN93cGFStGkfSY2ME4Rq01i7xV7NmTW1Pv/76a9o38umnWoP2yZZNA3f8+fj43PGSvov76mwjs+lzEefSvv90yC7lKhKxV5nGxWm9caPWvXtrnSOH1qB1hQpajx6t9alTdtmFpQavHKz5AD1u87gUrW/3c3XFCq29vLQuW1brw4ftu+0MQj7/9gds0cnkRJkVzpHeegsmTyYwOhovt8RF7eXlRWBg4B0vGVRnEFExUTKojXCaW33SK1Y0t4HnzIHOnR03T7oVgrYGMWbTGF6t/Sr9H+vv/ABmzjQD1FSoYAq2dGnnxyAyFUnujtarF/6LFhHk5oaPhwdKKXx8fAgKCsLf3/+O1SsXrIxfOT/G/T2O6Jjkq/OFSKvk+qQXLWqmCz9zBiZOTH990lNrzeE19F/WH79yfnzR5AvnB/Dll2YO9nr1TN/AwoWdH4PIdCS5O0PLlviHhHA0Z07iChXi6Pz5ySb2WwbXHczZa2eZ9u80JwYpMoN7zZO+fj306pU++6Sn1p7ze2g/uz2VClZiVvtZeLjZrZnR/cXFwRtvwODBZoa3lSshd27n7V9kapLcnaVePVMdlzWr6du6atVdV21UuhE1itbgsz8+IzYu1olBCleUsE+6lfOkO9v5a+dpMbMF2TyysbTLUnJldWKr9Ohoc7U+Zgy88oqZuzZbNuftX2R6ktydqVIl+PNP803avLmpA02GUooh9YawP2w/C/cudG6MwiXcq096aKjz50l3tqiYKNr83IbQ8FAWdV6ETx4f5+380iUznOzPP8Onn8JXX7luQYt0y4l1VAIw37QbNphquj59zGAWI0bccXOzbaW2lMtXjtG/j6ZtpbbWjKAlMpyDB81vxp9+MgOg5c8PL79sTjWrplN1Nm2bDOb3E7/zc/ufebzE487b+eHDZgKYI0cgONhcvQthAblyt0KuXGa+5hdegJEjoVs3U42XgLubO2898RZbQrew7sg6iwIVGcH16+5MnmxmXytfHkaPhmrVbs+TPnZs5knsACPWj4ifDKbjwx2dt+NNm6BOHTOn7S+/SGIXlpLkbhVPT5gwwfRBmjHDDMwdFpZolR5Ve1DEuwijfx9tUZAivUo4T3rbtk/Qp4/JKaNHp5950q0Q/G8wH6z/gJ5Vezp3Mpi5c6FhQ8iZ09x6e/pp5+1biGRIcreSUvDOO+YG6N9/w+OPw5498U9n88jG63VeZ83hNWwN3WphoCK9SG6e9GeeOcfvv5tTxxX6pKfWxmMb6bO4D/V96hPUMsg5t7K0NrVvHTqYGXQ2bTIDBghhMUnu6UGnTvDrrxAebqr1ErSkf6nWS+TOmluu3jOx+82T/tZb+1ymT3pqHQg7QJuf2+Cbx5f5neaTxd0JVRbXr4O/P7z3HnTvDmvXQsGCjt+vECkgyT29qFvXXL2XLm0a5Hz1FWhNrqy56P9Yf+btnsee83vuvx3hMjLqPOnOFhYZRvMZzVFKsbzrcvJlz+f4nZ45Y6rhZ86Ejz82LRilq5tIRyS5pyelSsFvv5lhKgcNghdfhOhoXq/zOtk9sxO48c7haoVrSTpP+q0+6StXunaf9NSKjommzc9tOH7lOAs7LaRsPicUzubNULMm7NwJ8+fDkCGZu9pEpEuS3NMbb2/zhfHuuzBpEjRoQMGrMQx4bAAzd81kf9h+qyMUdpaSPulNmkhX6aRudXnbeHwjU1pPoV6peo7f6dSpprFDliym4VybNo7fpxCpIMk9PXJzM5doc+aYq4OaNXkjSwOyumeVq3cXcuCAqWb38TF3Ytatg/794Z9/YMsWGDAA8jmhhjmjGh4ynOCdwQQ+E0jnRzo7dmcxMWYY2Z49zaD7f/+dufoXigxHknt61r69uTrIlo3CjdvwctZ6BP8bzMGLB62OTKRSRIQZZOapp0yL908+gapVze+40FAzx4jkjPubvH0yH234iBeqv+D4Lm9nz5quql9+aYaSXbUKChRw7D6FSCNJ7undo4+aq4SGDXnzwzV4xsKokBFWRyUeQMI+6UWKcEef9GXLzO+4zNYnPbXWHl5LwNIAnivzHN83/96xXd42bTL31//6yzSa+/prM0aFEOmcJPeMIH9+WLaMoq+/R8DmWKb+M40j/26wOipxHwnnSb/VJ93V5kl3tl3ndtFudjseKvAQczrMwdPdQYlWa4otWpT4/nqPHo7ZlxAOIMk9o3B3hxEjeLvvFDziIPCj52DxYqujEkncbZ70KVPg9GnXmifd2ULDQ2kW3AwvTy+WdV1G7mwOmj41PBz8/akwdqypjt+yxYznK0QGIsk9gynWticBD/dgSuUbHOjdCt56C27etDqsTC+5edLfecc0mlu/3rTD8va2OsqMKzw6nOYzmnMp6hLLui6jVO5SjtnRzp1Qqxb8/DOHX3zRjOMrrRpFBiTJPQMa2uITsmb14oOA8vD559Cggbl5K5wq6TzpP/yQuE/6yJFQrpzVUWZ8N2Nv0mFOB3ae3cmcDnOoXrS6/XeiNfz4oxkC+upVWLuW4/7+pueKEBmQnLkZUBHvIrz2+GvMzHaQf6d8Yq42qlaFBQusDs3lJdcn3c3N9Ek/fVr6pNub1pr+y/qz6tAqxrcYT9NyTe2/kytXoEsXM2jUE0+YapgGDey/HyGcSJJ7BvXWE2+RK2su3sv6O2zbZoYta9vWdJS+ft3q8FzOgQPm/nnCPukvv2zywNat0ifdUUZuGMnE7RMZ9tQwXqzxov13sHmzqXaZO9e0fly1CgoXtv9+hHAySe4ZVN7seXnribdYvG8xf2ULM02w33wTvv8eatc2V/MiTZL2SU9unvSqVa2O0nVN3j6Z90Pep2fVnnzU8CP7bjw21owJX68exMXBhg2mkYRUuQgXIck9A3utzmsU9CrIsHXDTHedzz4zN3zPnzeNgsaMMV9cIsWS65N+/rzMk+5sKw+upO+SvjQu25gJLSfYty/70aNm0pehQ83wsTt2mOp4IVyIJPcMzDuLN0OfGsraI2tZe3itWdikiblq9/ODN94wLbxOnrQ20Azg1ClzIZe0T/pvv8k86c62NXQr7We3p0rhKsztMNd+fdm1hunTTXXLjh1mUJqff4Y8eeyzfSHSEUnuGdxLtV6iZK6S/G/N/4jTtqv0ggVN47oJE8wIW48+CtOmmS83Ee9Wn/RmzcyEfEOH3u6TfuaM6ZNer570SXemQxcP0WxGMwp4FWBZ12XkzJrTPhs+fx46dDDzrj/6qBnAv0cPeXOFy5LknsFl88jGqEaj2HZ6GzN3zrz9hFKm9e+OHVC5svkia9PGjJOdySWdJ33nztvzpN/qky7zpDvfuWvnaDK9CbFxsazqtoqiOYvaZ8MLFsDDD8OSJeb+yvr1ULq0fbYtRDplSXJXSjVVSu1TSh1USg25x3qPKaVilVLtnRlfRtP10a7UKFqDoeuGEhUTlfjJ8uVNY6FPPzX34x9+2FRFOvkqPjg4GF9fX5555hl8fX0JDg526v5lnvT0LeJGBM2CmxEaHsrSrkupWKBi2jcaFgbduplGEiVLmm4Nb78tjeZEpuD05K6UcgfGAX5AZaCLUqryXdb7BFjl3AgzHjflxufPfc7xK8f5+q+v71zB3d2MZLdtG5QpY24mt2ljmnw7QXBwMAEBARw7dgytNceOHSMgIMDhCV7mSc8YbsTeoN3sduw4s4PZHWZTp0SdtG1Qa5g9GypVMj9kP/jA3J565BG7xCtERmDFlXtt4KDW+rDW+gYwC2iVzHqvAPOAc84MLqNqWLohLSq0IHBjIBciLyS/UuXK8Mcf5ip+1Srz+McfHX4VP2zYMCIjIxMti4yMZNiwYQ7Z36150kuVurNPusyTnr7E6Tj6LOrD6kOrCWoZRIsKLdK2wdBQc6XeqZMZlGDbNhg+XGZyE5mO0k6unrVVsTfVWr9oe9wdeFxrPTDBOsWBGcAzwI/AUq313GS2FQAEABQuXLjmrFmz7BZnREQE3hlsMPCj147ywpYXaFW8Fa+We/We62Y/eZKKn39Onn/+4VK1aux//XWul3LMeN3PPPMMyZ1nSinWrVtnl31cv+5OSEhBVqwows6deXBz09SufZGmTU/zxBNheHq6bmPCjHiughl9btyhccw7NY8XfF+gm0+31G8sNpZiS5ZQZuJE1M2bHO3Th5Pt26PTUC2TUcs1PZMytb+GDRtu1VrXuuMJrbVT/4AOwMQEj7sD3yRZZw5Qx/b/KUD7+223Zs2a2p5+/fVXu27PWfot6ac9Rnjovef33n/l2Fitf/hB6zx5tM6SRev339f6+nW7x+Tj46OBO/58fHzStN24OK03bNC6d2+tc+TQGrSuUEHrjz/W+tQp+8SeEWTUc3XUhlGaD9CDVgzScXFxqd/Qjh1aP/64OQEaNdJ6/367xJdRyzU9kzK1P2CLTiYnWlEtfxIomeBxCSDpzd9awCyl1FGgPfCdUqq1U6LL4D5s8CFenl4MWjUo2avlRNzcICDATC7evj2MGAFVqpgqezsKDAzEy8sr0TIvLy8CAwNTtb1Tp+49T/qQIdInPb2buG0iQ9cNxf9Rf75o8kXqBqm5etWM5VCzJhw+bPqw//KLaUQqRCZnRXL/GyivlCqtlMoCdAYSTUyutS6ttfbVWvsCc4H+WuuFTo80AyrsXZgP6n/AyoMrWbp/aQpfVBiCg2H1anP/vWlTc9/y2DG7xOTv709QUBA+Pj4opfDx8SEoKAh/f/8UbyPhPOmlSplx3osVS9wnXeZJzxjm75lPv6X98Cvnx+RWk3FTD/g1dGswmooV4csvzTCCe/eCv7+cAELYOD25a61jgIGYVvB7gNla6/+UUi8ppV5ydjyuaGDtgVQqUIlBqwbd2TXuXp57zkxEfmsCjUqVTB8xO0xE4+/vz9GjR1m3bh1Hjx5NcWJPbp70oUNNo7mQEOmTntGsObyGLvO68Hjxx5nTYc6Djz63fTvUr28GoylZEv76y/RrlBaSQiRiST93rfVyrXUFrXVZrXWgbdl4rfX4ZNbtpZNpTCfuztPdk6+afsXhS4f58s8vH+zFWbOaCTT27IHmzeH99+Ghh0y/MSc1vgwLg2++STxP+rPPmt8bt/qkyzzpGc+mk5toPas1FfNXZFnXZeTI8gC/ys6ehb59TRX87t0moW/aBI895riAhcjAZIQ6F/Vc2edo81AbRm4cycmrqRhbvlQpUw/+66/mqqhLFzMW619/2T9YTJ/0lStv90l/9dXbfdJPnzbdlRs3lj7pGdWuc7toFtyMIt5FWNVtFXmz503ZC69fN6PKlS9v7sG8/roZSrBvX9NmRAiRLPl0uLAxTcYQp+N465e3Ur+RBg1M5/AffzSNlurUMWN0HzhglxgTzpPu52f6pPfvb4b+lj7pruHQxUM0ntaY7J7Z+aX7LykbVjY21iTzChVMTVL9+vDff/DFFzLRixApIMndhfnm8eXtem8za9cs1hxek/oNububRksHDpgBQVasMAPgDByYqrHq7zdP+pdfmkb7IuM7ceUEjaY24kbsDVZ3W03pvPcZ011rWL7cjBHcu7eZd/fXX8248BUqOCdoIVyAJHcXN+TJIZTPV56Xl73M9ZtpbBiXM6cZyvNWtej48WY423fegYsX7/lSrc30qZ98UjF+nvRz52SedFd2NuIsz057lktRl1jVbRUPF3r43i8ICTG/+Jo3N78AZ80yt4EaNHBGuEK4FEnuLi6bRzbGtxjPwYsHCdyYun7ldyhSBL77zjS6a90aPvnEzLL14Ydw+XKiVRPOk/7UU7B+fcH4edL37pV50l3VxesXaTy9MSeunGBZ12XULFbz7iv/+afpqdGwoWkxOX68Obc6dZL76kKkknxyMoFnSj9Dj6o9+PT3T/nv3H/223D58qZ//L//mubsH3wAvr7EDH2fhZMvJpon/Vaf9Hnz/pB50l3c1eir+AX7sffCXhZ1XsSTpZ5MfsUNG8x588QTppHFmDHm1k+/flKFI0QaSXLPJD5/7nNyZc3FS8teIk7H2XfjjzwC8+axd8Y2dhRohMfHH9Gojw/NfxvCqFfPcPDg7T7p2bPbed8iXbk1deu209uY02EOz5V9LvEKWptuEQ0amEZyu3aZRnJHjpiW8NmzWxK3EK5GknsmUTBHQT5v/Dm/Hf+NH7f9aLftJuyTXqlrdeqcnMdbTf4lon5z+kd8ytvf+1L20352a10v0q/Im5G0nNmSTSc3MbPdTJ6v+PztJ2NiYOZMc6L4+Zl2G199ZZL64MEyEpEQdibJPRPpWbUnDXwb8OYvb3LiyolUbyfpPOmvvmpujX7zjWnt/tnKRykaMgu1bx/06gU//WRuurdtS+5//3XaYDjCeaJiomg9qzXrj65napuptK/c3jxx5Yrp/lC+PHTtasYRnjTJdKt89VW5UhfCQSS5ZyJKKX58/kdi4mLou6Tv/SeWSSK5edJfeskMEbt1q+kZl6hPevnypnHUsWOmRf369VR/7TWoVQumTTNf9CLDi4qJot3sdvxy+BcmtZpE10e7wqFDMGiQGSJ28GDz74IFpq96795yT10IB5PknsmUyVuGT5/9lFWHVjFp+6T7rp+0T/onn5ia1blzzVX6V19B1ar32UjhwhAYCCdOsG/wYDPqWI8e5gv/nXdMC2mRIUXHRNNudjuWH1jOhOY/0Ot4PjPxULlypkdFq1ZmNKING0zPCmn9LoRTyCctE3r5sZdp4NuA11e9zvErx+94XmvYuNH0Rb9bn/R27VJx8eXlxemWLc3V2+rVpsn8p5+avvLNm5sruxs37HOQwuGiY6JpO7utSew8z4vtAk0y37XLdIs8etTU0NS8Rzc4IYRDSHLPhNyUG5Oen0ScjuPFxS/GV8+fPGkmhKtQ4fY86Z063Z4n3W590pUy/ZoXLDAJ4N13zWxfbdtCiRLw5pvmB4BIt6LDL9N2bF2WH1hO0BJ48cMlpl3F/PnmPX3/fRnAQAgLSXLPpErnLc1nz33GL4d/ISDoB/z8zPjuw4ZB8eK350n/8UcHz5NesiSMGAHHj5sqgaeeMnX9jzxihiD94gtT/y+sFxcHv/7K9Rd70vqVAiy/tp0ffs9H3xbDTQO51auhTRvw8LA6UiEyPUnumdT27bB7Wj88jj/HxBOD2X5yD0OHkqhPulN7J3l4mKr5efPMsHZjx5plb75pruYbNTKN81Ixlr1Ig7g4U3UzaBCUKkVkk2doGRfMKt9YJpYbTMCKc/GDFwkh0g9J7plIwj7pNWpA0A9u+EX9RO7sOSg6oCvvDo+mbFmrowQKFYLXXoPNm2HfPnjvPZPwX37ZVPU2bAhff236SAv7u3kT1qwxXdVKlYInn4Tx44moXY1mIyvxq69mSpufeMH/C5mDV4h0SpK7i0tunnQ3t9vzpC8OLsq09pPYcXYHQ9cOtTrcO1WoYBpn7dljhrl9913Tuu+110xDvCpVzL2E334zA6WI1Dl3DqZPhy5doGBB0yZiwgTTbXH6dK6cOECTFpf4LWo/09pMo0fVHlZHLIS4B7k55qIOHDBd2KZONRe9+fObedJ7975zOtWWFVsy4LEBjNk0hiblmtC4bGNrgr4XpeDRR83fhx+aA1yyBBYvNs34R42C3LlN9X3jxvDMM6Y7lgxgn7zISDNhy9q1sGoVbNtmlhcqBO3bw/PPm3Hfvbw4f+08TaY3Yde5XcxqP+v2ADVCiHRLkrsLiYgwLdwnTzZd2dzcTJfjsWPNd/W9uq599txnhBwNoceCHvzz0j8U9i7stLhTpXx5MzjK4MFw6dLtJLVqlWmxDaZlYIMGpun/E0+YOegzaz/ry5fN9Km//24aVfz1l+l26O5uyiYw0Jws1aolKqPQ8FCenfosRy4fYWHnhTQr38yqIxBCPABJ7hncrXnSJ0+G2bPh2jVTkz16NHTvnvLeSNk9szOz3UxqT6xN1/ldWd1tNe5uGeR+at685mqzfXtTIPv3mwQWEmLuHQcHm/Vy54bHH4fatU3f6xo1TGt9V7u6j4qCnTvN4DFbt5pE/t9/pmzc3Mxxv/aaabtQrx7kypXsZo5cOsKz057l3LVzrPRfSX3f+k4+ECFEaklyz6BOnTJV7pMnmxpqb2/o3NkMOFO3bury1aOFH+X75t/Te1Fv3v/1fQIb2Wn+d2dSyvS3rljRTB2qtRkK9c8/4Y8/zL8ff2waIwAUKGCq+h95xPw9/LD5dVSgQPpP+jduwJEj5P/9d3Nsu3aZpL537+32B/nzw2OPmQELnnjC/LDx9r7vpnee3UmT6U2IiolibY+11C5e28EHI4SwJ0nuGUh0tLnFPHmyqX2OizOzZg4bZi5a7dF1rVe1Xvx+/HdG/TaKOiXq0LJiy7Rv1EpKmXvv5cqZqgwww9/++6+5qt22zSTFyZPNfY1bcuc2Vf9ly5oBAEqVMn/Fi5vhdAsVAk9Px8WtNYSHm8EGzpwxIwwdO3b77+BB01sgNpZHb73Gx8f8QGnZ0tRM1Kxplj3gj5Tfj/9Oi5kt8PL0YkPvDTxS6BG7H54QwrEkuWcAO3aY3DN9Oly8aLp9v/OOmXCtXDn77++bZt+w7cw2eizswdaArZTJW8b+O7FS9uymev7xx28vi4szA+ns3m2qQvbvN//+/ffdh8XNn9/85cljbg3kyWN+YXl5mb+sWc097Vt/Wpsag5gY83f9umnYFhlpflhcumT+Ll6ECxfM88nt08fHVK136gQVKrAtPJwa3bvftXr9QSzdv5QOczpQKncpVndbjU8enzRvUwjhfJLc06mwMJgxw8yOuWOHyROtW5tq90aNHNu9OJtHNuZ2mEuNoBq0m92O3/v8jpenl+N2mB64uZmBWJIbjCUuznQVO3bM9B88c8YMpnPmTOKEfOTI7WQdGWnufd9rf7d+BOTIYf7y5jWJu1o1k8SLFjWD+xcpYmoMSpVKVD0THBzMsGHDOH78OKU++4zAwED8/f1TXQQ/7fiJFxa/QPWi1VnedTkFcxRM9baEENaS5J6OxMbCL7+YhL5okblYrFnT9Env0iXJdKoOVjpvaYLbBtNiRgt6LuzJz+1/xk1l0pbmbm63k+yDioszb2xsrKke9/Aw20vj/fzg4GACAgKIjIwE4NixYwQEBAA8cILXWjNyw0jeD3mfZ8s8y/yO88mZNWea4hNCWCuTflunLwcOmPvmPj7g52fmSX/5ZXPFvmULDBjg3MR+S7Pyzfjsuc+Yu3suH4Z86PwAXIGbm7k3ny3b7Wp6OzTUGzZsWHxivyUyMpJhw4Y90HZuxt4kYEkA74e8T4+qPVjWdZkkdiFcgFy5W+RufdK/+sq0h3rg6VQdZHDdwew+v5sRG0ZQqWAlOj/S2eqQBHD8+J1T9d5reXLCo8PpNLcTKw6u4N2n3mVEwxGo9N5DQAiRIpLcnchefdKdSSnF9y2+58DFA/Re1JsyectIt6h0oFSpUhw7dizZ5Slx4soJWsxswa5zu/ihxQ8E1Aywd4hCCAtJtbwTJDdPeufOJtHbdZ50B8ninoV5HedR1LsoLWe25NDFQ1aHlOkFBgbi5ZW4kaOXlxeBgfcfm2Dzqc3Unlibo5ePsrzrcknsQrggSe4OEh1tknjCedKLFbs9T/rEiWZwsIxSC1owR0FW+K8gNi6WxtMbczZCpl61kr+/P0FBQfj4+KCUwsfHh6CgoPs2ppvz3xzqT6lPNo9s/NHnD5qUa+KkiIUQziTJ3UF++snMxLZrF/HzpK9fb8E86XZUsUBFlnVdxpmIM/gF+3E1+qrVIWVq/v7+HD16lHXr1nH06NF7JvY4Hcf7v75Px7kdqVG0Bptf3MzDhR52YrRCCGeSe+4O0rGj6TLt6D7pzvZ4iceZ13EeLWe2pM3PbVjedTlZPbJaHZa4h6vRV+m+oDuL9y2md7XefNf8O7J5ZLM6LCGEA8mVu4PkyWNmHnWlxH5L03JNmfT8JNYdWUfneZ25GXvT6pDEXRwIO0CdiXVYtn8ZXzf9mh+f/1ESuxCZgCR3kSrdq3bnG79vWLh3oST4dGrh3oXUmlCLc9fO8Uv3X3jl8Vekq5sQmYQkd5FqA2sPZGyTsczfMx//+f7ExMVYHZIAYuJi+N8v/6PNz22okL8CWwO20rB0Q6vDEkI4kdxzF2nyWp3XiNWxvLH6DdyUG9PbTsfDTU4rq5wOP03neZ3ZcGwDL9d6mS+bfCltIoTIhORbWKTZ4LqDiY2L5X9r/kdUTBSz2s+S+7oWWHFgBT0X9iTiRgTT2kyjW5VuVockhLCIVMsLu3ir3lt84/cNi/ctpun0plyJumJ1SJnGjbgbDF41mGYzmlHEuwhbArZIYhcik5Mrd2E3A2sPJH/2/PRY2IOGPzVkhf8KCnsXtjosl7b3wl4Gbh/IgYgDDHhsAJ83/lxqTYQQcuUu7KvLo11Y0mUJey/s5cnJT7Lvwj6rQ3JJcTqOr//6muo/VOds1FkWdFrAt82+lcQuhAAkuQsHaFquKWt7rOVK1BXq/FiHXw79YnVILuX4leM8O/VZXlv5Go1KN2JSrUm0fqi11WEJIdIRSe7CIeqWrMvmvpspmaskfsF+fLv5W7TWVoeVocXpOMZvGc8j3z3C36F/M6HlBJZ0WUL+rPmtDk0Ikc7IPXfhML55fPm9z+90W9CNV1a8wr9n/6WdVzurw8qQ9l7YS98lffnt+G80Kt2ICS0nUDpvaavDEkKkU3LlLhwqZ9acLOi0gHeefIcJ2ybQf3t/uQ//AKJiohixfgRVx1flv3P/MbnVZH7p/oskdiHEPUlyFw7nptwY1WgUy7suJ+xGGDWDahL8b7DVYaVrWmuW7FvCw989zPCQ4bR5qA17BuyhV7VeMoSsEOK+JLkLp/Er78eEmhOoUbQG3RZ0o/uC7ly8ftHqsNKd/WH7aTGzBc/Pep6s7llZ030Ns9rPkm6FQogUk+QunKpg1oKs67mOD+p/wKxds6g8rjLz98y3Oqx04WzEWQYsG0DlcZXZeGwjXzT+gn9e+odGZRpZHZoQIoOR5C6czsPNg+ENhrOl7xaK5SxGu9nt6DS3E6HhoVaHZonw6HA+Wv8R5b4pxw9bf6BfzX4ceOUAg+sOxtPd0+rwhBAZkCR3YZmqRary14t/EfhMIAv3LqTCNxUYtXEUUTFRVofmFFejrzJq4yh8v/Ll/ZD3aVy2Mf/1/49xzcdJFbwQIk0kuQtLebp7MvSpoezuv5vnyj7HsHXDqDSuEnN3z3XZfvEXIi8wYv0IfMf6MmzdMOqWqMtfL/7FvI7zqFigotXhCSFcgCR3kS6UzVeWBZ0WsKb7GryzeNNhTgeq/1Cd+XvmE6fjrA7PLvac30O/Jf0o+WVJhocM58lST/J3379Z2nUptYvXtjo8IYQLkeQu0pVGZRqxvd92fmr9E5E3I2k3ux3Vxldj1q5Z3Ii9YXV4Dyw6Jpo5/82hyfQmVP6uMj/98xPdHu3Grpd3sbjLYmoVq2V1iEIIFyQj1Il0x8PNgx5Ve+D/qD+zds1i5MaRdJnXhSLeRehboy8BNQMokauE1WHeldaaHWd2MO3faUz9Zyph18MokasEIxqM4KVaL1EwR0GrQxRCuDhJ7iLdcndzx7+KP10e7cKqg6sY9/c4Rm4YyaiNo2hctjEdH+5I64dakydbHqtDRWvN9jPbmfPfHObsnsOhS4fwdPOk1UOteKH6CzxX5jnc3dytDlMIkUlIchfpnptyw6+8H37l/Thy6QhBW4OYuWsmvRf1pt/SfjQu2xi/cn40Kt2ICvkrOG0Et7MRZ1lzeA2rD6/ml0O/cDriNO7KnWdKP8Pb9d6mTaU2FPAq4JRYhBAiIUuSu1KqKfAV4A5M1FqPTvK8P/C27WEE8LLW+h/nRinSo9J5S/Pxsx8zqtEoNp/azOz/ZjNvzzyW7l8KQPGcxWlYuiE1itSgetHqVC1clbzZ86Zpn1przl47y74L+9h2eht/nfqLzac2c+TyEQDyZ8/Ps2WepUnZJrSs2FISuhDCck5P7kopd2Ac8BxwEvhbKbVYa707wWpHgPpa60tKKT8gCHjc2bGK9EspxeMlHufxEo/zeePPOXTpEGsPr2XtkbWsObyG6f9Oj1+3qHdRfPL4UCp3KUrlKkV+r/x4Z/HGO4s3OTxzEKfjuBF7g5txN7l+8zrnI89z7to5zl47y/Erx9kftp+r0Vfjt1cyV0keL/E4/R/rT0PfhlQvWh03JW1ThRDphxVX7rWBg1rrwwBKqVlAKyA+uWut/0iw/iYg/baeEpZTSlEuXznK5StHv1r9AFNlvuPMDraf2c7+sP0cv3Kc7ae3s2jvIqJjo++9PRT5vfJTKEchiuUsRvcq3amQvwIV81ekSuEqFM1Z1BmHJYQQqWZFci8OnEjw+CT3vip/AVjh0IiEyynsXZgm5ZrQpFyTRMu11kTFRBFxI4KIGxFcu3kNN+VGFvcseLp5ks0jG/m98uPhJs1RhBAZlxXfYMm1dkp2KDKlVENMcn/yLs8HAAEAhQsXJiQkxE4hQkREhF23JwwpV/uTMnUMKVf7kzJ1HiuS+0mgZILHJYA7ZgxRSlUBJgJ+Wuuw5DaktQ7C3I+nVq1aukGDBnYLMiQkBHtuTxhSrvYnZeoYUq72J2XqPFa0AvobKK+UKq2UygJ0BhYnXEEpVQqYD3TXWu+3IEYhhBAiw3L6lbvWOkYpNRBYhekKN0lr/Z9S6iXb8+OB94H8wHe2PssxWmsZp1MIIYRIAUtaDWmtlwPLkywbn+D/LwIvOjsuIYQQwhVI51whhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhBDCxUhyF0IIIVyMJHchhMMFBwfj6+uLm5sbvr6+BAcHWx2SEC7Nw+oAhBCuLTg4mICAACIjIwE4duwYAQEBAPj7+1sZmhAuS67chRAONWzYsPjEfktkZCTDhg2zKCIhXJ8kdyGEQx0/fvyBlgsh0k6SuxDCoUqVKvVAy4UQaSfJXQjhUIGBgXh5eSVa5uXlRWBgoEURCeH6JLkLIRzK39+foKAgfHx8UErh4+NDUFCQNKYTwoGktbwQwuH8/f0lmQvhRHLlLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi5HkLoQQQrgYSe5CCCGEi1Faa6tjsAul1HngmB03WQC4YMftCUPK1f6kTB1DytX+pEztz0drXTDpQpdJ7vamlNqita5ldRyuRsrV/qRMHUPK1f6kTJ1HquWFEEIIFyPJXQghhHAxktzvLsjqAFyUlKv9SZk6hpSr/UmZOonccxdCCCFcjFy5CyGEEC5GknsylFJNlVL7lFIHlVJDrI4no1NKlVRK/aqU2qOU+k8p9ZrVMbkSpZS7Umq7Umqp1bG4AqVUHqXUXKXUXts5W9fqmFyBUup12+d/l1JqplIqm9UxuTJJ7kkopdyBcYAfUBnoopSqbG1UGV4M8IbWuhJQBxggZWpXrwF7rA7ChXwFrNRaPwRURco2zZRSxYFXgVpa60cAd6CztVG5Nknud6oNHNRaH9Za3wBmAa0sjilD01qf1lpvs/0/HPNlWdzaqFyDUqoE0ByYaHUsrkAplQt4GvgRQGt9Q2t92dKgXIcHkF0p5QF4AaEWx+PSJLnfqThwIsHjk0gishullC9QHfjL4lBcxVjgf0CcxXG4ijLAeWCy7VbHRKVUDquDyui01qeAz4HjwGngitZ6tbVRuTZJ7ndSySyTLgV2oJTyBuYBg7TWV62OJ6NTSrUAzmmtt1odiwvxAGoA32utqwPXAGl3k0ZKqbyYGtDSQDEgh1Kqm7VRuTZJ7nc6CZRM8LgEUn2UZkopT0xiD9Zaz7c6HhdRD3heKXUUc/voGaXUdGtDyvBOAie11rdqluZikr1Im2eBI1rr81rrm8B84AmLY3Jpktzv9DdQXilVWimVBdPoY7HFMWVoSimFuYe5R2s9xup4XIXW+h2tdQmttS/mPF2ntZaroTTQWp8BTiilKtoWNQJ2WxiSqzgO1FFKedm+DxohDRUdysPqANIbrXWMUmogsArTonOS1vo/i8PK6OoB3YGdSqkdtmVDtdbLrQtJiLt6BQi2/bg/DPS2OJ4MT2v9l1JqLrAN03tmOzJanUPJCHVCCCGEi5FqeSGEEMLFSHIXQgghXIwkdyGEEMLFSHIXQgghXIwkdyGEEMLFSHIXQgghXIwkdyGEEMLFSHIXQgghXIwkdyHEA1FKlVNK3VRKfZhk+fdKqXClVC2rYhNCGJLchRAPRGt9EDN//OtKqQIASqn3gT5AG631FivjE0LI8LNCiFRQShUBDgHfAXsx44R30VrPtjQwIQQgE8cIIVJBa31GKTUWeAPzPfKqJHYh0g+plhdCpNYBICvwp9Z6nNXBCCFuk+QuhHhgSqlngB+AP4F6SqmqFockhEhAkrsQ4oEopWoACzGN6hoAx4FRFoYkhEhCkrsQIsWUUuWAFcBq4BWt9Q3gQ6CZUuppS4MTQsST1vJCiBSxtZD/A3Ol3kRrHW1b7g7sAi5prZ+wMEQhhI0kdyGEEMLFSLW8EEII4WIkuQshhBAuRpK7EEII4WIkuQshhBAuRpK7EEII4WIkuQshhBAuRpK7EEII4WIkuQshhBAuRpK7EEII4WL+DyRcTP0gwzWYAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# and plot again\n",
"def plot_raw_data(xi, yi, ax):\n",
" \"\"\"plot x vs y on axes ax, \n",
" add axes labels and turn on grid\n",
" \"\"\"\n",
" ax.plot(xi, yi, 'ko', label='raw data')\n",
" ax.set_xlabel('$x$', fontsize=16)\n",
" ax.set_ylabel('$y$', fontsize=16)\n",
" ax.grid(True)\n",
" \n",
" \n",
"# set up figure\n",
"fig = plt.figure(figsize=(8, 6))\n",
"ax1 = fig.add_subplot(111)\n",
"ax1.margins(0.1)\n",
"\n",
"x = np.linspace(0.4, 9.1, 100)\n",
"\n",
"ax1.plot(x, p0(x), 'k', label='Constant')\n",
"ax1.plot(x, p1(x), 'b', label='Linear')\n",
"ax1.plot(x, p2(x), 'r', label='Quadratic')\n",
"ax1.plot(x, p3(x), 'g', label='Cubic')\n",
"\n",
"# Overlay raw data\n",
"plot_raw_data(xi, yi, ax1)\n",
"\n",
"ax1.legend(loc='best', fontsize = 12)\n",
"ax1.set_title('Polynomial approximations of differing degree', fontsize=16);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.4600000000000001\n",
"0.33298899237933954\n",
"0.1994782421425496\n",
"0.15730343662323767\n"
]
}
],
"source": [
"# we use the square of the difference to ensure each contribution\n",
"# to the total error is positive, otherwise errors of different signs\n",
"# could/would cancel out giving a false estimate of how good our approximation is\n",
"\n",
"\n",
"def sqr_error(p, xi, yi):\n",
" \"\"\"\"function to evaluate the sum of square of errors\"\"\"\n",
" # first compute the square of the differences\n",
" diff2 = (p(xi)-yi)**2\n",
" # and return their sum\n",
" return diff2.sum()\n",
"\n",
"\n",
"print(sqr_error(p0, xi, yi))\n",
"print(sqr_error(p1, xi, yi))\n",
"print(sqr_error(p2, xi, yi))\n",
"print(sqr_error(p3, xi, yi))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.46\n",
"0.3329889923793396\n",
"0.1994782421425496\n",
"0.15730343662323765\n"
]
}
],
"source": [
"# Alternatively we could just use the l2 norm\n",
"# [you can get different norms if you use pass arguments to the function, \n",
"# but the 2-norm is the default]\n",
"# To match the definition of the error we implemented above we need to take the square\n",
"print(sl.norm(p0(xi) - yi)**2)\n",
"print(sl.norm(p1(xi) - yi)**2)\n",
"print(sl.norm(p2(xi) - yi)**2)\n",
"print(sl.norm(p3(xi) - yi)**2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework - Degree of approximation (from Numerical Methods)\n",
"\n",
"Extend the previous question above by fitting and plotting polynomials of increasing degree past cubic. At what *degree* does the resulting polynomial approximation fit the data exactly?\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of points to fit: 6\n",
"poly_coeffs = \n",
"[[ 0.5 0. 0. 0. 0. 0. ]\n",
" [ 0.0508044 0.26714649 0. 0. 0. 0. ]\n",
" [ 0.02013603 -0.13983999 0.55279339 0. 0. 0. ]\n",
" [-0.00552147 0.09889271 -0.43193108 0.75909819 0. 0. ]\n",
" [-0.00420655 0.07403681 -0.38492428 0.59251888 0.27906056 0. ]\n",
" [-0.00301599 0.06536037 -0.49614427 1.59623195 -2.08266478 1.20030166]]\n",
"square of the difference between the data and the polynomial of degree 0 = 4.60000000e-01\n",
"square of the difference between the data and the polynomial of degree 1 = 3.32988992e-01\n",
"square of the difference between the data and the polynomial of degree 2 = 1.99478242e-01\n",
"square of the difference between the data and the polynomial of degree 3 = 1.57303437e-01\n",
"square of the difference between the data and the polynomial of degree 4 = 4.69232378e-02\n",
"square of the difference between the data and the polynomial of degree 5 = 7.23216068e-26\n"
]
}
],
"source": [
"print(\"Number of points to fit: \", np.size(xi))\n",
"# in this example we have 6 pieces of information, to fit a polynomial\n",
"# that exactly goes through these points we need 6 unknowns, or free \n",
"# parameters, to choose in the polynomial. [Too few and we won't be able\n",
"# to fit the data exaclty, and too many would just be a waste. Cf. over-\n",
"# and under-determined systems.]\n",
"# A 5th order polynomial has 6 free parameters (all the powers up to 5,\n",
"# including 0).\n",
"# So calling polyfit to fit a polynomial of degree 5 (size(x)-1) should\n",
"# fit the data exactly (repeat below with size(x)-2 etc, and size(x) and \n",
"# above to convince yourself of this).\n",
"# Let's check the errors up to this point:\n",
"\n",
"xi = np.array([0.5, 2.0, 4.0, 5.0, 7.0, 9.0])\n",
"yi = np.array([0.5, 0.4, 0.3, 0.1, 0.9, 0.8])\n",
"\n",
"# Let's set up some space to store all the polynomial coefficients\n",
"# there are some redundancies here, and we have assumed we will only \n",
"# consider polynomials up to degree N\n",
"N = 6\n",
"poly_coeffs = np.zeros((N, N))\n",
"\n",
"for i in range(N):\n",
" poly_coeffs[i, :(i+1)] = np.polyfit(xi, yi, i)\n",
"\n",
"print('poly_coeffs = \\n{}'.format(poly_coeffs))\n",
"\n",
"# and now compute the errors\n",
"for i in range(N):\n",
" p = np.poly1d(poly_coeffs[i, :(i+1)])\n",
" print('square of the difference between the data and the '\n",
" 'polynomial of degree {0:1d} = {1:.8e}'.format(i, sqr_error(p, xi, yi)))\n",
" \n"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.13"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": true
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}