{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget\n",
"from mayavi import mlab"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"import scipy as sp\n",
"import scipy.linalg\n",
"import sympy as sy\n",
"sy.init_printing() \n",
"np.set_printoptions(precision=3)\n",
"np.set_printoptions(suppress=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mlab.init_notebook(backend='x3d')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start from plotting basics in Python environment, in the meanwhile refresh the system of linear equations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Visualisation of A System of Two Linear Equations "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a linear system of two equations:\n",
"\\begin{align}\n",
"x+y&=6\\\\\n",
"x-y&=-4\n",
"\\end{align}\n",
"Easy to solve: $(x, y)^T = (1, 5)^T$. Let's plot the linear system."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6a6cda854af743ebbc5838943dbf4d16",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-5, 5, 100)\n",
"y1 = -x + 6\n",
"y2 = x + 4\n",
"\n",
"fig, ax = plt.subplots(figsize = (12, 7))\n",
"ax.scatter(1, 5, s = 200, zorder=5, color = 'r', alpha = .8) \n",
"\n",
"ax.plot(x, y1, lw =3, label = '$x+y=6$')\n",
"ax.plot(x, y2, lw =3, label = '$x-y=-4$')\n",
"ax.plot([1, 1], [0, 5], ls = '--', color = 'b', alpha = .5)\n",
"ax.plot([-5, 1], [5, 5], ls = '--', color = 'b', alpha = .5)\n",
"ax.set_xlim([-5, 5])\n",
"ax.set_ylim([0, 12])\n",
"\n",
"ax.legend()\n",
"s = '$(1,5)$'\n",
"ax.text(1, 5.5, s, fontsize = 20)\n",
"ax.set_title('Solution of $x+y=6$, $x-y=-4$', size = 22)\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to Draw a Plane "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before drawing a plane, let's refresh the logic of Matplotlib 3D plotting. This should be familiar to you if you are a MATLAB user. \n",
"\n",
"First, create meshgrids."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"x, y = [-1, 0, 1], [-1, 0, 1]\n",
"X, Y = np.meshgrid(x, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mathematically, meshgrids are the coordinates of Cartesian product. To illustrate, we can plot all the coordinates of these meshgrids"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "375bc50022f44bc89ee0463add9f8606",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize = (12, 7))\n",
"ax.scatter(X, Y, s = 200, color = 'red')\n",
"ax.axis([-2, 3.01, -2.01, 2])\n",
"ax.spines['left'].set_position('zero') # alternative position is 'center'\n",
"ax.spines['right'].set_color('none')\n",
"ax.spines['bottom'].set_position('zero')\n",
"ax.spines['top'].set_color('none')\n",
"ax.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try a more complicated meshgrid."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9d8b829bf13249ce8e14e2bc73b46c82",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x, y = np.arange(-3, 4, 1), np.arange(-3, 4, 1)\n",
"X, Y = np.meshgrid(x, y)\n",
"\n",
"fig, ax = plt.subplots(figsize = (12, 12))\n",
"ax.scatter(X, Y, s = 200, color = 'red', zorder = 3)\n",
"ax.axis([-5, 5, -5, 5])\n",
"\n",
"ax.spines['left'].set_position('zero') # alternative position is 'center'\n",
"ax.spines['right'].set_color('none')\n",
"ax.spines['bottom'].set_position('zero')\n",
"ax.spines['top'].set_color('none')\n",
"ax.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now consider the function $z = f(x, y)$, $z$ is in the $3rd$ dimension. Though Matplotlib is not meant for delicate plotting of 3D graphics, basic 3D plotting is still acceptable. \n",
"\n",
"For example, we define a simple plane as \n",
"$$z= x + y$$\n",
"Then plot $z$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9a42f621a1404ea5965c21db201dc2b3",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Z = X + Y\n",
"fig = plt.figure(figsize = (9,9))\n",
"ax = fig.add_subplot(111, projection = '3d')\n",
"ax.scatter(X, Y, Z, s = 100, label = '$z=x+y$')\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or we can plot it as a surface, Matplotlib will automatically interpolate values among the Cartesian coordinates such that the graph will look like a surface."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "dd024b438d1e4be4bc39f52632f4f071",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize = (9, 9))\n",
"ax = fig.add_subplot(111, projection = '3d')\n",
"ax.plot_surface(X, Y, Z, cmap ='viridis') # MATLAB default color map\n",
"ax.set_xlabel('x-axis')\n",
"ax.set_ylabel('y-axis')\n",
"ax.set_zlabel('z-axis')\n",
"ax.set_title('$z=x+y$', size = 18)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Visualisation of A System of Three Linear Equations "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have reviewed on plotting planes, now we are ready to plot several planes all together."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider this system of linear equations\n",
"\\begin{align}\n",
"x_1- 2x_2+x_3&=0\\\\\n",
"2x_2-8x_3&=8\\\\\n",
"-4x_1+5x_2+9x_3&=-9\n",
"\\end{align}\n",
"And solution is $(x_1, x_2, x_3)^T = (29, 16, 3)^T$. Let's reproduce the system visually."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "ab9ad0b1a8f84be4a9feae9f9a22039c",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x1 = np.linspace(25, 35, 20)\n",
"x2 = np.linspace(10, 20, 20)\n",
"X1, X2 = np.meshgrid(x1, x2)\n",
"\n",
"fig = plt.figure(figsize = (9, 9))\n",
"ax = fig.add_subplot(111, projection = '3d')\n",
"\n",
"X3 = 2*X2 - X1\n",
"ax.plot_surface(X1, X2, X3, cmap ='viridis', alpha = 1) \n",
"\n",
"X3 = .25*X2 - 1\n",
"ax.plot_surface(X1, X2, X3, cmap ='summer', alpha = 1)\n",
"\n",
"X3 = -5/9*X2 + 4/9*X1 - 1\n",
"ax.plot_surface(X1, X2, X3, cmap ='spring', alpha = 1)\n",
"\n",
"ax.scatter(29, 16, 3, s = 200, color = 'black')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are certain there is a solution, however the graph does not show the intersection of planes. The problem originates from Matplotlib's rendering algorithm, which is not designed for drawing genuine 3D graphics. It merely projects 3D objects onto 2D dimension to imitate 3D features.\n",
"\n",
"Mayavi is much professional in rendering 3D graphics, we give an example here. If not installed, run ```conda install -c anaconda mayavi```."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mlab.clf()\n",
"X1, X2 = np.mgrid[-10:10:21*1j, -5:10:21*1j]\n",
"X3 = 6 - X1 - X2\n",
"mlab.mesh(X1, X2, X3,colormap=\"spring\")\n",
"X3 = 3 - 2*X1 + X2\n",
"mlab.mesh(X1, X2, X3,colormap=\"winter\")\n",
"X3 = 3*X1 + 2*X2 -4\n",
"mlab.mesh(X1, X2, X3,colormap=\"summer\")\n",
"mlab.axes()\n",
"mlab.outline()\n",
"mlab.points3d(1, 2, 3, color = (.8, 0.2, .2), )\n",
"mlab.title('A System of Linear Equations')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualisation of An Inconsistent System "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's visualise the linear system that does not have a solution."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"x+y+z&=1\\\\\n",
"x-y-2z&=2\\\\\n",
"2x-z&=1\n",
"\\end{align}\n",
"Rearrange the system to solve for $z$:\n",
"\n",
"\\begin{align}\n",
"z&=1-x-y\\\\\n",
"z&=\\frac{x}{2}-\\frac{y}{2}+1\\\\\n",
"z&=2x-1\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mlab.clf()\n",
"X, Y = np.mgrid[-5:5:21*1j, -5:5:21*1j]\n",
"Z = 1 - X - Y\n",
"mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
"\n",
"Z = X/2 - Y/2 + 1\n",
"mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
"\n",
"Z = 2*X - 1\n",
"mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
"mlab.axes()\n",
"mlab.outline()\n",
"mlab.title('A Inconsistent System of Linear Equations')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualisation of A System With Infinite Numbers of Solutions "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our system of equations is given\n",
"\n",
"\\begin{align}\n",
"y-z=&4\\\\\n",
"2x+y+2z=&4\\\\\n",
"2x+2y+z=&8\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Rearrange to solve for $z$\n",
"\n",
"\\begin{align}\n",
"z=&y-4\\\\\n",
"z=&2-x-\\frac{y}{2}\\\\\n",
"z=&8-2x-2y\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mlab.clf()\n",
"X, Y = np.mgrid[-2:2:21*1j, 2:6:21*1j]\n",
"Z = Y - 4\n",
"mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
"\n",
"Z = 2 - X - Y/2\n",
"mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
"\n",
"Z = 8 - 2*X - 2*Y\n",
"mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
"mlab.axes()\n",
"mlab.outline()\n",
"mlab.title('A System of Linear Equations With Infinite Number of Solutions')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution of the system is $(x,y,z)=(-3z/2,z+4,z)^T$, where $z$ is a **free variable**. \n",
"\n",
"The solution is an infinite line in $\\mathbb{R}^3$, to visualise the solution requires setting a range of $x$ and $y$, for instance we can set\n",
"\n",
"\\begin{align}\n",
"-2 \\leq x \\leq 2\\\\\n",
"2 \\leq y \\leq 6\n",
"\\end{align}\n",
"\n",
"which means\n",
"\n",
"\\begin{align}\n",
"-2\\leq -\\frac32z\\leq 2\\\\\n",
"2\\leq z+4 \\leq 6\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can pick one inequality to set the range of $z$, e.g. second inequality: $-2 \\leq z \\leq 2$. \n",
"\n",
"Then plot the planes and the solutions together."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mlab.clf()\n",
"X, Y = np.mgrid[-2:2:21*1j, 2:6:21*1j]\n",
"Z = Y - 4\n",
"mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
"\n",
"Z = 2 - X - Y/2\n",
"mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
"\n",
"Z = 8 - 2*X - 2*Y\n",
"mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
"\n",
"ZL = np.linspace(-2, 2, 20) # ZL means Z for line, we have chosen the range [-2, 2]\n",
"X = -3*ZL/2\n",
"Y = ZL + 4\n",
"\n",
"mlab.plot3d(X, Y, ZL)\n",
"\n",
"mlab.axes()\n",
"mlab.outline()\n",
"mlab.title('A System of Linear Equations With Infinite Number of Solutions')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Reduced Row Echelon Form "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For easy demonstration, we will be using SymPy frequently in lectures. SymPy is a very power symbolic computation library, we will see its basic features as the lectures move forward."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define a SymPy matrix:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}5 & 0 & 11 & 3\\\\7 & 23 & -3 & 7\\\\12 & 11 & 3 & -4\\end{matrix}\\right]$"
],
"text/plain": [
"⎡5 0 11 3 ⎤\n",
"⎢ ⎥\n",
"⎢7 23 -3 7 ⎥\n",
"⎢ ⎥\n",
"⎣12 11 3 -4⎦"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = sy.Matrix([[5, 0, 11, 3], [7, 23, -3, 7], [12, 11, 3, -4]]); M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Think of it as an **augmented matrix** which combines coefficients of linear system. With row operations, we can solve the system quickly. Let's turn it into a **row reduced echelon form**."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & - \\frac{2165}{1679}\\\\0 & 1 & 0 & \\frac{1358}{1679}\\\\0 & 0 & 1 & \\frac{1442}{1679}\\end{matrix}\\right], \\ \\left( 0, \\ 1, \\ 2\\right)\\right)$"
],
"text/plain": [
"⎛⎡ -2165 ⎤ ⎞\n",
"⎜⎢1 0 0 ──────⎥ ⎟\n",
"⎜⎢ 1679 ⎥ ⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢ 1358 ⎥ ⎟\n",
"⎜⎢0 1 0 ──── ⎥, (0, 1, 2)⎟\n",
"⎜⎢ 1679 ⎥ ⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢ 1442 ⎥ ⎟\n",
"⎜⎢0 0 1 ──── ⎥ ⎟\n",
"⎝⎣ 1679 ⎦ ⎠"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M_rref = M.rref(); M_rref # .rref() is the SymPy method for row reduced echelon form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Take out the first element in the big parentheses, i.e. the rref matrix."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1, 0, 0, -2165/1679],\n",
" [0, 1, 0, 1358/1679],\n",
" [0, 0, 1, 1442/1679]], dtype=object)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M_rref = np.array(M_rref[0]);M_rref"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you don't like fractions, convert it into float type."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1. , 0. , 0. , -1.289],\n",
" [ 0. , 1. , 0. , 0.809],\n",
" [ 0. , 0. , 1. , 0.859]])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M_rref.astype(float)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last column of the rref matrix is the solution of the system."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: rref and Visualisation "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's use ```.rref()``` method to compute a solution of a system then visualise it. Consider the system:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"3x+6y+2z&=-13\\\\\n",
"x+2y+z&=-5\\\\\n",
"-5x-10y-2z&=19\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Extract the augmented matrix into a SymPy matrix:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}3 & 6 & 2 & -13\\\\1 & 2 & 1 & -5\\\\-5 & -10 & -2 & 19\\end{matrix}\\right]$"
],
"text/plain": [
"⎡3 6 2 -13⎤\n",
"⎢ ⎥\n",
"⎢1 2 1 -5 ⎥\n",
"⎢ ⎥\n",
"⎣-5 -10 -2 19 ⎦"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = sy.Matrix([[3, 6, 2, -13], [1, 2, 1, -5], [-5, -10, -2, 19]]);A"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}1 & 2 & 0 & -3\\\\0 & 0 & 1 & -2\\\\0 & 0 & 0 & 0\\end{matrix}\\right], \\ \\left( 0, \\ 2\\right)\\right)$"
],
"text/plain": [
"⎛⎡1 2 0 -3⎤ ⎞\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢0 0 1 -2⎥, (0, 2)⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎝⎣0 0 0 0 ⎦ ⎠"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A_rref = A.rref(); A_rref"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In case you are wondering what's $(0, 2)$: they are the column number of pivot columns, in the augmented matrix above the pivot columns resides on the $0$th and $2$nd column."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because it's not a rank matrix, therefore solutions is in general form\n",
"\\begin{align}\n",
"x + 2y & = -3\\\\\n",
"z &= -2\\\\\n",
"y &= free\n",
"\\end{align}\n",
"Let's pick 3 different values of $y$, for instance $(3, 5, 7)$, to calculate $3$ special solutions:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ -9, 3, -2],\n",
" [-13, 5, -2],\n",
" [-17, 7, -2]])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"point1 = (-2*3-3, 3, -2)\n",
"point2 = (-2*5-3, 5, -2)\n",
"point3 = (-2*7-3, 7, -2)\n",
"special_solution = np.array([point1, point2, point3]); special_solution # each row is a special solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can visualise the general solution, and the 3 specific solutions altogether."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"y = np.linspace(2, 8, 20) # y is the free variable\n",
"x = -3 - 2*y\n",
"z = np.full((len(y), ), -2) # z is a constant"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "4145126158fe4b3d8111e8b470633547",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize = (12,9))\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"ax.plot(x, y, z, lw = 3, color = 'red')\n",
"ax.scatter(special_solution[:,0], special_solution[:,1], special_solution[:,2], s = 200)\n",
"ax.set_title('General Solution and Special Solution of the Linear Sytem', size= 16)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: A Symbolic Solution "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a system where all right-hand side values are indeterminate:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align}\n",
"x + 2y - 3z &= a\\\\\n",
"4x - y + 8z &= b\\\\\n",
"2x - 6y - 4z &= c\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define $a, b, c$ as SymPy objects, then extract the augmented matrix"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 2 & -3 & a\\\\4 & -1 & 8 & b\\\\2 & -6 & -4 & c\\end{matrix}\\right]$"
],
"text/plain": [
"⎡1 2 -3 a⎤\n",
"⎢ ⎥\n",
"⎢4 -1 8 b⎥\n",
"⎢ ⎥\n",
"⎣2 -6 -4 c⎦"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a, b, c = sy.symbols('a, b, c', real = True)\n",
"A = sy.Matrix([[1, 2, -3, a], [4, -1, 8, b], [2, -6, -4, c]]); A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can immediately achieve the symbolic solution by using ```.rref()``` method."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & \\frac{2 a}{7} + \\frac{b}{7} + \\frac{c}{14}\\\\0 & 1 & 0 & \\frac{16 a}{91} + \\frac{b}{91} - \\frac{10 c}{91}\\\\0 & 0 & 1 & - \\frac{11 a}{91} + \\frac{5 b}{91} - \\frac{9 c}{182}\\end{matrix}\\right], \\ \\left( 0, \\ 1, \\ 2\\right)\\right)$"
],
"text/plain": [
"⎛⎡ 2⋅a b c ⎤ ⎞\n",
"⎜⎢1 0 0 ─── + ─ + ── ⎥ ⎟\n",
"⎜⎢ 7 7 14 ⎥ ⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢ 16⋅a b 10⋅c ⎥ ⎟\n",
"⎜⎢0 1 0 ──── + ── - ──── ⎥, (0, 1, 2)⎟\n",
"⎜⎢ 91 91 91 ⎥ ⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢ 11⋅a 5⋅b 9⋅c⎥ ⎟\n",
"⎜⎢0 0 1 - ──── + ─── - ───⎥ ⎟\n",
"⎝⎣ 91 91 182⎦ ⎠"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A_rref = A.rref(); A_rref"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, we can substitute values of $a$, $b$ and $c$ to get a specific solution. "
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & \\frac{31}{14}\\\\0 & 1 & 0 & - \\frac{16}{91}\\\\0 & 0 & 1 & - \\frac{69}{182}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 31 ⎤\n",
"⎢1 0 0 ── ⎥\n",
"⎢ 14 ⎥\n",
"⎢ ⎥\n",
"⎢ -16 ⎥\n",
"⎢0 1 0 ────⎥\n",
"⎢ 91 ⎥\n",
"⎢ ⎥\n",
"⎢ -69 ⎥\n",
"⎢0 0 1 ────⎥\n",
"⎣ 182 ⎦"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vDict = {a: 3, b: 6, c: 7}\n",
"A_rref = A_rref[0].subs(vDict);A_rref # define a dictionary for special values to substitute in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Polynomials "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider this question : How to find a cubic polynomial that passes through each of these points $(1,3)$,$(2, -2)$ ,$(3, -5)$, and $(4, 0)$.\n",
"\n",
"The form of cubic polynomial is \n",
"\n",
"\\begin{align}\n",
"y=a_0+a_1x+a_2x^2+a_3x^3\n",
"\\end{align}\n",
"We substitute all the points:\n",
"\n",
"\\begin{align}\n",
"(x,y)&=(1,3)\\qquad\\longrightarrow\\qquad \\ 2=a_0+3a_1+9a_2 +27a_3 \\\\\n",
"(x,y)&=(2,-2)\\qquad\\longrightarrow\\qquad 3=a_0+a_1+a_2+a_3\\\\\n",
"(x,y)&=(3,-5)\\qquad\\longrightarrow\\qquad 2=a_0-4a_1+16a_2-64a_3\\\\\n",
"(x,y)&=(4,0)\\qquad\\longrightarrow\\qquad -2=a_0+2a_1+4a_2+8a_3\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It turns to be a linear system, the rest should be familiar already."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The augmented matrix is "
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 3\\\\1 & 2 & 4 & 8 & -2\\\\1 & 3 & 9 & 27 & -5\\\\1 & 4 & 16 & 64 & 0\\end{matrix}\\right]$"
],
"text/plain": [
"⎡1 1 1 1 3 ⎤\n",
"⎢ ⎥\n",
"⎢1 2 4 8 -2⎥\n",
"⎢ ⎥\n",
"⎢1 3 9 27 -5⎥\n",
"⎢ ⎥\n",
"⎣1 4 16 64 0 ⎦"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = sy.Matrix([[1, 1, 1, 1, 3], [1, 2, 4, 8, -2], [1, 3, 9, 27, -5], [1, 4, 16, 64, 0]]); A"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & 0 & 4\\\\0 & 1 & 0 & 0 & 3\\\\0 & 0 & 1 & 0 & -5\\\\0 & 0 & 0 & 1 & 1\\end{matrix}\\right], \\ \\left( 0, \\ 1, \\ 2, \\ 3\\right)\\right)$"
],
"text/plain": [
"⎛⎡1 0 0 0 4 ⎤ ⎞\n",
"⎜⎢ ⎥ ⎟\n",
"⎜⎢0 1 0 0 3 ⎥ ⎟\n",
"⎜⎢ ⎥, (0, 1, 2, 3)⎟\n",
"⎜⎢0 0 1 0 -5⎥ ⎟\n",
"⎜⎢ ⎥ ⎟\n",
"⎝⎣0 0 0 1 1 ⎦ ⎠"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A_rref = A.rref(); A_rref"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1, 0, 0, 0, 4],\n",
" [0, 1, 0, 0, 3],\n",
" [0, 0, 1, 0, -5],\n",
" [0, 0, 0, 1, 1]], dtype=object)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A_rref = np.array(A_rref[0]); A_rref"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last column is the solution, i.e. the coefficients of the cubic polynomial."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 4., 3., -5., 1.])"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"poly_coef = A_rref.astype(float)[:,-1]; poly_coef"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Cubic polynomial form is:\n",
"\\begin{align}\n",
"y = 4 + 3x - 5x^2 + x^3\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we have the specific form of the cubic polynomial, we can plot it"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"x = np.linspace(-5, 5, 40)\n",
"y = poly_coef[0] + poly_coef[1]*x + poly_coef[2]*x**2 + poly_coef[3]*x**3"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAHWCAYAAABucBCJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFz0lEQVR4nO3deXhU1f3H8fchCVsg7CiCikaJiIKiQgWU4Ar8tCqLIoobFlDcgVZr60prrQItooCCBRRFFi2IgKVI3NCiuIAbCIiCsikEiREIyfn9cWaYhCRkm5ubufN5Pc883HPnzsyXy/LJnHPuucZai4iIiMS+an4XICIiItGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgIhKqBtjnjXGbDPGfJZvX0NjzGJjzNehXxsU89ruxpjVxpi1xpi7o1GPiIhIPIrWN/UpQPeD9t0NLLHWHg8sCbULMMYkAE8CPYATgSuNMSdGqSYREZG4EpVQt9a+Bew4aPclwNTQ9lTg0iJe2gFYa61db63dB8wIvU5ERETKyMsx9cOstZsBQr82LeKY5sDGfO1NoX0iIiJSRok+f74pYl+R69YaYwYBgwBq1qx52lFHHeVlXXEvLy+PatU0j9JrOs/e0zn2XlDPce3vviNhzx4Afm3enP3Jyb7VsmbNmh+ttU1KOs7LUN9qjGlmrd1sjGkGbCvimE3AkfnaLYAfinoza+3TwNMAaWlpdvXq1dGuV/LJyMggPT3d7zICT+fZezrH3gvkOc7NhZSUSPvDD+Hww30rxxjzbWmO8/JHq3nAtaHta4G5RRzzAXC8MeYYY0x1oF/odSIiIv5ZvRqys932EUf4GuhlEa1L2l4E3gPSjDGbjDEDgb8B5xtjvgbOD7UxxhxhjFkAYK3dD9wCvA58Ccy01n4ejZpERETKbcWKyHb79v7VUUZR6X631l5ZzFPnFnHsD0DPfO0FwIJo1CEiIhIV+UP9tNP8q6OMgjezQUREpKI++iiyrVAXERGJUXl58PHHkXYMdb8r1EVERPL7+mvIynLbhx3mJsrFCIW6iIhIfgdPkjNFLalSNSnURURE8ovRSXKgUBcRESkoRifJgUJdREQkIi+vYKjH0CQ5UKiLiIhErF8PP//sths3hiOPPPTxVYxCXUREJCyGJ8mBQl1ERCQihifJgUJdREQkIoYnyYFCXURExLE2pifJgUJdRETE2bABdu502w0aQMuWflZTLgp1ERERiPlJcqBQFxERcWJ8khwo1EVERJwYnyQHCnURERE3Se7g7vcYpFAXERHZuBF++slt16sHqan+1lNOCnUREZH839JPPTUmJ8mBQl1ERCQQk+RAoS4iIhKISXKgUBcRkXgXkElyoFAXEZF49/33sG2b265TB44/3t96KkChLiIi8S1/1/upp0K12I3G2K1cREQkGgIySQ4U6iIiEu8CMkkOFOoiIhLvAjJJDhTqIiISzzZvdg+A2rUhLc3feipIoS4iIvErf9f7KadAQoJvpUSDQl1EROJXgCbJgUJdRETiWYAmyYFCXURE4lmAJsmBQl1EROLVtm2waZPbrlkTWrf2t54oUKiLiEh8yt/13q4dJCb6V0uUKNRFRCQ+BWySHCjURUQkXgVskhwo1EVEJF4FbJIcKNRFRCQe/fQTfPut265RA9q08beeKFGoi4hI/Mnf9d62LSQl+VdLFCnURUQk/gSw6x0U6iIiEo8COEkOFOoiIhKP9E1dREQkAHbuhPXr3XZSEpx0kr/1RJFCXURE4svHH0e2Tz7ZzX4PCIW6iIjEl4B2vYNCXURE4k0Al4cN8zTUjTFpxphP8j1+NsbccdAx6caYXfmOuc/LmkREJM7ln/kesG/qnt6Sxlq7GjgFwBiTAHwPvFLEoW9bay/yshYRERF27YKvv3bbiYlu4ZkAqczu93OBddbabyvxM0VERCI++SSy3aaNu496gFTmzWP7AS8W89yZxphPgR+A4dbazw8+wBgzCBgE0KRJEzIyMryqU4CsrCyd40qg8+w9nWPvxdI5bjFzJseFtjc3a8bqGKm7tIy11vsPMaY6LrDbWGu3HvRcCpBnrc0yxvQE/mmtPf5Q75eWlmZXr17tXcFCRkYG6enpfpcReDrP3tM59l5MneOrroIXXnDb48bB0KH+1lNKxpgV1trTSzqusrrfewAfHRzoANban621WaHtBUCSMaZxJdUlIiLxJMCT5KDyQv1Kiul6N8Ycbowxoe0OoZp+qqS6REQkXuzeDeFe3mrVoF07f+vxgOdj6saY2sD5wOB8+4YAWGsnAH2Am4wx+4FfgX62MsYEREQkvnz6KYTj5cQToXZtf+vxgOehbq3NBhodtG9Cvu1xwDiv6xARkTgX4JXkwrSinIiIxIcAryQXplAXEZH4EPBJcqBQFxGReLBjB3zxhdtOSIBTTvG1HK8o1EVEJPjefjsySe6006BOHX/r8YhCXUREgi//ynGxslBOOSjURUQk+BTqIiIiAbBjh7tGHdx4eufO/tbjIYW6iIgE28Hj6Skp/tbjIYW6iIgEW5x0vYNCXUREgk6hLiIiEgBxNJ4OCnUREQmyOBpPB4W6iIgEWRx1vYNCXUREgkyhLiIiEgBxNp4OCnUREQmqOBtPB4W6iIgEVZx1vYNCXUREgkqhLiIiEgBxOJ4OCnUREQmiOBxPB4W6iIgEURx2vYNCXUREgkihLiIiEgBxOp4OCnUREQmaOB1PB4W6iIgETZx2vYNCXUREgkahLiIiEgBxPJ4OCnUREQmSOB5PB4W6iIgESRx3vYNCXUREgkShLiIiEgBxPp4OCnUREQmKOB9PB4W6iIgERZx3vYNCXUREguLNNyPbCnUREZEYtXMnfPKJ247T8XRQqIuISBBoPB1QqIuISBBoPB1QqIuISBAo1AGFuoiIxDqNpx+gUBcRkdim8fQDFOoiIhLb1PV+gEJdRERim0L9AIW6iIjELo2nF6BQFxGR2KXx9AIU6iIiErvU9V6A56FujNlgjFlljPnEGPNhEc8bY8xYY8xaY8xKY0x7r2sSEZGAUKgXkFhJn9PNWvtjMc/1AI4PPToC40O/ioiIFE/j6YVUhe73S4Bp1nkfqG+MaeZ3USIiUsVpPL2Qygh1C/zHGLPCGDOoiOebAxvztTeF9omIiBRPXe+FVEb3e2dr7Q/GmKbAYmPMV9bat/I9b4p4jT14R+gHgkEATZo0ISP/H6ZEXVZWls5xJdB59p7Osff8OsenzZ9P3dD2yoYN2aE/Z4y1hfLTuw8z5gEgy1r7eL59E4EMa+2LofZqIN1au7m490lLS7OrV6/2uty4lpGRQbp+8vWczrP3dI6958s53rkTGjVy3e8JCbBjR6C7340xK6y1p5d0nKfd78aYZGNM3fA2cAHw2UGHzQOuCc2C/w2w61CBLiIiovH0onnd/X4Y8IoxJvxZL1hrFxljhgBYaycAC4CewFogG7je45pERCTWaTy9SJ6GurV2PdCuiP0T8m1bYKiXdYiISMAo1ItUFS5pExERKT1dn14shbqIiMQWjacXS6EuIiKxRV3vxVKoi4hIbFGoF6uy1n4PDmthyxb46ivIynJjOQ0b+l2ViEh80Hj6ISnUi7N/P6xf78L7yy8L/rprV+S4xETo1g169YJLL4XDD/etZBGRwNN4+iEp1LOyYPXqSGCHw/vrryEnp+TX798Pixe7x803u58ae/d2IX/UUd7XLyIST9T1fkjxEerWwtatRX/r3rix5NcfLCUFTjgB8vLgw3y3iLcW3nnHPe68E04/3YV7797QqlX0fj8iIvFKoX5IwQ11a+Hpp2HqVBfgmZllf4/mzV14t25d8NdmzcCE7kOzcSO88grMmVOwWwhc4H/4Ifzxj9CmTeQbfNu2kdeLiEjpaDy9RMEM9e3b4brrYMGCko9NTITjjisc3mlppRurOfJIuO0299i6FebOhZdfhiVLXNd82Oefu8dDD0FqaiTgzzgDqukiBBGREuX/4tS+vcbTixC8UM/IgKuugh9+KLi/bl0X2AeHd2oqJCVF57MPOwwGDXKPnTth/nz3Df7112HPnshx69bB3//uHi1awGWXuZDv0sX99CkiIoWp671EwQn1/fvh4YfdI38X+F13uccRR1Rul3eDBjBggHtkZcGiRS7g58937bBNm+CJJ9yjSRM3g75XLzjnHKhevfLqFRGp6hTqJQpGv+/GjS4EH3ooEuiNG8Nrr8GoUW5s3M8x7Dp1oE8fePFFNzTw6qtw/fWFr2/fvh2eeQZ69HDj9s88U/AHFBGReJV/PL1aNdezKYXEfqjPmwennOLGWsK6dYNPP4WePX0rq1g1a8JFF8Gzz7pFbP77X7jppsLXt+/Y4brxe/WCH3/0p1YRkapC16eXSuyG+t69cPvtcMklLgDB/fT28MPumvEjjvC3vtJISoJzz4WnnoLvv3eXwt11lxtnD/v3v+Hkk924vIhIvFLXe6nEZKhX27cPzjwTxo6N7GzRwv2h/+lPsTnZrFo1d3nGqFFuMZyh+W4xv2ULdO/ufoj59Vf/ahQR8YtCvVRiMtSTv/0WPv44suOSS1x3+1ln+VdUNNWuDePGuTkBTZtG9o8d6y6BW7myUsux1tKuXTumTp0KwNq1axk8eDDt2rUjISGB9HL+A5s9ezadOnWiUaNG1KxZk7S0NEaOHMm+ffsOHDN06FAGDhwYjd+GiMQqjaeXWkyG+oFxlerV3azxV14J5k1VevaEVavcGHzY55+7YB892q1oVwlmzpzJzp076d+/f6iEz1mwYAGtWrWiVQVWyvvpp5/o1q0bkyZNYuHChdxwww385S9/4a677jpwzIgRI5g+fTpr166t8O9DRGKUxtNLLXYvaWvVCl56yU2SC7KmTd1kwIkT3Xj7r7/Cvn0wbJhbXGfqVDe730Njx45lwIABJIWu57/44ou55JJLAOjTpw8/lnMi3+DBgwu0u3Xrxs8//8yTTz7JE088gTGGli1b0qVLF8aPH8+oUaMq9hsRkdikrvdSi8lv6jkpKbBiRfADPcwYGDIEPvrIraIUtmSJm0Q3e7ZnH7127VqWLVtGnz59Duyr5uEKeI0aNSrQ/Q7Qu3dvpk+fTl4l9UyISBWjUC+1mAz1PYcf7q79jjcnnADvvQf33BO57n7nTujb1133vnt31D9yyZIlJCcn065du6i/d1hubi7Z2dm88847jB07lptuugmTb12BTp06sXXrVlatWuVZDSJSRWk8vUxiMtTjWvXq8Ne/up9c89/adcoU13Px3ntR/bgVK1bQunVrT7+dJycnk5yczFlnnUXXrl157LHHCjzfpk0bEhISWL58uWc1iEgVpfH0MlGox6qzz3Yz/kOT1wBYv95dAfDAAwVvJlMBW7ZsoXHjxlF5r+IsW7aMt99+m1GjRjF37lxuueWWAs8nJiZSv359tmzZ4mkdIlIFqeu9TBTqsax+fZg+3T3q1XP7cnPhwQdduK9bV+GP2LNnDzVq1Kjw+xxK+/bt6dKlC3fddRdjx45l/PjxrDuo9ho1arAn/01xRCQ+KNTLRKEeBP37u2/tZ58d2ff++647/l//qtD68Q0bNiSzPPeiL6f2oYmA33zzTYH9mZmZNAziZYsiUjyNp5eZQj0ojj4a3ngDHnnE3SMe3N3gbrjBTaT76adyvW1aWlqhgPXSu+++C8AxxxxzYN/27dvJzs6u0DXxIhKDNJ5eZgr1IElIgLvvdt/S09Ii++fMgbZt3c1jyqhz58589913bN++/cC+7OxsZs+ezezZs/n+++/Zvn37gXZ2dvaB4zIyMjDGkJG/+yyf7t278/jjj7Nw4UL+85//cP/99zNs2DCuuOIKUlNTDxz34YcfYoyhU6dOZa5fRGKYut7LLHYXn5HinXaau6Z92DCYMMHt++EHOP98uPNON3u+Zs1SvVV6ejoNGzZk0aJFDBgwAIBt27bRt2/fAseF29988w0tW7YEOBDwTfMvdZvPGWecwZQpU9iwYQOJiYkce+yxPPLIIwwZMqTAcYsWLaJr1640atSoVDWLSEAo1MtM39SDqnZtGD/e3bu9SZPI/jFjoGNHKGWXevXq1bn66quZMWPGgX0tW7bEWlvkIxzoAP/73/9IT0/nxBNPLPK9H374YT777DOysrLIzMzko48+4tZbbz2wch24a9jnzJmj9d9F4o3G08tFoR50F13k1o/Pf2/5lSuha1d3CVwpjBgxgoyMDNasWVOmj162bFmBddzLY9asWdSqVYt+/fpV6H1EJMZoPL1cFOrx4LDDYP58d9/28OVpGze67qxSXPbWokULJk+ezObNm8v0sYsXL+biiy8uR8ER1lomT55MYqJGikTiirrey0X/U8YLY+Cmm+C44+C3v4U9eyLBvnSp238Ifn1TvvLKK335XBHxmUK9XPRNPd6cf74bZw9PlNu0yf2D0a1NRaSq0Hh6uSnU49F557nu+Fq1XPv77yE9nXVLNnDzzW7o6pxzupKSAjffHJWF6URESk/j6eWmUI9X555bINgXfn8ybc9vyqRnLLt3g7WG3bth0iR3ifvChT7XKyLxQ13v5aZQj2fnnAPz57Ouxon0YTbZtjY5+02BQ3JyIDsb+vTRN3YRqSQK9XJTqMe7c85h1IX/IYekQx6Wk+MucRcR8ZTG0ytEoS48v7Q5OVQ/5DE5OfDcc5VUkIjEL42nV4hCXcjKiu5xIiLlpq73ClGoC3XqRPc4EZFyy3/jKYV6mSnUhauvhqRDD6mTlGgJ3c9FRMQbX3/tlrUGt/rlWWf5W08MUqgLw4aVItT3/8qd/1e2td9FRMpk9uzIdvfuULeuf7XEKIW6kJrq/i3Vrl043JPYR21+YTa9Sb22C3z2mT9Fikjw5Q/1Pn38qyOGKdQFgB493M3bBg1yk02NsaSkwKBeP7EyuRM9WATbt7tr2xXsIhJt69fDRx+57erVoYI3g4pXCnU5IDUVxo2DXbvgjTfeZNcuGDenGamLJ0S6wbZvh27dIuNeIiLRMGdOZPuCC6BePf9qiWGehrox5khjzFJjzJfGmM+NMbcXcUy6MWaXMeaT0OM+L2uScjjzTPjPfyLXi/74owv2lSv9rUtEgkNd71Hh9Tf1/cAwa21r4DfAUGPMiUUc97a19pTQ4yGPa5Ly+M1vCgb7Tz+5rvhPP/W3LhGJfd9+C8uXu+3ERHd7aCkXT0PdWrvZWvtRaHs38CXQ3MvPFA917Fg42M89N7Kko4hIeeTvej/vPGjQwL9aYlyljakbY1oCpwL/K+LpM40xnxpjFhpj2lRWTVIOHTvC4sWR8a5wsH/8sb91iUjsyt/13revf3UEgLHhNXa9/BBj6gBvAn+x1r580HMpQJ61NssY0xP4p7X2+CLeYxAwCKBJkyanzZw50/O641lWVhZ1DrGEXN2vvqLd8OEk/vILADl16/Lp44+T1apVZZUYCCWdZ6k4nWPvVeQc19i+nTMvvxwAW60a7778Mvs1Sa6Qbt26rbDWnl7ScZ6HujEmCZgPvG6tHV2K4zcAp1trfyzumLS0NLt69eroFSmFZGRkkF7SEo0ffgjnnw+Zma7doIFb4rF9e6/LC4xSnWepEJ1j71XoHI8dC7eH5lCff74b4pNCjDGlCnWvZ78bYDLwZXGBbow5PHQcxpgOoZp+8rIuiZLTT3chHh7/2rnTdcXrcjcRKS3Neo8qr8fUOwMDgHPyXbLW0xgzxBgzJHRMH+AzY8ynwFign62MMQGJjtNOKxjsmZlu0YitW30tS0RiwObN8M47brtaNbj0Ul/LCYJEL9/cWvsOYEo4Zhwwzss6xGPt27tg79rV3Z/122/hssvgjTegZk2/qxORqurllyP3Tu/aFZo29beeANCKchId7dvDjBnup22A996DgQMj/2BFRA6mWe9Rp1CX6Pm//4NRoyLtF16AkSP9q0dEqq6tW+Gtt9y2Ma53TypMoS7RdfvtMHhwpH3ffaDLD0XkYP/+N+Tlue2zzoLDD/e1nKBQqEt0GQNPPOFmwYdde21kCUgREYBZsyLbmvUeNQp1ib6kJPcPNrwQzZ49bi3n777zty4RqRq2b4eMjEi7Vy/fSgkahbp4o0EDmD8fGjZ07a1b3aVuu3f7W5eI+G/uXMjNddudO0Nz3RIkWhTq4p3jj3c3akgMXTm5ciVcdVXkH7OIxCctOOMZhbp4Kz0dJkyItF99Fe6+27dyRMRnO3bAkiWRtrreo0qhLt4bOBBGjIi0H38cJk3yrx4R8c/cubB/v9vu2BGOOsrfegJGoS6V45FH3GS5sJtugqVL/atHRPyhrndPKdSlciQkwPTp0K6da+/fD717w9df+1uXiFSezExYvDjSVqhHnUJdKk+dOm5MPbzIxM6dcNFFboxNRILv1VchJ8dtn346tGzpazlBpFCXynXkkTBvXuRGL2vWuDWfw//QRSS41PXuOYW6VL4zzoBp0yLtN96AW27RzV9Eguznn+H11yPt3r39qyXAFOrij7594eGHI+2nn4Z//MO3ckTEY/Pnw969bvuUU+C443wtJ6gU6uKfe+91i9GEDRvm/uGLSPCo671SKNTFP8a469U7dXJta+HKK93KcyISHFlZsHBhpK17p3tGoS7+qlkTXnklMgs2K8utEb9li69liUgULVjgbuwEcPLJkZs9SdQp1MV/TZu6S13q1nXt776DSy+FX3/1tSwRiRLdZrXSKNSlajjpJHjpJagW+iv5v//BDTdoRrxIrPvlF/dNPUyh7imFulQdPXrAmDGR9owZ8NBD/tUjIhW3aBFkZ7vt1q3hxBP9rSfgFOpStdx6q1sXPuyBB+DFF30rR0QqKP+sd02Q85xCXaoWY+Cf/4Tzz4/su/56eP99/2oSkfL59deCl6mq691zCnWpepKSYOZMOOEE1967Fy65BL791t+6RKRsXn/dXdECbsb7SSf5W08cUKhL1VS/vvsJv1Ej1962zV3qFv4PQkSqvoMXnDHGv1rihEJdqq7UVHj5ZffNHWDVKhg4UDPiRWLB3r3u5k1h6nqvFAp1qdrOPhsmToy0Z850Y+4iUrUtXgy7d7vtY491672L5xTqUvVdfz0MGRJpDx8Ob7/tXz0iUrKDZ72r671SKNQlNvzjH9Chg9vOzYXLL4fNm30tSUSKsW8fzJ0baavrvdIo1CU21KjhfvJv3Ni1t2xxwZ6T429dIlLYkiWQmem2jz4aTjvN13LiiUJdYseRR7pV5sJLyb7zDvzhD/7WJCKFada7bxTqElvOPRdGjoy0x4xxa8aLSNWQkwP//nekra73SqVQl9jzhz+4xWjCBg6EL77wrx4RicjIgB073PaRR0LHjr6WE28U6hJ7qlWDqVPhuONc+5dfoFcv+Plnf+sSkYK3We3dW13vlUyhLrGpXj23ME2tWq69erW79E0L04j4Z/9+eOWVSFtd75VOoS6x6+ST4ZlnIu2XX4ZRo/yrRyTevfUW/Pij227WDM4809964pBCXWLbVVfBLbdE2n/4gxvTE5HKl3/We+/ekStVpNLojEvsGzUq8o0gLw+uuAK+/97fmkTiTW6u6y0LU9e7LxTqEvuqV3eTc5o2de1t29yylPv2+VuXSDx55x3YutVtH3YYdOlS7KHWWtq1a8fUqVMLPff9999Tp04djDFkRfGujFOmTMEYU+gxYcKEAscNHTqUgQMHRu1zK1ui3wWIREXz5u569fPOc98Y3nvPrRE/dqzflYnEh/xd7716QUJCsYfOnDmTnTt30r9//0LPjRgxgjp16vDLL794USVvvPEGtcITbIFjjz220OefcMIJ3HPPPRwXvsImhuibugRHejo88kik/cQTMH26b+WIxI28PJgzJ9Iuoet97NixDBgwgKTwbZVD3n77bRYtWsTw4cO9qBKAM844g9/85jcHHk3DPXwhLVu2pEuXLowfP96zGrykUJdgGT7cfUsIGzTI3YddRDyT8vnnkRssNW7sbplcjLVr17Js2TL6HBT8ubm53Hrrrdx33300Dt/jwSe9e/dm+vTp5OXl+VpHeSjUJViMgX/9C9LSXDs724X8rl3+1iUSYE3eeivSuOwySCx+ZHfJkiUkJyfTrl27AvsnTJjAnj17GDp0qFdlApCamkpiYiJpaWlMnDixyGM6derE1q1bWRWDXwgU6hI8KSluFm5ysmuvXQvXXuu6CEUkuvLyaPLmm5F2376HPHzFihW0bt2aavkud/vpp5/485//zOjRowt1yUdLs2bNePjhh3nuued49dVX6dixI0OGDGHMmDGFjm3Tpg0JCQksX77ck1q8pFCXYDrxRJg8OdKeOxf+/nf/6hEJmHXr4OabIaVuHrW3byGFTG6uMYl1R6Yf8nVbtmwp1L1+77330rFjR3r27OlZvRdeeCF/+tOfuOCCC+jRowfTpk3j8ssvZ+TIkYW62RMTE6lfvz5btmzxrB6vKNQluK64Au64I9K+9153n2cRqZCFC6FtW5g0CXZnJ2Kpxm7qMWnftbQ9LYmFC4t/7Z49e6hRo8aB9ueff86zzz7L/fffT2ZmJpmZmWRnZwOwa9cufv31V89+H3369GHHjh1s2LCh0HM1atRgz549nn22VzwPdWNMd2PMamPMWmPM3UU8b4wxY0PPrzTGtPe6Jokjf/975HrZvDzo1w82bvS3JpEYtm6dm9yene3usppfjk0kO9s9v25d0a9v2LAhmZmZB9pff/01OTk5nHnmmTRo0IAGDRocGFdv0aIFt956q0e/kwhTxE1nMjMzadiwoeefHW2eXqdujEkAngTOBzYBHxhj5llr898nswdwfOjRERgf+lWk4pKSYOZMaN8etmxx61L36ePWqM73bUFESmfUqMJhfrCcHBgzBsaNK/xcWloa77333oF2ly5dWLp0aYFjFi1axKOPPsqCBQsKXUceTXPmzKFx48YcffTRBfZv376d7OxsWrVq5dlne8Xrb+odgLXW2vXW2n3ADOCSg465BJhmnfeB+saYZh7XJfGkWTMX7OHFMJYvhzvv9LcmkRj1/POlC/Xnniv6uc6dO/Pdd9+xfft2ABo3bkx6enqBxwknnADAWWedRVr4ShYgIyMDYwwZxdzfYdq0aSQmJvLtt98Weq537948+uijLFy4kPnz5zNgwABeeukl7rvvvgKT9gA+/PBDjDF06tTp0L/RKsjrFeWaA/n7OjdR+Ft4Ucc0BzbnP8gYMwgYBNCkSZNi/1AlOrKysgJ3jlsMHsxxTz3lGuPH82W9emy98EJfawriea5qdI6jKyurK1DyPdJ377ZkZLxZaL8xhpSUFMaMGcMFF1xQ5Gu/+uorwC1Gk3/1t/fffx+ADRs2FPln+sUXX5Cbm8t7773HN998U+C5mjVr8uSTT7Jt2zastRx99NHcc889nHzyyYXea9KkSbRr1y4mL2nDWuvZA+gLTMrXHgA8cdAxrwFd8rWXAKcd6n1btWplxVtLly71u4Toy8uz9vLLrXV3Xbe2Zk1rP/7Y15ICeZ6rGJ3j6KpbN/JP6FCPlJTi3+O2226zPXv2LPNn33fffTY9Pb0C1Zds//79tnnz5va5557z9HPKCvjQliJ3ve5+3wQcma/dAvihHMeIVJwx7jK31q1de88ed3vInTv9rUskhlx9tZuqcihJSTBgQPHPjxgxgoyMDNasWVOmz162bBl33XVXmV5TVrNmzaJWrVr069fP08/xiteh/gFwvDHmGGNMdaAfMO+gY+YB14Rmwf8G2GWt3XzwG4lERZ06bmGaOnVce/1697+PFqYRKZVhwyApyR7ymKSkQ09badGiBZMnT2bz5rL9V7948WIuvvjiMr2mrKy1TJ48mcRDrIpXlXka6tba/cAtwOvAl8BMa+3nxpghxpghocMWAOuBtcAzwM1e1iTCCSfAlCmR9muvwV//6ls5IrEkNRVmP7aB2vxCEgVvb5yUBLVruxu2paYe+n369etH165dPay0fK688krOPsTa9VWd59epW2sXWGtbWWtTrbV/Ce2bYK2dENq21tqhoedPttZ+6HVNIvTu7W7+EnbfffD66/7VIxJDemz5FytpyyCeJiUpG2MsKSnu/kkrV0KPHn5XGL+0opzEr0cecbdrBTe3p39/KGJlKRHJZ+9emDSJVNYzjlvZ9cIC3njjTXbtctell/QNXbylUJf4lZgIM2bAEUe49o4dbmGaGFwaUqTSTJ8euc1qs2bg8Ri3lI1CXeLbYYfBrFmRW0WuWAG33eZvTSJVVV4ePP54pH3HHVqZsYpRqIt06gSjR0fazzwDzz7rXz0iVdVrr8GXX7rtunVh8GB/65FCFOoiALfc4sbUw26+GT76yL96RKqi/LcvHjwY6tXzrxYpkkJdBNzCNE8/DSed5Np797oZ8jt2+FuXSFXx3nvwzjtuOykJbr/d33qkSAp1kbDkZLcwTUqKa2/YAFddpYVpRAAeeyyy3b8/tGjhXy1SLIW6SH7HHw/TpkXaixbBQw/5V49IVbBmDfz735F2/jUepEpRqIsc7JJL4O67I+0HH4QFC/yrR8Rvo0a5tRwA/u//IsNUUuUo1EWK8vDDcO65kfZVV7l14kXizdatMHVqpD1ihH+1SIkU6iJFSUyEF1+MjBtmZrqJc7/+6mtZIpXuiSfcxFGADh0ghtdFjwcKdZHiNGni7kwRvs/kJ5+4S93soe9QJRIYWVnw5JOR9u9/764UkSpLoS5yKB07wj//GWlPmeIWpxGJB5MmuV4qgOOOg0sv9bMaKQWFukhJhgyBa66JtG+9FT74wL96RCpDTg6MGRNpDxsGCQn+1SOlolAXKYkxMH48tG3r2vv2ufH1H3/0ty4RL82cCd9957abNIFrr/W3HikVhbpIadSu7RamCS+LuXGjW4AjN9ffukS8YG3BJWFvvRVq1fKvHik1hbpIaaWmwvPPR9qLF8P99/tXj4hXFi+GlSvddu3aboKoxASFukhZXHQR/OlPkfZf/gKvvupfPSJeyP8t/cYboVEj/2qRMlGoi5TVAw/ABRdE2gMGwNq1vpUjElUffQRLlrjthAS4805/65EyUaiLlFVCAkyfDkcd5dq7drmJc9nZ/tYlEg35b9xy+eXQsqVvpUjZKdRFyqNxY7cwTfXqrr1ypbu/tBamkVj2zTdu1nuYloSNOQp1kfI64wwYNy7Sfv55d+mbSKwaPTpyq+HzzoNTT/W3HikzhbpIRdx4I1x/faR9xx3w/vu+lSNSbj/+CJMnR9q//71/tUi5KdRFKsIYtzZ2+BtNTg706QPbtvlbl0hZPfVU5IZFp5zivqlLzFGoi1RUrVowZw40aODa338P/frB/v3+1iVSWtnZ7m5sYSNG6MYtMUqhLhINxxzjZsSH/yNcurTg9ewiVdnUqZFlj48+Gvr29bceKTeFuki09OhRcIW5Rx+FV17xrx6R0sjNhccfj7Tvuityu2GJOQp1kWj6859duIddey2sWeNfPSIlefllWL/ebTdoADfc4G89UiEKdZFoqlbNXdoWXrBj927o1Qt+/tnXskSKZG3BxWaGDoU6dfyrRypMoS4SbQ0buolzNWq49uefw5VX6o5uUvW8+SZ88IHbrlHD3Y1NYppCXcQL7dvD009H2gsWwPDh/tUjUpT8N2657jpo2tS3UiQ6FOoiXrnmGrjnnkj7H/+AiRN9K0ekgFWrYOFCt20MDBvmbz0SFQp1ES+NHOnG1MOGDoX//te/ekTC8s9479ULjj/ev1okahTqIl6qVg2mTXPd8QC5uay7bDg3X7mTlBQ455yupKTAzTfDunX+lipxZONGeOGFSFs3bgkMhbqI15KTYd48OOIIFtKdtlnvMmlGMrt3g7WG3bth0iRo2zbSGyriqX/+M7Li4dlnQ8eO/tYjUaNQF6kMzZuz7qlF9GE22SSTQ/UCT+fkuJU6+/TRN3bxWGZmwbkdunFLoCjURSrJqNdPJieh1iGPycmBMWMqqSCJTxMmQFaW2z7xxIKLJUnMU6iLVJLnn4ec3EP/k8vJgeeeq6SCJP7s3eu63sNGjHDzPiQw9KcpgWStpV27dkydOhWAWbNm8dvf/pbmzZtTp04dTjvtNF588cVKrSn85ciZBfwWaA7UAU4DXix03NChQxk4cGBllShB9/zzsGWL2z7iCOjf3996JOoU6hJIM2fOZOfOnfQP/ac1evRo6tSpw5gxY5g3bx7dunWjf//+PJH/dpMeK7j65mhcmI8B5gHdgP7AEwWOGzFiBNOnT2ft2rWVVqcEVF5ewSVh77gDqlcv9nCJTYl+FyDihbFjxzJgwACSQnebevXVV2ncuPGB58855xx++OEHRo8eza2VtDTm1Ve7We45OQCvAo3zPXsO8AMwmgEDIvW0bNmSLl26MH78eEaNGlUpdUpAzZ8Pq1e77ZQUGDTI33rEE/qmLoGzdu1ali1bRp8+fQ7syx/oYaeeeirbtm2rtLqGDct/R8vC9cCpwDbuHPJrgb29e/dm+vTp5OXleVyhBFr+JWGHDIF69fyrRTyjUJfAWbJkCcnJybRr1+6Qxy1btowTTzyxkqqC1FSYPRtq1y58u+ok9pHAWxxPIqkPXuO6SkM6derE1q1bWbVqVaXVKgGzbBm8+67bTkqC22/3tx7xjEJdAmfFihW0bt2aaoeY1btkyRLmzp3L0KFDK7Eyd/XQypWu5zMlBYyxpKRAj5Onksd8/sjPLvnvv//Aa9q0aUNCQgLLly+v1FolQPKPpV99tZskJ4GkMXUJnC1bthTZ3R62YcMG+vfvzyWXXMJ1111X7s/ZtWsXmzdvLvG4E044oUA7NRXGjXOPjIw3admyJR07/olLjj2W69avdweNHAlpaXD11SQmJlK/fn22hGcti5TF6tUwd26krbsFBppnoW6MeQy4GNgHrAOut9ZmFnHcBmA3kAvst9ae7lVNEh/27NlD7dq1i3xux44d9OjRg6OOOornn3++Qp8za9Ysfve735V4nLW22Od+/vnnSD3//S/06weLFrknBw6EY46Bzp2pUaMGe/bsqVC9EqdGjYLw38GLLnILzkhgedn9vhg4yVrbFlgD3HOIY7tZa09RoEs0NGzYkMzMzEL7s7Ozueiii9i3bx+vvfYaycnJFfqcG2+8EWttiY/iZGdn88c//jFST716MGNG5D/dffvgssvgm2/IzMykYcOGFapX4tCWLRBaqwHQkrBxwLNQt9b+x1obumMA7wMtvPoskfzS0tL45ptvCuzbv38/ffv25euvv2bhwoU0bdrUp+oK1rNp06aC9dSr5y49Cg8fbN/O9u7dyc7OplWrVv4VLLFp7Fj3wyHAb34DXbr4W494rrImyt0AFHf/KQv8xxizwhijCyelwjp37sx3333H9u3bD+y7+eabWbBgAX/+85/ZsWMH77///oHH3r17DxyXkZGBMYaMjAxPawzXM2DAgML1HHEE/PvfBxYG+XDNGgzQqUMHT2uSgNm9G8aPj7RHjABj/KtHKkWFxtSNMf8FDi/iqXuttXNDx9wL7AemF/M2na21PxhjmgKLjTFfWWvfKuKzBgGDAJo0aeL5f7rxLisrK2bPsTGGlJQUxowZwwUXXADAvHnzALi9iEt5XnzxRQ4/3P01fv/99wE3mc7L33+4nnHjxjFu3Lgi6zls+HBa//WvLAK6Ar/efDMZlbRQTpDE8t/lijjq+ec5NjQMld2iBcvr1QOPzkO8nuMqqTRjguV9ANcC7wG1S3n8A8Dwko5r1aqVFW8tXbrU7xIq5LbbbrM9e/Ys8+vuu+8+m56e7kFFRSvpPO+/5x7bHOxzbqqTtU8+WTmFBUis/10ul48+sjYpyf2dAWsnTvT04+LyHFcy4ENbihz1rPvdGNMd+APwW2ttdjHHJBtj6oa3gQuAz7yqSeLHiBEjyMjIYM2aNWV63bJly7jrrrs8qqrsZp10ErXq1KFfeMdtt8HixX6WJFVddra7UYtbjxg6dIDrr/e3Jqk0Xo6pjwPq4rrUPzHGTAAwxhxhjFkQOuYw4B1jzKfAcuA1a+0iD2uSONGiRQsmT55cquvI81u8eDEXX3yxR1WVnTWGyXPmkHh66MKQ3Fzo2xe+/NLfwqTqGj4cvvrKbScnw/TphZcwlMDy7Dp1a+1xxez/AegZ2l4PHHotT5Fy6tevX8kHVXFXXnml2zj5ZPeNa9Mm2LXLXW/8v/9FZsmLALz6asHJcWPHwnFF/lcsAaVlYkViQbNmMG+eWzgeYP166NUL8s3clzi3ZQvccEOk3bu3ut3jkEJdJFaceqrrSg1flvT223DjjQVu/iJxKi8PrrsOfvzRtZs3h6ef1iVscUihLhJLLr0UHn000n7+ebecbG6ubyVJFTBuHLz+uts2BqZNA61AGJcU6iKxZvhwd5u3sClTFOzxbNWqgsu/Dh8O55zjXz3iK4W6SKwxxk2Gyj9+OnWqayvY48uePe7ytfDcilNPhYcf9rcm8ZVCXSQWVasGzzzjxtTDpk1zE6MU7PHj7rvhs9DSHrVquTkXNWr4W5P4SqEuEquqVYOJEyH/7V+fe85NmFKwB9+iRfDPf0bao0dD69b+1SNVgkJdJJZVqwYTJhQcY3/+ebj2WgV7kG3f7n54C7v4Yhg82LdypOpQqIvEumrV3Bh7/v/Up0+Ha66B/fuLf53EJmvdxMitW137sMNg8mRdviaAQl0kGKpVg6eegptuiux74QUYMEDBHjQTJ7qV48KmTIEmTXwrR6oWhbpIUFSrBk8+WTDYZ8xQsAfJV19B/hsO3XYbdO/uXz1S5SjURYLEGBfsQ4dG9s2YAVddpWCPdfv2ucvXfv3VtU86qeBCRCIo1EWCxxh44gm45ZbIvpkzC96OU2LPn/8MH3/stmvUcMMrNWv6W5NUOQp1kSAyxt2h69ZbI/tmzVKwx6o33oDHHou0H33U3blP5CAKdZGgMsZdx3z77ZF9s2fDlVcq2GPJjh3uSgZrXfvCCwv+sCaSj0JdJMiMgTFj4I47IvvmzIF+/RTsscBatwbB99+7duPG8K9/uUmRIkXQ3wyRoDPGrTZ2552RfS+/DFdc4SZf+cxaS7t27Zg6dSoAs2fPplOnTjRq1IiaNWuSlpbGyJEj2edTrVOmTMEYU+gxYcKEqLxu6NChDBw4sLg3cT+EhU2eDM2aVfB3JEGW6HcBIlIJjIFRoyIBD/DKK3D55W4SXfXqvpU2c+ZMdu7cSf/+/QH46aef6NatGyNGjKB+/fosX76cBx54gC1btjBu3Djf6nzjjTeoVavWgfaxxx4bldeNGDGCE044gXvuuYfjjjsu8sTatQW72YcMgd/+tnzFS9xQqIvEC2Pg8ccjAQ8wd67vwT527FgGDBhAUlISAIMPWu60W7du/Pzzzzz55JM88cQTGJ9WTjvjjDOoU6dO1F/XsmVLunTpwvjx4xkV/nPJyXGXIf7yi2unpUX+zEQOQd3vIvHEGDeLesSIyL65c6FPn8jtOyvR2rVrWbZsGX369DnkcY0aNfKt+70y9O7dm+nTp5OXl+d2PPQQLF/utpOS3OVrtWv7V6DEDIW6SLwxxl0S9fvfR/a9+qovwb5kyRKSk5Np165doedyc3PJzs7mnXfeYezYsdx0002+fUsHSE1NJTExkbS0NCZOnBjV13Xq1ImtW7eyatUqePtt+OtfI0+OHAnt21e0fIkT6n4XiUfGwN/+Fgl4gPnzoXdvNzGrku7JvWLFClq3bk21ImZzJycnszf0Q8Y111zDY/mv065EzZo14+GHH6ZDhw7k5uby4osvMmTIELKzs7kz/+TDCryuTZs2JCQksDwjg3ajR0P4G3u3bjB8uJe/PQkYhbpIvDIGHnnEXR71yCNu32uvQa9eLtgrYbWyLVu20Lhx4yKfW7ZsGdnZ2SxfvpyHHnqIW265haeeeqpcn7Nr1y42b95c4nEnnHBCoX0XXnghF1544YF2jx492Lt3LyNHjuT2228v8geSsr4uMTGR+vXrs2XyZPjuO7ezQQOYOlWXr0mZKNRF4pkx8Je/uF/DXb4LFrhgf/llz4N9z5491C5mrLh9qMu5S5cuNG7cmGuvvZZhw4aRmppa5s+ZNWsWv/vd70o8zoYXeClBnz59mDlzJhs2bCj1LPiSXlcjL489q1ZFdkycCEceWer3FgGNqYuIMW7c9t57I/sWLoRzz4VvvvH0oxs2bEhmZmaJx4UD/pty1nPjjTdirS3xUVblHeMv9LoNG8jcuZOG4fb110PfvuV6b4lvCnURccH+8MPupiFhy5ZBu3bw3HORJUqjLC0trVRB/e677wJwzDHHeFJHWc2ZM4fGjRtz9NFHV/x1+/ez/YoryAZaAaSmuuV9RcpB3e8i4hgDDz4IdevCPfdAbi7s3u3WHX/tNRg/3o3zRlHnzp156KGH2L59O02aNAGge/funHfeeQcmj7377ruMGjWKK664okDXe0ZGBt26dWPp0qWkp6dHta78evfuTYcOHWjbti25ubm89NJLvPTSS4wdO7bAuPi0adO44YYbWLduHUcffXSpX8ff/saHy5djgE7VqsH06e7PQKQcFOoiEmGMu4b97LPh6qvdqmYAL73kvrlPmwZRDND09HQaNmzIokWLGDBgAOAWa5kyZQobNmwgMTGRY489lkceeYQhQ4YUeG12djYATZs2jVo9RUlLS+PZZ59l48aNWGs58cQTmTZt2oF6w/Ly8sjNzT3QjV+q102fDg88wCKgK9DowQehY0dPfz8ScKUZZ6pqj1atWlnx1tKlS/0uIS5U6fO8e7e1Awda6zrf3cMYa3//e2v37o3ax9x22222Z8+eZX7dfffdZ9PT00s8rkqe46wsa6+/3lqw+8E2B/tcq1bW7t/vd2XlUiXPccAAH9pS5KPG1EWkaHXqwKRJbhZ8w9AULmvh73+H3/wGvvwyKh8zYsQIMjIyWLNmTZlet2zZMu66666o1FCpPv0UTj/d3W0NmAXUSkqi38KFkJDgb20S8xTqInJol10Gq1bBBRdE9n38sVvl7MknKzyJrkWLFkyePLlU15Hnt3jxYi6++OIKfXalshaeesp1r3/1VWR3585MnjePxDJcGidSHIW6iJTsiCPcZW7/+Edktbk9e+CWW+Cii2Dr1gq9fb9+/ejatWvF66yqdu50q/UNHRpZird2bZgyhSvffpuzu3f3tz4JDIW6iJROtWpw++3wwQdw8smR/QsWuPb8+f7VVpUtWwannOJudRvWrh2sWAHXXusmJ4pEiUJdRMrm5JPdHcTyj2dv3w4XXww33QShWelxLzfXrdJ39tmRpV/BfVt//30oYklakYpSqItI2dWs6e7vvXix65oPmzDBjbWvWOFfbVXBli1w4YVulb7cXLevQQM36XDcuEpZV1/ik0JdRMrvvPNg5Uo3Xhy2erWbHf+3v0UCLZ68/rrrXl+yJLKvc2f45BM36VDEQwp1EamYRo1g1ix49ll3GRzA/v1uVbpzzoFvv/W3vsqSkwN/+AN07w7btrl9xrhv6xkZcNRRvpYn8UGhLiIVZ4y7Ccknn7hv6WFvveW+tb7wgm+lVYpvvoGzznLX8Icdfrgbnhg5EhK1eKdUDoW6iERPaiq8/TY88EBkIZVdu+Cqq9yjFHdkizmzZrnZ7f/7X2TfhRe6RWbOPde3siQ+KdRFJLoSE+H++124519Q5YUX3Lf2t97yr7Zoys6GwYPh8svh55/dvsREeOwxd5mfx2vSixRFoS4i3jjzTNcdf/31kX3ffeduCNO9Ozz9dGTsOdZ8/jl06OB+D2HHHAPvvgvDh7tr+kV8oL95IuKdunXdBLpZsyK3bbXWzRAfPBiaNYOuXWHsWNi40d9aS8NaeOYZOOMMF+xhl1/uls7t0MG/2kRQqItIZejTx60f37Nnwf15ea47/vbb3ezwjh3h0Ucjt3ytSnbtgn79YNAg+PVXt69WLRfyM2ZAvXr+1ieCQl1EKkvz5vDaa+4St3/8w80WP3iJ1OXL4e674fjjoW1bePBB98NABW8aUyG//OKuOT/1VJg5M7K/TRu3ZO6NN2qpV6kyFOoiUrmOOsp9M3/rLdi82a1Cd8EFhS/7WrXKzaJv2xbS0lzYL19eqoBftw5uvhlSUuCcc7qSkuLa69YV8wJr3fj+m2/CxIlwxx1u3L9lS3ft/XnnucvWwgYNcrW0aVPOkyDiDc8unjTGPAD8Dtge2vVHa+2CIo7rDvwTSAAmWWv/5lVNIlLFHHaYG1sfPBh27HA3hZkzx425h+9mBvD1165b/tFH4cgj3cpsvXu7ldoOugf5woWutz8nxz3AsHu3uzX81KmW2U9soUfTFe72p19+Gfl1586S601JcW/Ut29UT4NItHi9IsIYa+3jxT1pjEkAngTOBzYBHxhj5llrv/C4LhGpaho2hGuucY+sLHdZ2Msvuy77rKzIcRs3uol1Y8e6y8YuvRR69YJu3Vi3sTp9+hR9TxkX8oY+A1NYye2ksr50dSUkuOvvO3Z0wwHHHBOV366IF/xe5qgDsNZaux7AGDMDuARQqIvEszp13Izyyy93921fvNgF/Ny5Bb9Rb9vmLit7+mmoX59RyVPIye4BVC/2rXNIYgx3Mo5bCz6RnOzunNa6dcFfjzsOqhf/fiJVidehfosx5hrgQ2CYtfbg/q3mQP7rWDYBHT2uSURiSc2a7rauF1/svm6/+abron/lFdi6NXJcZibPZ6aTc4hAB8ihOs8lXMe4QV8UDO8WLTThTWKesRWYVWqM+S9weBFP3Qu8D/wIWOBhoJm19oaDXt8XuNBae2OoPQDoYK096EdoMMYMAgYBNGnS5LSZ+WehStRlZWVRJ3xzDvGMznMF5OZS74svaPzWWzR5+21qbt1KNXKxpZj/a4zljTferIQi44P+HnuvW7duK6y1p5d0XIVCvbSMMS2B+dbakw7afybwgLX2wlD7HgBr7SOHer+0tDS7evVqj6oVgIyMDNLT0/0uI/B0nqPEWvj4Y1K6nMzuX5NKPDwlxV12LtGhv8feM8aUKtQ9u6TNGNMsX/My4LMiDvsAON4Yc4wxpjrQD5jnVU0iElDGQPv2XH1dEkklZHpSEgwYUDlliVQ2L69T/7sxZpUxZiXQDbgTwBhzhDFmAYC1dj9wC/A68CUw01r7eXFvKCJyKMOGUapQv/POyqlHpLJ5NlHOWlvkz8LW2h+AnvnaC4BC16+LiJRVairMnn3wdepOUpJ7zJ7tjhMJIq0oJyKB0qMHrFzpFn1LSXGT4lJSXHvlSve8SFAp1EUkcFJTYdw4NxnujTfeZNcu19Y3dAk6hbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgEj06o2NMS8BaaFmfSDTWntKEcdtAHYDucB+a+3pXtUkIiISZJ6FurX2ivC2MWYUsOsQh3ez1v7oVS0iIiLxwLNQDzPGGOBy4ByvP0tERCSeGWuttx9gzNnA6OK61Y0x3wA7AQtMtNY+Xcxxg4BBAE2aNDlt5syZHlUsAFlZWdSpU8fvMgJP59l7Osfe0zn2Xrdu3VaUZni6QqFujPkvcHgRT91rrZ0bOmY8sNZaO6qY9zjCWvuDMaYpsBi41Vr71qE+Ny0tza5evbrcdUvJMjIySE9P97uMwNN59p7Osfd0jr1njClVqFeo+91ae14JRSQCvYDTDvEeP4R+3WaMeQXoABwy1EVERKQwry9pOw/4ylq7qagnjTHJxpi64W3gAuAzj2sSEREJJK9DvR/wYv4dxpgjjDELQs3DgHeMMZ8Cy4HXrLWLPK5JREQkkDyd/W6tva6IfT8APUPb64F2XtYgIiISL7SinIiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgERIVC3RjT1xjzuTEmzxhz+kHP3WOMWWuMWW2MubCY1zc0xiw2xnwd+rVBReoRERGJZxX9pv4Z0At4K/9OY8yJQD+gDdAdeMoYk1DE6+8GllhrjweWhNoiIiJSDhUKdWvtl9ba1UU8dQkww1q711r7DbAW6FDMcVND21OBSytSj4iISDzzaky9ObAxX3tTaN/BDrPWbgYI/drUo3pEREQCL7GkA4wx/wUOL+Kpe621c4t7WRH7bFkKK6KOQcCgUHOvMeaziryflKgx8KPfRcQBnWfv6Rx7T+fYe2mlOajEULfWnleOD98EHJmv3QL4oYjjthpjmllrNxtjmgHbDlHH08DTAMaYD621pxd3rFScznHl0Hn2ns6x93SOvWeM+bA0x3nV/T4P6GeMqWGMOQY4HlhezHHXhravBYr75i8iIiIlqOglbZcZYzYBZwKvGWNeB7DWfg7MBL4AFgFDrbW5oddMynf529+A840xXwPnh9oiIiJSDiV2vx+KtfYV4JVinvsL8Jci9t+Yb/sn4NxyfPTT5XiNlI3OceXQefaezrH3dI69V6pzbKyt0Pw1ERERqSK0TKyIiEhAxFSoG2O6h5adXWuM0epzHjDGPGuM2aZLBr1jjDnSGLPUGPNlaJnl2/2uKWiMMTWNMcuNMZ+GzvGDftcUVMaYBGPMx8aY+X7XElTGmA3GmFXGmE9KmgUfM93voWVm1+Am1G0CPgCutNZ+4WthAWOMORvIAqZZa0/yu54gCl2+2cxa+5Expi6wArhUf5ejxxhjgGRrbZYxJgl4B7jdWvu+z6UFjjHmLuB0IMVae5Hf9QSRMWYDcLq1tsS1AGLpm3oHYK21dr21dh8wA7fMrESRtfYtYIffdQSZtXaztfaj0PZu4EuKXnFRysk6WaFmUugRG99gYogxpgXwf8Akv2sRJ5ZCvbRLz4rEDGNMS+BU4H8+lxI4oW7hT3CLWi221uocR98/gN8DeT7XEXQW+I8xZkVoddVixVKoR33pWRE/GWPqAHOAO6y1P/tdT9BYa3OttafgVrTsYIzRcFIUGWMuArZZa1f4XUsc6GytbQ/0AIaGhkmLFEuhXtqlZ0WqvNA47xxgurX2Zb/rCTJrbSaQgbsNtERPZ+C3ofHeGcA5xpjn/S0pmKy1P4R+3YZbG6aou54CsRXqHwDHG2OOMcZUx92vfZ7PNYmUWWgS12TgS2vtaL/rCSJjTBNjTP3Qdi3gPOArX4sKGGvtPdbaFtbalrj/j9+w1l7tc1mBY4xJDk2oxRiTDFwAFHt1UsyEurV2P3AL8DpuYtHM0HK0EkXGmBeB94A0Y8wmY8xAv2sKoM7AANw3m09Cj55+FxUwzYClxpiVuC8Ei621uuRKYtFhwDvGmE9x91B5zVq7qLiDY+aSNhERETm0mPmmLiIiIoemUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCD+HybkWTH7fzhLAAAAAElFTkSuQmCC\n",
"text/plain": [
"