{ "metadata": { "name": "", "signature": "sha256:a5146b953e6ce4107eb6016b9e4a4a34b2153b7ea3de0a799042c5e879c6ce01" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "

# Introduction to NumPy & SciPy - the core numerical scientific packages

\n", "

\n", "

\n", "

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

This tutorial is prepared by ACM Student Chapter of King Abdullah University of Science and Technology (KAUST).
\n", "Material has been adapted from the following:

\n", "\n", "

- David Ketcheson's [Introduction to NumPy and Matplotlib](http://nbviewer.ipython.org/github/ketch/AMCS252/blob/master/2_Introduction_to_Numpy%20and_Matplotlib.ipynb)

\n", "

- Software Carpentry's [Numerical analysis with NumPy](http://nbviewer.ipython.org/github/geocarpentry/2014-01-30-mit/blob/gh-pages/lessons/vignettes/Numerical%20analysis%20with%20NumPy.ipynb)

\n", "

- J.R. Johansson's [Numpy - multidimensional data arrays](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-2-Numpy.ipynb) and [SciPy - Library of scientific algorithms for Python](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-3-Scipy.ipynb#SciPy---Library-of-scientific-algorithms-for-Python)

\n", "

- Official [Tentative NumPy tutorial](http://www.scipy.org/Tentative_NumPy_Tutorial)

\n", "\n", "**Prerequisites:** Introduction to Python\n", "

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## NumPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python includes a package for numerical computation called numpy, which is an essential tool for scientific computing. \n", "It is based on multidimensional arrays (vectors and matrices) and provides a large library of functions and operators that work efficiently on these objects.\n", "In this way, any numerical algorithm is expressed as operations on arrays and in many cases runs almost as quickly as the equivalent C code.\n", "\n", "Like MATLAB, NumPy relies on BLAS and LAPACK for efficient linear algebra computations.\n", "\n", "NumPy is used in:
\n", "Numerical Analysis, Linear algebra, Solution of differential equations, Statistical analysis and many more...\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why to use NumPy?\n", "\n", "Consider the following list and suppose we want to perform some simple algebraic operations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "price = [5.99, 10.25, 2.0, 40.99, 5.60, 63.49]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "price*2" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 3, "text": [ "[5.99, 10.25, 2.0, 40.99, 5.6, 63.49, 5.99, 10.25, 2.0, 40.99, 5.6, 63.49]" ] } ], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "price+0.5" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "can only concatenate list (not \"float\") to list", "output_type": "pyerr", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprice\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mTypeError\u001b[0m: can only concatenate list (not \"float\") to list" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clearly, if we wanted to double the value of each element and add half, this cannot be done easily. Instead, we would have to use a loop:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "for i in range(len(price)):\n", " price[i] = price[i]*2\n", " price[i] = price[i]+0.5\n", "print price" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[12.48, 21.0, 4.5, 82.48, 11.7, 127.48]\n" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "But, there should be a more convinient way..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Getting started\n", "\n", "We start, by importing the NumPy module. Import statements like this are the typical way of getting access to functions in Python. The most general way to import a module is the following:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from numpy import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above imports all functions included in the package and can call them simply by their name." ] }, { "cell_type": "code", "collapsed": false, "input": [ "sin(3*pi/2)*exp(2)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "-7.3890560989306504" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead we can import the numpy module and tell Python that we want to refer to numpy by the short abbreviation *np*. This is recommended since we keep track of which module the functions belong to." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Arrays\n", "\n", "NumPy's main object is the homogeneous multidimensional array. \n", "It is a table of elements - all of the same type - indexed by a tuple of positive integers.\n", "\n", "NumPy's array class is called *ndarray* and the package includes basic utility functions for creating and manipulate arrays in a similar fashion to MATLAB.\n", "For example we can create a one-dimensional array by:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "price = np.array([5.99, 10.25, 2.0, 40.99, 5.60, 63.49])\n", "print type(price)\n", "price = price*2+0.5\n", "print price" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "[ 12.48 21. 4.5 82.48 11.7 127.48]\n" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are other ways to create one-dimensional arrays:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.linspace(-1,1,5)\n", "print x" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[-1. -0.5 0. 0.5 1. ]\n" ] } ], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "y = np.arange(0,1,0.2)\n", "print y" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0. 0.2 0.4 0.6 0.8]\n" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Exercise 7\n", "\n", "Change the inputs of the *linspace* and *arange* function until you understand exactly what they do. Using both ways, create a vector of starting from -1.5 to 4.5 inclusive with a step of 0.5. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# type your solution here" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Solution" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#%load solutions/1D_array.py" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can always find out what a function does by typing its name, followed by a question mark:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.arange?" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resize the help window (to get it out of your way) by dragging the divider. Double click on the divider to close it." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Multidimensional arrays\n", "\n", "We can create multidimensional arrays, like the following:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = np.array([[1,2.4,-13],[4.1,5,0],[7.2,8,9]])\n", "print \"A = \", A\n", "print type(A)\n", "print A.shape\n", "print A.dtype" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "A = [[ 1. 2.4 -13. ]\n", " [ 4.1 5. 0. ]\n", " [ 7.2 8. 9. ]]\n", "\n", "(3, 3)\n", "float64\n" ] } ], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see all different attributes of the array object type its name followed by a fullstop and press *tap*. Note that some attributes are basically functions and require parentheses (with or without arguments), e.g. *conjugate* and others don't, e.g. *size*:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = np.array([[1,2+4j,3j],[4j,5,6-5j],[7j,8,9+3j]]) # dtype=complex\n", "print A.conjugate()\n", "print A.size" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 1.-0.j 2.-4.j 0.-3.j]\n", " [ 0.-4.j 5.-0.j 6.+5.j]\n", " [ 0.-7.j 8.-0.j 9.-3.j]]\n", "9\n" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also define the shape of the array:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.random.uniform(0,1,size=(5,5)) # 5 x 5 matrix with elements from a uniform distribution" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 17, "text": [ "array([[ 0.18894904, 0.4016313 , 0.30167486, 0.04452608, 0.29880108],\n", " [ 0.85862765, 0.18898635, 0.46421297, 0.6303044 , 0.56185042],\n", " [ 0.78056999, 0.82705839, 0.63247392, 0.07986216, 0.64999473],\n", " [ 0.07852346, 0.27204531, 0.85521149, 0.94891731, 0.7810846 ],\n", " [ 0.47172365, 0.57284027, 0.57618836, 0.30150177, 0.54945194]])" ] } ], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Special arrays\n", "\n", "The function *zeros* creates an array full of zeros, the function *ones* creates an array full of ones and *eye* an identity array. Function *diag* is used for a diagonal array with specified diagonal elements. A matrix with random entries is generated by the *random* function." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print np.zeros((3,4))\n", "print np.ones((2,3))\n", "print np.eye(3)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0. 0. 0. 0.]\n", " [ 0. 0. 0. 0.]\n", " [ 0. 0. 0. 0.]]\n", "[[ 1. 1. 1.]\n", " [ 1. 1. 1.]]\n", "[[ 1. 0. 0.]\n", " [ 0. 1. 0.]\n", " [ 0. 0. 1.]]\n" ] } ], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.diag(x,0)\n", "print np.diag([2,3,4])\n", "print np.random.random((2,2))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[-1. 0. 0. 0. 0. ]\n", " [ 0. -0.5 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0.5 0. ]\n", " [ 0. 0. 0. 0. 1. ]]\n", "[[2 0 0]\n", " [0 3 0]\n", " [0 0 4]]\n", "[[ 0.89073408 0.92386659]\n", " [ 0.43499535 0.1673178 ]]\n" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that there are two pairs of parentheses used when calling the *zeros*, *ones* and *random* command, because we are passing as argument (e.g. the Python tuple *(3,4)*), which is the shape of the resulting array." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reshape\n", "\n", "Manipulation of arrays can be done in many ways: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = np.arange(1,10)\n", "print \"A = \", A\n", "A = A.reshape((3,3))\n", "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "A = [1 2 3 4 5 6 7 8 9]\n", "[[1 2 3]\n", " [4 5 6]\n", " [7 8 9]]\n" ] } ], "prompt_number": 20 }, { "cell_type": "code", "collapsed": false, "input": [ "B = arange(15).reshape((3,5))\n", "print \"B = \", B\n", "print B.ravel() # flatten the array" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "B = [[ 0 1 2 3 4]\n", " [ 5 6 7 8 9]\n", " [10 11 12 13 14]]\n", "[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14]\n" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Indexing\n", "\n", "We can slice numpy arrays just as other programming languages, but the indexing is a little different. See the examples below:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print A\n", "print A[0,0]\n", "print A[1:3,0:2]\n", "print A[1:3,[0,1]]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1 2 3]\n", " [4 5 6]\n", " [7 8 9]]\n", "1\n", "[[4 5]\n", " [7 8]]\n", "[[4 5]\n", " [7 8]]\n" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the fundamental difference from MATLAB is that the indexing starts from zero and is exclusive of the second index. Hence, in the above example 1:3 refers to the second and third row and 0:2 to the first and second column." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can index in some clever ways too: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1 2 3]\n", " [4 5 6]\n", " [7 8 9]]\n" ] } ], "prompt_number": 23 }, { "cell_type": "code", "collapsed": false, "input": [ "print A[0,:] # row selection\n", "print A[:,0] # column selection" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[1 2 3]\n", "[1 4 7]\n" ] } ], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "print A[:,1:] # all rows, from 2nd column until the end\n", "print A[-1,:] # last row, all columns" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[2 3]\n", " [5 6]\n", " [8 9]]\n", "[7 8 9]\n" ] } ], "prompt_number": 25 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can change array's elements:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A[:,2] = [.3,.4,9.12]\n", "print A\n", "print A.dtype" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1 2 0]\n", " [4 5 0]\n", " [7 8 9]]\n", "int64\n" ] } ], "prompt_number": 26 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the elements of A remain integers. That is because A has type *int*. The following does the desired change:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = A.astype(float)\n", "A[:,2] = [.3,.4,9.12]\n", "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 1. 2. 0.3 ]\n", " [ 4. 5. 0.4 ]\n", " [ 7. 8. 9.12]]\n" ] } ], "prompt_number": 27 }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is an easy way to reverse a vector:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "a = arange(10)\n", "print a\n", "print a[: :-1] # reverse a" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[0 1 2 3 4 5 6 7 8 9]\n", "[9 8 7 6 5 4 3 2 1 0]\n" ] } ], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also index an array using boolean arrays:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "b = a>5\n", "print b\n", "a[b] = 0 # All elements of 'a' higher than 5 become 0\n", "print a" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[False False False False False False True True True True]\n", "[0 1 2 3 4 5 0 0 0 0]\n" ] } ], "prompt_number": 29 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Array and matrix multiplication" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print a**2\n", "print 10*sin(np.pi/a[1:5]) # avoiding division with zero!\n", "print A*A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0 1 4 9 16 25 0 0 0 0]\n", "[ 1.22464680e-15 1.00000000e+01 8.66025404e+00 7.07106781e+00]\n", "[[ 1. 4. 0.09 ]\n", " [ 16. 25. 0.16 ]\n", " [ 49. 64. 83.1744]]\n" ] } ], "prompt_number": 30 }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, NumPy uses component-wise multiplication. For a matrix-matrix (or matrix-vector) multiplication, we use the *dot* function: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.dot(A,A)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 31, "text": [ "array([[ 11.1 , 14.4 , 3.836 ],\n", " [ 26.8 , 36.2 , 6.848 ],\n", " [ 102.84 , 126.96 , 88.4744]])" ] } ], "prompt_number": 31 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other basic operations" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = array(A,dtype='int')\n", "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1 2 0]\n", " [4 5 0]\n", " [7 8 9]]\n" ] } ], "prompt_number": 32 }, { "cell_type": "code", "collapsed": false, "input": [ "print A.sum() # sum of all elements\n", "print A.min()\n", "print A.max()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "36\n", "0\n", "9\n" ] } ], "prompt_number": 33 }, { "cell_type": "code", "collapsed": false, "input": [ "print A.sum(axis=0) # sum of each column\n", "print A.min(axis=1) # min of each row\n", "print A.cumsum(axis=1) # cumulative sum along each row" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[12 15 9]\n", "[0 0 7]\n", "[[ 1 3 3]\n", " [ 4 9 9]\n", " [ 7 15 24]]\n" ] } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Python *axis* refers to the dimensions of arrays; *axis=0* refers to the colums and *axis=1* to the rows." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stacking and splitting\n", "\n", "The *vstack* and *hstack* commands can be used to stack arrays and *hsplit* to split." ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = floor(10*random.random((2,2)))\n", "print A\n", "B = floor(10*random.random((2,2)))\n", "print B" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 4. 0.]\n", " [ 6. 7.]]\n", "[[ 1. 2.]\n", " [ 5. 3.]]\n" ] } ], "prompt_number": 35 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.vstack((A,B))\n", "print np.hstack((B,A))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 4. 0.]\n", " [ 6. 7.]\n", " [ 1. 2.]\n", " [ 5. 3.]]\n", "[[ 1. 2. 4. 0.]\n", " [ 5. 3. 6. 7.]]\n" ] } ], "prompt_number": 36 }, { "cell_type": "code", "collapsed": false, "input": [ "A = floor(10*random.random((2,12)))\n", "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 6. 7. 5. 2. 3. 6. 2. 4. 7. 5. 2. 9.]\n", " [ 8. 4. 4. 5. 7. 0. 5. 3. 3. 3. 6. 3.]]\n" ] } ], "prompt_number": 37 }, { "cell_type": "code", "collapsed": false, "input": [ "print np.hsplit(A,3) # Split A into 3 arrays\n", "print np.hsplit(A,(3,4)) # Split A after the third and the fourth column" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[array([[ 6., 7., 5., 2.],\n", " [ 8., 4., 4., 5.]]), array([[ 3., 6., 2., 4.],\n", " [ 7., 0., 5., 3.]]), array([[ 7., 5., 2., 9.],\n", " [ 3., 3., 6., 3.]])]\n", "[array([[ 6., 7., 5.],\n", " [ 8., 4., 4.]]), array([[ 2.],\n", " [ 5.]]), array([[ 3., 6., 2., 4., 7., 5., 2., 9.],\n", " [ 7., 0., 5., 3., 3., 3., 6., 3.]])]\n" ] } ], "prompt_number": 38 }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Exercise 8\n", "\n", "1. Create a vector with values ranging from 10 to 99\n", "2. Find indices of non-zero elements from [1,2,0,0,4,0]. Hint: use np.nonzero\n", "3. Declare a 8x8 matrix and fill it with a checkerboard pattern of zeros and ones, without a loop.\n", "4. Create a 5x3 matrix whose $(i,j)$ entry is equal to $i\\times j^2$ without using a loop. Hint: use np.fromfunction.\n", "5. Create a 4x4 array with numbers from a normal distribution with mean = 3 and standard deviation = 1.5. Hint: use np.random.normal. Sum the values of each row and store the result in a new array." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# type your solution here" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 39 }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Solution" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#%load solutions/array_operations.py" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 40 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's work with some real data. Ensure that the stockholm_temp.txt in the same directory with this notebook. \n", "This file contains the [daily average temperatures](http://bolin.su.se/data/stockholm/homogenized_daily_mean_temperatures.php)\n", "according to observations of Stockholm for the years 1756 - 2012.\n", "The following code checks if the file exists, opens it, prints the first lines and stores the data in an array." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import os\n", "filename = 'data/stockholm_temp.txt'\n", "if os.path.exists(filename):\n", " # look at the first 3 lines\n", " with open(filename,'r') as f:\n", " print '\\n'.join(f.readlines()[:3])\n", " # alternative look at the first head lines\n", " !head 'data/stockholm_temp.txt'\n", " \n", " data = np.loadtxt(filename) # alternative we can use getfromtxt command\n", " \n", "else:\n", " print 'Weather data does not exist, please download \"stockholm_temp.txt\" from Github.'" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "1.756000000000000000e+03 1.000000000000000000e+00 1.000000000000000000e+00 -8.699999999999999289e+00\n", "\n", "1.756000000000000000e+03 1.000000000000000000e+00 2.000000000000000000e+00 -9.199999999999999289e+00\n", "\n", "1.756000000000000000e+03 1.000000000000000000e+00 3.000000000000000000e+00 -8.599999999999999645e+00\n", "\n", "1.756000000000000000e+03 1.000000000000000000e+00 1.000000000000000000e+00 -8.699999999999999289e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 2.000000000000000000e+00 -9.199999999999999289e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 3.000000000000000000e+00 -8.599999999999999645e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 4.000000000000000000e+00 -7.700000000000000178e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 5.000000000000000000e+00 -7.200000000000000178e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 6.000000000000000000e+00 -1.600000000000000089e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 7.000000000000000000e+00 6.999999999999999556e-01\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 8.000000000000000000e+00 1.300000000000000044e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 9.000000000000000000e+00 2.399999999999999911e+00\r\n", "1.756000000000000000e+03 1.000000000000000000e+00 1.000000000000000000e+01 8.000000000000000444e-01\r\n" ] } ], "prompt_number": 41 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the data are in the form\n", "*Year\n", "- Month\n", "- Day\n", "- Temperature*\n", "we would like to express them in a better way. We define a dtype, that is the type for each column: int, int, int, float:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "dt = np.dtype([('Year', 'int16'), ('Month', 'int8'), ('Day', 'int8'), ('Temp', 'float64')])\n", "data = np.loadtxt(filename,dtype=dt)\n", "data[:10] # first 10 entries of our data" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 42, "text": [ "array([(1756, 1, 1, -8.7), (1756, 1, 2, -9.2), (1756, 1, 3, -8.6),\n", " (1756, 1, 4, -7.7), (1756, 1, 5, -7.2), (1756, 1, 6, -1.6),\n", " (1756, 1, 7, 0.7), (1756, 1, 8, 1.3), (1756, 1, 9, 2.4),\n", " (1756, 1, 10, 0.8)], \n", " dtype=[('Year', '

\n", "

\n", "

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## SciPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "NumPy is complemented by SciPy; a library that adds more MATLAB-like functionality.\n", "\n", "SciPy is a collection of mathematical algorithms and convenience functions that provides high-level commands and classes for Nowdays, everything from parallel programming to web and data-base subroutines and classes have been made available to Python All of this power is also accessible through the the mathematical libraries in SciPy NumPy and SciPy.\n", "\n", "SciPy is used in:
\n", "Numerical Analysis, Optimization, Linear algebra, Image and signal processing, Statistical analysis, Data processing, Data transformation and query and many more..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Getting started" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As with NumPy, we can access all packages in SciPy:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 53 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Special functions\n", "\n", "SciPy provides implementations of a large set of special functions, usually used in mathematics, physics and engineering. \n", "Available functions include airy, elliptic, bessel, gamma, beta, hypergeometric, parabolic cylinder, mathieu, spheroidal wave, struve, and kelvin.\n", "\n", "For a list of functions see SciPy's documention at http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special.\n", "or try" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy import special\n", "#help(special)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 54 }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example we can look at the Bessel functions.\n", "\n", "Bessel functions are a family of solutions to Bessel\u2019s differential equation with real or complex order alpha:\n", "\\begin{align}\n", " x^2 \\frac{d^2 y}{d x^2} + x \\frac{d y}{d x} + (x^2 \u2212\\alpha^2)y = 0\n", "\\end{align}" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.special import jn, jn_zeros\n", "\n", "x = 0.0\n", "\n", "n = 0 # order\n", "# Bessel function of first kind, zero order \n", "print \"J_%d(%f) = %f\" % (n, x, jn(n, x))\n", "\n", "n = 1 # order\n", "# Bessel function of first kind, first order \n", "print \"J_%d(%f) = %f\" % (n, x, jn(n, x))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "J_0(0.000000) = 1.000000\n", "J_1(0.000000) = 0.000000\n" ] } ], "prompt_number": 55 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next two lines allow plotting inside the ipython notebook. More during the next sessions..." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# The next two\n", "%pylab inline --no-import-all\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "prompt_number": 56 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now plot the Bessel functions of first kind for different values of order:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = linspace(0, 12, 100)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(x,zeros(len(x)),'--k')\n", "for n in range(4):\n", " ax.plot(x, jn(n, x), label=r\"$J_%d(x)$\" % n)\n", "ax.legend();" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXlcVNUbxp/JJZdUkEVc0yT3JTVNW9wyyzAtW0zNyjbN\nXDLLTEtbLXHL3FL7aVnumOa+awqIILIIioigILKq7AzL3Of3xxETZWCGuXcG4Xw/n/nkcM89553b\nzHPPfc973ldHEhKJRCIp/9xnawMkEolEYh2k4EskEkkFQQq+RCKRVBCk4EskEkkFQQq+RCKRVBCk\n4EskEkkFobKtDShAp9PJ+FCJRCIxE5I6U9uWqRk+yXL5mjlzps1tkJ9Pfj75+crfy1zKlOBLJBKJ\nRDuk4EskEkkFQQq+Fejdu7etTdAU+fnubeTnqzjoSuMH0gKdTseyYotEIpHcC+h0OtCMRdsyE6Uj\nkUgqLjqdyZpVYVFjQiwFXyKRlAnkE75x1LohSh++RCKRVBCk4EskEkkFQQq+RCKRVBCk4EskEkkF\nQQq+RCKRaERUVFSJbeLi4pCVlWUFa6TgSyQSiSZERkbCx8enxHZOTk5wd3e3gkVS8CUSiUQTli9f\njmHDhpXYrnLlynBzc8OaNWs0t0kKvkQikZjB+vXr4ebmhm7dusHLy6vINkFBQWjUqJHJfXbt2hUH\nDx5Uy0SjSMGXSCQSMxg2bBhq166NSZMm4Yknniiyzc6dO9G3b1+z+nVyckJERIQaJhpFCr5EIpGY\nAUkcOXKkWEH38/NDmzZtzOq3Y8eO8Pf3t9S8YpGpFSQSicQMgoKC4OjoiHr16hltk5WVdVc6hO3b\nt6NSpUo4fvw42rdvj71792L69Olo1aoVAMDe3h7h4eGa2i4FXyKRlHnUyq2mRrqeQ4cOoV+/fsW2\nMRgMhd5HR0ejTZs2cHV1xYwZMzB16lTUqVMHTZo0udWmevXqyM3NtdzAYrBI8HU63SoAbgASSbY3\n0uYXAAMAZAF4m2SAJWNKJJKKR1nKq3bw4EGMHTsWN27cwIoVK+Ds7IwOHTqgS5cut9pUrlxYWguE\nPSEhAbVq1YKdnR0GDhxYqE1qairq1q2rqe2W+vBXA3jO2EGdTvc8AFeSDwP4AMAyC8eTSCQSq5Kd\nnY0GDRogMjIS8fHxCA0NRd++fbF69Wr06dMHI0eOxPz58wud4+LigoyMjFvvw8LCEBQUhN27d6Nn\nz54AxMLu7cTFxcHV1VXTz2KR4JM8DuBGMU0GAfjjZtuTAOx0Op1xx5dEIpGUMapWrYqxY8ciODgY\n8+bNw969e1GzZk1ERUWhfv36qFy5Mq5fv17onF69esHX1/fW+/3792Pnzp0gCb1ej61bt8LZ2bnQ\nOYGBgUajftRCax9+QwAxt72/AqARgASNx5VIJBJVqFSpEr788ksAwIsvvnjr74qioFKlSgDuzlc/\nZMgQzJ0791Ykz4QJE4odQ6/Xo3bt2qhWrZqapt+FNRZt71xuMeqNc3D4Gvb2wIMPAh991BtDhvTW\n1jKJRCIpJS1btkRCQgLq1q2L2rVrFzpmZ2cHR0dHJCcnw9HRscS+NmzYgNGjR5fY7ujRozh69Ghp\nTba8pq1Op2sKYEdRi7Y6ne5XAEdJbrj5PgxAL5J3zfB1Oh3PniUiIoDt2wEPD+Dpp4Hx44FevSwy\nUSKRlHFu1ma1tRlmce3aNaxatQp16tRB+/bt0aNHj0LHSeK3337D+++/X2w/MTExOH36NAYPHmy0\njbHrY25NW60F/3kA40g+r9PpugP4mWR3I/0UKmKelgasXQu4uwvBX7AAsLe3yFSJRFJGuRcF35qo\nJfgWLdrqdLr1ALwBtNTpdDE6ne4dnU43WqfTjQYAkrsBROp0uggAywGMNbXv2rWBDz8EzpwBHngA\naN8e2LXLEmslEomkYmPxDF8t7pzh38nRo8CoUcAbbwDffqveRgyJRGJ75Ay/eMqMS0ctShJ8AEhK\nAp57DujWDVi8GLi5QC6RSO5xpOAXT5lw6VgbJyfgyBEgLAwYPhzQeBeyRCKRlCvuKcEHhG9/zx4h\n9q++CtyRskIikUgkRrjnBB8AqlUDNm0CMjKAyZNtbY1EIpHcG9yTgg8AVaoAW7YA+/cLf75EIpFI\niueeTo9sZydCNR9/HGjWDHBzs7VFEolEUna5p6J0jOHjA7zwAuDtDTz8sMqGSSQSzZFROsVTIaN0\njNG9OzBzJjBiBJCXZ2trJBKJRBAVFVVim7i4OGRlZVnBmnIi+ADw0UcibHPmTFtbIpFIJEBkZCR8\nfHxKbOfk5AR3d3crWFSOBF+nA1avBn7/XezKlUgkEluyfPlyDBs2rMR2lStXhpubG9asWaO5TeVG\n8AHA2Rn43/+AN98E7qhHIJFIJKqwfv16uLm5oVu3bvDy8iqyTVBQEBo1amRyn127dsXBgwfVMtEo\n5UrwAWDAAGDQIOCzz2xtiUQiKY8MGzYMtWvXxqRJk4xWqNq5c+et4iem4uTkhIiICDVMNEq5E3wA\nmDVLxOcfP25rSyQSSXmDJI4cOVKsoPv5+aFNmzZm9duxY0f4+/tbal6x3NNx+MaoXRv4+WdgzBgg\nIACoWtXWFkkkkvJCUFAQHB0dUa+e8fLcWVlZd5U93L59OypVqoTjx4+jffv22Lt3L6ZPn45WrVoB\nAOzt7REeHq6p7eVS8AFgyBBg1Spg/nxg6lRbWyORSCxB9406+dA50/JY/0OHDqFfv37FtjHckeQr\nOjoabdq0gaurK2bMmIGpU6eiTp06aNKkya021atXR67GGSHLreDrdCLlQteuwNChYieuRCK5N1FD\nqNXi4MGDGDt2LFJTU3Ho0CGcP38eX3zxRaE2lSsXltYCYU9ISECtWrVgZ2eHgQMHFmqTmpqKunXr\namp7ufThF9CsGfDpp0AJBeMlEonEKNnZ2WjQoAEiIyMRHx+P0NBQPP3006hTpw66dOlS5KzcxcUF\nGRkZt96HhYUhKCgIu3fvRs+ePQGIhd3biYuLg6urq6afpVwLPgB88glw9ixw+LCtLZFIJPciVatW\nxdixYxEcHIx58+Zh7969qFGjRrHn9OrVC76+vrfe79+/Hzt37gRJ6PV6bN26Fc7OzoXOCQwMNBr1\noxbl1qVTQNWqwI8/ijBNPz/gvnJ/i5NIJGpSqVIlfPnllwCAF1980aRzhgwZgrlz596K5JlQgptB\nr9ejdu3aqFatmmXGlkCFkL9XXwUqVwbWr7e1JRKJpDxhLOGbnZ0dHB0dkZycbFI/GzZswOjRo9U0\nrUgqhODrdMDcucD06YBeb2trJBJJeSAjIwNbtmyBv78/QkJC7jo+ceJEbN26tcR+YmJiYG9vj5Yt\nW2phZiHKRXpkU3nxReCJJ+QuXImkrCHTIxePWumRK5Tgnz8PPPmk+K/G0U8SicQMpOAXj8yHXwpa\nthSz/PnzbW2JRCKRWJ8KNcMHgKgo4NFHgfBwwMFB8+EkEokJyBl+8cgZfilp1kykXViwwNaWSCQS\niXWpcDN8ALh0CejSRc7yJZKygpzhF4+c4VtA06bAyy/LWb5EIqlYVMgZPiBn+RJJWULO8ItHzvAt\npGCW//PPtrZEIpFIrEOFneEDQEQE0KOHiNx54AGrDi2RSG5DzvCLR87wVcDVFejTB1i50taWSCSS\n8khUVFSJbeLi4pCVlWUFayq44APA55+Lxdu8PFtbIpFIyhORkZHw8fEpsZ2TkxPc3d2tYJEUfHTp\nArRoITNpSiQSdVm+fDmGDRtWYrvKlSvDzc0Na9as0dymCi/4ADBlCuDuDiiKrS2RSCRlnfXr18PN\nzQ3dunWDl5dXkW2CgoLQqFEjk/vs2rUrDh48qJaJRpGCD+CZZ0ShlD17bG2JRCIp6wwbNgy1a9fG\npEmTjFao2rlz563iJ6bi5OSEiIgINUw0ihR8iHz5U6YAs2fb2hKJRFLWIYkjR44UK+h+fn5o06aN\nWf127NgR/v7+lppXLFLwb/LKK8Dly8CpU7a2RCKRlGWCgoLg6OiIevXqGW2TlZUFna5wtOT27dux\na9cuTJ06FWvXrsXIkSMRFhZ267i9vT2uXLmimd2AFPxbVK4MjB8PLFxoa0skEsld6HTqvFTg0KFD\n6NevX7FtDAZDoffR0dFo06YN3NzccODAAbi5uWHo0KFo0qTJrTbVq1dHbm6uKjYao9wXMTeHd98F\nHnoIiIsD6te3tTUSieQWZWhT1sGDBzF27FhERETgzJkzCA4OxgsvvIDOnTvfalO5cmFpLRD2hIQE\n1KpVC3Z2dhg4cGChNqmpqaircWUmOcO/DXt7YPhwYNkyW1sikUjKCtnZ2WjQoAEiIyMRHx+P0NBQ\n9O3bFzt27EDDhg3xySefYO7cuYXOcXFxQUZGxq33YWFhCAoKwu7du9GzZ08AYmH3duLi4uDq6qrp\nZ5GCfwcTJgDLl8ti5xKJRFC1alWMHTsWwcHBmDdvHvbu3YuaNWti0qRJ6NatG2JiYtCsWbNC5/Tq\n1Qu+vr633u/fvx87d+4ESej1emzduhXOzs6FzgkMDDQa9aMWFTqXjjHc3ERitXfesbUlEknF4F7O\npfPDDz9g0qRJqFGjxq2/paSkYO7cufj+++9N6kOv12PatGmYb6T+qsyloyETJ4osmvfo908ikViJ\n7du3Y8KECYiNjS30dzs7Ozg6OiI5OdmkfjZs2IDRo0drYWIhpOAXwTPPAAYDcPSorS2RSCRlla1b\nt+K7777DkCFDsGnTpruOT5w4EVu3bi2xn5iYGNjb26Nly5ZamFkI6dIxwtKlwOHDgIeHrS2RSMo/\n97JLxxqo5dKRgm+E9HTgwQeB4GDAlJQYChXEpMYgLiMOde6vA/vq9qhbvS6qVqqqvbESyT2OFPzi\nkYJvBcaNE+UPv/nm7mMkcTL2JP4I/APeV7xx4doF2Fe3R8NaDZGWk4Yb+hu4kX0D7eu1xwDXARjg\nOgCPNXoMle+TWx8kkjuRgl88UvCtQGio8OdfuiSSqwFAVl4Wlvguwf8C/geCeLvj2+jfvD9aOLRA\nrftrFTo/z5CHE1dOYM+FPdh1YRfSc9Mx/anpeLPjm3LmL5HchhT84pGCbyX69AE+/BB49VViW9g2\nTNo3Cd0adsOk7pPQvVH3u/JlFIdXtBe++fcbhF8Lx1c9v8KoTqNwn06um0skUvCLRwq+ldi8GZj7\n2yXYjRyNK2lXsHjAYvRp1seiPr1jvDF5/2TUrFITv7/4OxrVNj1vtkRSHpGCXzwyDt9K2D1yBKce\n6Y421foicHSgxWIPAI83fhzHRx1Hrwd7ocuKLtgculkFSyUSiaR45AzfCCSxyHcRZh2fhf4Za1Ez\n4WlNcuz4xvrijb/fwPMPP4/5z86XLh5JhUTO8IunzLh0dDrdcwB+BlAJwG8kZ99xvDeAfwBE3vzT\nFpJ37TcuS4KvUMGHOz+ET6wPtg3dhqpZzdC+PRAdDTzwgPrjpehTMHjDYLg84II1L67B/ZXvV38Q\niaQMIwW/eMqES0en01UCsBjAcwDaABim0+laF9H0X5Kdbr5MSy5hIxQqGL1jNM4ln4PXO15oZt8M\nDRsCPXtqV+jcrpod9r2xDwbFgOfWPodUfao2A0kkkgqNpf6DbgAiSF4imQdgA4DBRbRTp/KAxpDE\n2F1jcS75HHYN34UHqv43nR89WmTR1Ipqlath4ysb0c6pHfqu6StFXyIpB0RFRZXYJi4uDllZWVaw\nxnLBbwgg5rb3V27+7XYI4HGdThek0+l263Q68wo9WgmSGL9nPIISgrB7xO67Yur79weSkgAtS05W\nuq8SfhnwC3o06oFBGwYhOy9bu8EkEommREZGwsfHp8R2Tk5OcHd3t4JFlle8MsXpdhpAY5JZOp1u\nAIBtAFoU1fDrr7++9e/evXujd+/eFppnOu5e7vCO8caRt46g9v217zpeqRLw/vvAihXazvR1Oh1+\nGfALRm4didc8XsPfr/2NKpWqaDegRCLRhOXLl2P27NkltqtcuTLc3NywZs0avPnmm8W2PXr0KI5a\nkNXRokVbnU7XHcDXJJ+7+f4LAMqdC7d3nBMFoAvJ63f83WaLtjvO78CYXWNw8r2TxcbEX70KtG0r\nFm9r1TLaTBXyDHl4aeNLsKtmhzUvrSk5eic7G4iJAWJjhaEZGWKFuWZNoE4doF07wMlJW6MlklJy\nLy3arl+/Hn/99ReSkpKwYMGCIouWBAUF4dixYxg/frzJ/b755ptYs2ZNkcfUWrS1dIZ/CsDDOp2u\nKYCrAIYCGHaHQfUAJJKkTqfrBnGTuX5nR7YiNDEU72x/BzuG7ShxA1SDBmLn7bp1wqevJVUqVcHm\nVzej35/98N2/32Fm75mFG+TmAkeOiBzO//4LBAUJAwteDzwAZGUJ4b9xQ2SBq10b6NwZ6NcPeO01\n4I6KOxKJpGSGDRuG7du344033jBaoWrnzp148cUXzerXyckJERER2pY5JGnRC8AAAOcBRAD44ubf\nRgMYffPfHwEIARAIwBtAdyP90NokZyaz+cLm/CPwD5PP2buX7NSJVBQNDbuNuPQ4NprfiNvObRN/\nOHuWnDyZdHYme/QgZ8wgDx0iMzOL70hRyIgIcsMGcsQIsk4dsn9/8q+/yLw87T+IRFIMtvj9lxZF\nUVivXj3Gx8cbbTN48GAqZorEH3/8wQ0bNhR5zNj1ufl3k/W6wm68IokX1r+Alg4tMe/ZeSafpyhA\n8+Yi5cKjj2po4G34xvpixo/94XGmDR4IjwLeekvUX2xR5FKIaWRlATt2iMT/sbHAV18BI0YAlWU2\nT4n1uZdcOoGBgXjjjTcQEhJitE3//v2xf//+Qn/bvn07KlWqhOPHj6N9+/bYu3cvpk+fjlatWgEA\nduzYgfDwcEyePPmu/sqKS+eeZYnfEiRkJuDvoX+bdd599wHvvgv89puVBD8wEN2+/BYe/pXwU++L\n+HRnCOxqqeCLr1EDGDpUvI4eFTmgv/sO+OUX4PnnLe9fIlERnUrl56hCIMihQ4fQr1+/YtsYDIZC\n76Ojo9GmTRu4urpixowZmDp1KurUqYMmTZrcalO9enXk5uZabF+xmPM4oOULVnykO5Nwhg6zHRie\nHF6q869cIe3tyYwMlQ27nfR0cuJEsl49cuFCMjub43aN40sbXjL7UdFk9u0jH3qIHD6cTEzUZgyJ\npAis+fu3lOeee47bt2/npUuXuGnTJs6aNYunTp0q1KZ///5FnhsfH89evXoVeczDw4O//vprkceM\nXR+Y6dKpcIlb9Pl6DN8yHLP7zcbDDg+Xqo+GDYEnngCKKGOpDrt3i6ialBSRlH/CBKBaNcztPxdR\nKVFYeXqlNuP27y8Wd11cgPbtgS1btBlHIrmHyM7ORoMGDRAZGYn4+HiEhobi6aefhpeXFxwcHPDw\nww8jPDy80DkuLi7IyMi49T4sLAxBQUHYvXs3evbsCUAs7N5OXFyctgu2qIAunakHp6KFQwu80+kd\ni/p5/31g9mxg1CiVDANEaOWECcChQ8DKlaL6ym3cX/l+rH95PZ5a/RSeavIUWjsVlcXCQmrWBObN\n+8/dc+IE8NNPRfr2b+TlwSs1FQEZGYjS6xGZnY3onBxkKwrySeQpCmpWqoT6Vaui/v33o8n996Nr\nrVp4rHZttKxRA/eZUUtAIrEVVatWxdixYxEcHAwvLy/s3bsXNWrUwPDhwxEVFYX9+/fj22+/LXRO\nr1694Ovri759+wIA9u/fj/T0dNSvXx96vR5bt25Fw4aF96gGBgbivffe0/SzVKhFW69oL7y6+VWE\njA1B3ep1LeorPx9o0gQ4cEDE5ltMWJgIlWzXTuzsKibQf6X/SizxW4KT753UNtHatWvA8OEiBHTj\nRihOTjiRlobNiYk4nJKCKL0ej9Wqha61a6N5tWpoVr06Hrz/ftSoVAlVdDpU1umQYTAgLjcXcbm5\niMzOhm96Ok6mpeFGfj4G1K2L15yc8FzduqhWqZJ2n0NS5rmXFm3vxMfHB9u3b8esWbNu/S0lJQVz\n587F99+bljpMr9dj2rRpmD9/fpHH1Vq0tbnvvuAFjX14+jw9Wy9uzU0hm1Trc9o0ctIkFTpat450\ndCRXrDAp3lNRFA7ZOIQf7/lYhcFLID+fYT/8wE+mTGGjf/9l25Mn+U1UFH1TU5lrMJS626t6PZdd\nucI+AQG0O36co86d45n0dBUNl9xLaP3714IpU6YwNDSUp0+f5pAhQ+46vmDBAiYlJZnU1+rVqxkW\nFmb0uLHrAzN9+DYX+luGaPw//OsjX/OFdS+ouuAZESF0Wq8vZQcGAzl9OtmsGRkUZNap17KuscG8\nBjx26VgpBy/BNEXh7uRkPhcURGdPT37xzz8M6dKF9PFRfaz4nBz+cOkSXby8OCAoiIevX1d9DEnZ\n5l4UfG9vb/7zzz+cNWsWQ0JC7jquKApXrFhRYj/R0dHctm1bsW3UEvwK4dI5l3QOT61+CgGjA9C4\nTmNV+376aeHPf/11M0/MyhLx9HFxwNatpUp78E/YP/j0wKcIGhOEGlVqmH1+UZDE9mvXMCMqCvfp\ndJjYsCFed3YWLpedO8Wixdq1YoFXZfQGA9YmJuKn6Gi4Vq+Oec2bo03NmqqPIyl73MsuHWtQZgqg\nqIVWgq9QQc/VPfF6u9cxrts41ftfvx5YtUr48k0mMREYOBBo2VIE9N9fej/8sC3D0LBWQ8ztP7fU\nfRRw4Pp1TIuKQp6i4PtmzeDm4HB3kXZPT+Dll8U6g5lbx00lV1GwNDYWP0RHY5izM75p2hT2VWQC\nufKMFPzikT58E/nf6f/xsZWPMd+Qr0n/2dmkgwMZGWniCZcukS1akF99pUp+hqTMJLrMdeGJmBOl\n7iMiK4svBAfT1ceHmxISaCjJrlOnSCcncs+eUo9pCkk5ORxz/jwbeXtz77Vrmo4lsS1a/f7LC8au\nD6QP/z9S9al0metCv1g/1fu+nfHjhX6XyNmzZOPG5M8/qzr+xpCNbL24NfV55i0mZOXn86vISDoc\nP85Zly5Rb84irJeXWMA4csQ8Y0vBwevX2cTbm6PDwpgu8/6US6TgF49agl+uXTpTDkxBclYyVg1e\npWq/dxIUJDw0ly6JvPlF4u8PuLkB7u5ACTmvzYUkXtr4EjrX74wZvWaYdI5nSgreOX8eHR94APOb\nN0fjatXMH/jwYRGrv2MH0L27+eebQWp+Pj6OiIBXaiq2tmuHttK3X66QLp3ikS6dEjiffJ4Osx0Y\nlx6nar/GePTRYjwcvr4iu+XWrZqNfznlMh1mO/DCtQvFtsvIz+eE8HDW9/Li32qkT9i1S6R/OH/e\n8r5M4I+4ODp6enJTQoJVxpNYB7V//+UNY9cHMrWCYPL+yfj8ic/h8oCLVcYrSKh2F76+Ymb/22+a\nLXICQJM6TfD5E59j3O5xRmdKfmlpeOTUKVzPz8eZrl3xkhoFUZ5/XiRdc3MDkpMt768E3nRxwb4O\nHTAlMhKfX7wIg5wVSiSmY87dQcsXVLzD77mwh66/uJrt07aElBSRYr7QpNnHRyxu7thhFRty83PZ\nbmm7uzaX5SsKZ126RGctZ8ZTp5KPPy5Wsa1AUk4O+wQEcMiZM8zK12ZBXmI9IMqlylcxL2PXjRXZ\nh29QDOj4a0fMenoWBrUcpIJlpvPWW0DHjsAnnwA4fRp47jkRszlwoNVs8Iz2xOser+PsR2dR+/7a\nuJqTgxHnzoEk/mzdunS+elNQFJGGARAlwe7T/uExR1HwdlgYruTkYHu7djJ0U1LhMNeHX+5cOn8G\n/wn76vZ4ocULVh+7wK3DkFDh6li+3KpiDwBPNnkS/Zv3x9dHv8bB69fxqL8/+trZ4dAjj2gn9oAQ\n+N9/By5fFsnWrMD9992Hta1bo1utWngyIAAxer1VxpVI7lXK1Qxfn69Hi0UtsOGVDXi88eMqWWY6\nJNC/2QXszOqD+392/2/Ga2XiMxLR/J+ZqPngMGxo2x597e2tN/jVq0DXruLON2CA1YadFxODxbGx\nONKxI5pWr261cSUSW1KhZ/hLfJegc/3ONhF7ANDFRMMj9RlsaPm1zcT+Rl4eRkXGo34jN7SJ/hl9\n7Oysa0CDBsDGjcDbbwMXL1pt2MmNG2Nyo0boHRiIqOxsq40rkdxLlBvBT9GnYLbXbMx6elbJjbUg\nKQno3x+6jyfi45D3cFvtA6sRnJGBrv7+aF2jBs480R8JN85hZ/jOkk9UmyefBGbMAF56CcjMtNqw\n4xo1wpQmTdA7MBAX7yHRJ4ErV4AzZ4Djx4Fdu4Bjx4CICKtePkkFoNy4dKYdmoaEjAT8b/D/VLTK\nRNLSgL59xSLt999j0CARgfmOZTVWzGJjYiLGXbiAha6uGF6vHgBgX8Q+jNszDiEfhmibN78oSJFo\nLT8f+PNPwIrFTpZfvYpZly/jWKdOeFDLdYtSQgIBAcDevYC3N+DjI+rLODoCdnZA7driK3X1qsit\n5+AA9OkjXk8/DTz4oK0/gaSsUCE3XsWlx7Hu7LqMTokudR+lJjub7NOHHD36Vm6cbdtEhKI1MCgK\nv7h4kU1PnGBAWtpdxweuG8jZnrOtY8ydZGaSbdqQq1ZZfeifY2L4sI8PE3JyrD62Ma5dI3/5hezY\nkWzaVNRS8PAgY2ONn6Mo5IULolTC66+LbBZ9+pCbNpFl6KNJbAQqYi6dj/d8zIl7Jpb6/FKTn0++\n9BL52mvi3zfJzSVdXETqHC1JzcvjwOBg9jp9mklGfv3hyeFW3XF8FyEhQqVCQ60+9IzISD7i58cU\nG+ffSUsTuZbs7YVoHzwoSiGUhpwccsMGsndv8R1zd7fa1gdJGcRcwb/nXTpX06+i3dJ2CB0bivq1\n6mtgmRFIYMwYIDJS5Im/I8Xx1KmAwQDMmaPN8BFZWRgcEoJednZY6OqKKsXEvX+y7xNk5mZi+QvL\ntTGmJFatAubPF7uOa6iTt98USGJiRAQCMjKwr0MH1LByGcX8fGDFCrERuX9/8d8mTdTrPzQU+Oor\nseXj++9FnIAW2x+upF3Bvoh9CIwPRGRKJC5ev4i4jDjooEOl+yqhyn1V0LhOY7RybIWWDi3RpX4X\n9G7aG9UulvmnAAAgAElEQVSryGgpralw+fAn7JmAyvdVxvxni64FqRkzZ4rVtSNHiqw/Gx4O9OwJ\nxMQAau8HOnLjBoadPYuZTZviwzsKIRfF9ezraLW4FY68dQRtndUowGsmJDByJFC9uijObkUUEm+e\nO4csRcHmtm1RyUprCRcuiI9co4aoCd+pk3ZjeXoCn34q9r79/jvQpo3lfV64dgG/nf4NuyN2Iy49\nDs80fwbdGnRD87rN0dy+ORrUagAAMNCAPEMeLqdexvnk8ziXfA7eMd4IjA9Er6a9MKjFILze7nXU\nut94jWZJ6alQgn8l7Qo6LOuAcx+dQ70H6mlkWREsXQosWAB4eQHOzkab9ewJfPwxMGSIekMvv3oV\nM6KisL5NG7Pi63/2+Rn7L+7H7hG71TPGHNLTheq5u6t7QUwgV1HwbHAwOj3wAOa7umo6Fim2IEyb\nJuYEH31knfVqUtxLp00TT5eTJhWTudUIChXsjdiLRb6L4H/VH6MeGYUhrYfg0QaPotJ95nV2I/sG\n9l3cB4+zHjgcdRgjO4zEuG7j8LDDw+YZZQokEB0tniBjYsRKd1wckJcnJhnVqgH29kCrVuJu2KpV\nkZO0e5EKJfgf7foINarUwJz+GvlNisLDA5g4UcTPPfRQsU3XrBEh6bt2WT5svqLgk4sXceDGDexo\n1w6uZrpGcg25aLu0LZY+vxTPNH/GcoNKg48PMHiwCFFp0MCqQ9/Iy8MTAQH4sEEDjG/USJMxMjPF\n9oOICFEFUo2ZtrlERgobAJHhwtSPeuzyMXyy7xMYaMDExybi9Xavo1pldSKcYlJjsOzUMvx2+jf0\nadYHPz79Ix6yL/63UyIZGcDff4sfl5eX8J917w40bQrUrw+4uABVqwJ6PZCdDVy/DoSFAWfPisfv\n1q3Fbng3N+DRR62SCkQLKkyUTnRKNO1/smdChhXT5B4+LJKhBQSY1Dwzk6xbl4yJsWzY67m5fCYw\nkM8GBvJGbm6p+9lydgvbL22vWfUvk5g5k3zmmdKvWlpAVFYW63t58Z+kJNX7jo0lO3cm33rLgqL2\nKmEwkLNmkfXrk8dKqHF/8fpFDtk4hE0WNOG64HVUVKjCZoyMnAx+e/Rb1p1dl5P2TuK1LDOrmCmK\nWPEePlxkKhw4kFy9moyIMK96XG6uKNzz2WciiqxJE/L778k49QIbFEVhRFYWV129yikREXwlJISd\n/fzYyNub9b286OzpyXqenuzg68sBQUF8PyyMC2Ni6JeaylwzfhuoKFE6H+36iJ/t/8yscywiIECI\n/eHDZp324Yfkd9+VftjwzEy28PHhxPBw5lkokoqi8KlVT/E3/98s6sci8vLIxx4jFy60yfC+qal0\n8vRkUHq6an0GBIhCZj/8oErVStXYu1eUYViy5G67FEXhMr9ldJjtwO///Z5ZuVlWsys+PZ6jd4xm\nvTn17srsWiSKQm7fTnbrRrZuTS5adEdaWgs5dYp87z3Szk6EUZUyoizPYODO5GS+dfYsG98U9mGh\noZx16RI3JCTQLzWVl7OzeUWvZ5xez1i9nqfT0rg9KYlLr1zhB2FhbOfry5r//sunAwL4a2wsE0uI\nva0Qgn817Srtf7JnfHq8yedYxMWLZIMG5ObNZp/q7y9irkuj1QeuXaOzpyeXFxeobSY+MT5sOK8h\nM3MzVevTbC5cEKGaISE2GX5dfDybnjhR4o/JFI4cER9l40bL7dKCCxfIdu3IMWP+ixyOS4/j82uf\nZ5flXXgu6ZzNbDsRc4ItF7Xka5tfY2KGEQE/coTs1ElsXti8Wdsnwxs3yNmzxcTuzTdNLlR9LiOD\nn1y4wHqenuzh789FMTEMy8ws9dNSSl4etyQmcmhICGsfO8b+gYHcmpjI/CL6qxCC/+m+Tzl+93iT\n21tEXBzZvDm5dGmpu+jUiTxwwPT2iqJwUUwM63l68sj166Ue1xivbHqFs47NUr1fs1i+nOzSRTxe\n24AvLl5kz9OnmWOBgBw6JMT+0CEVDdOAtDTy6afJIUPIvWGH6TLXhdMPTWdOvu13bmXlZnHyvsl0\nmevCfRH7/jsQF0e+8YZ4dNq82bqPTikp5IwZwh/7+efCN1sEgenpfCUkhE6envzi4kWGGWlnCZn5\n+VwXH89up06x+YkTXBQTw4zb9vyUe8FPzkym/U/21tlVm5JCPvII+c03FnWzZAk5dKhpbXMMBn4Q\nFsa2J0/yYpY2j9kFm7GSMtX3ZZuMopDPPmvxtS0tBkXhoOBgfhAWVqqZ2MGDYiJ49KgGxmmAXk92\nHr2EVaY5858zZe8OdSTqCOvPrc/vjn5Lw4rl4k46ZQqpouvNbOLihIvnoYfIff/djM5nZnJQcDDr\ne3lxXnQ0062wsU9RFHqmpPClM2dY38uLv8bGMs9gKP+C/9Xhr/jeP++ZfcHMJjub7NWLHDfO4tnF\njRtijamktcKEnBw+dfo0XwgOZqrGX6KPdn1km93JtxMTI1Tz9GmbDJ+Wl8e2J0+a7TI7dEiYXdKC\naFkhNz+XY3aMYevFrfnG+Ag+8giZnGxrq+7m6uVQHu7iwItNajHF97itzfmP3bvJpk2ZMWoUp549\nS4fjx+l++bLNKq35paayT0AAW508Wb4FPyU7hQ6zHRhxLaJUF8pkcnPJwYPJYcNU8xmOHEnOn2/8\n+Km0NDbx9uZXkZE0WOHxNSEjgXVn1+XF6xc1H6tY1qwRTmYbhbaEZ2bSydOTJ1JSTGp/+vS9NbPP\nzM3ks38+y+fXPs9UfSoVRXgpOnYsY6J/7BjZpAnzx4/j5H8+YotFLbT/nZvBjitX2HjnTr4xaxZj\nT5ywtTlUFIW7k5PLt+DPOjaLI7aMKNUFMhmDgRwxgnz+eVWzU/37rwgwKErL/4qPp6OnJz3UjDww\ngW+OfsOhm030NWmFooib6xdf2MyEf5KS2Mjbm/El/P+OjBRr91u2WMkwC0nJTuGTq57km1vfZJ7h\nvyfGMif6ixeT9eoVqv281HcpXea60Dva24aGkel5eXw/LIxNT5wQ62lbtojQp59+sklo8Z2UW8HP\nys1ivTn1GJKgYWSHoohwhl69SJX954pCtmxJenr+97ccg4Hjw8PZ/MQJBtvAV5mek876c+vzVOwp\nq49diLg48SM6ZTs7voqMZK/Tp43GQCcmkg8/LLTpXiAxI5Gdfu3EcbvG0aDc/ZnKhOjn5oq45TZt\nRCTcHewK30Undyd6hHrYwDjyREoKm584wVHnzhV2sV6+LNLhvviiWBG3IeVW8Jf4LuGg9YNKfWFK\nRFHEIlHXrmRqqiZDzJsnor1IMlav5xP+/nwhONiizVSWstR3Kfut6Wez8W/x559k+/Y2y/mbrygc\nEBTETy5cuOtYdjbZvTs5bZoNDCsFCRkJbL24Nacfml7sgrSikJ9+KsLbrT7fuH6d7NtXPEkX83s7\nffU0G8xrwJX+K61mWkGUnLOnJ7cYe+rOySE/+IBs21bEvtqIcin4eYY8Nvu5Gb2ivSy6OMUyY4bw\nJWs43UlKEou32y5fZwMvL35/6ZJV/PXFkZufS9dfXLk/Yr9N7aCiiJ2TM2bYzIRrubl80Nubf9/2\nI1cUsXv21VfL1qYqY1zLusYOyzrwq8NfmdReUch33xVhm1ZbRomNFb+1iRMLpRU3RnhyOB9c8CDn\nexezCKYSmfn5HHn2LDv4+jLClKf8ZcvE0+nBg5rbVhTlUvDXn1nPp1Y9ZdGFKZZvvhGPlQnapmnI\nVxS2mxPF2vu9eOCamdvKNWRjyEZ2Xt65yEd/q3LlilgRDQy0mQknb+7ELfixz58v3B4ZGTYzyWRS\n9ansuqIrP9n7iVmhpvn55Msvi5fmgScXLpDNmoncD2bYGJ0SzRaLWvDrI19rlv4hJjubnfz8ODw0\nlJnmXIijR4Xor1mjiV3FUe4EX1EUdlzWkTvP77T44hTJ99+TrVqR8dru2o3T69k3IICPHA5gqyf0\nZWq2aFAMfHTFo1x/Zr2tTRHVsTp1stmGLJJcfOUKH/Hz4479+XRxIS9dspkpJpOZm8mnVj3FMTvG\nlEoQ9XqyXz+RYUCz72ZQkFj1Xr68VKfHp8ezw7IOnHpgquqiH5CWxkbe3nS/fLl0fZ89K7bUWzm/\nRrkT/L0X9rLd0nba3NW//55s0YK8elX9vm/jn6Qkunh5cUZkJPMMCl1dyTIQ2VWIQ5GH2Hxhc9vv\nvizYkPXjjzY0QeHAk6Gs9mXYPRF+mW/I56D1gzhiywiLntLS00UCuG+/VdG4AgICRCSOhTkokjKT\n2H5p+xLXJ8xh77VrdPL05CZLn/CvXhWTldGjrfCoJCh3gt/7995cE6jyo5KiiDDANm2KLyhqIRn5\n+Rx9M6TL87Y4b3d3ctQozYYtNf3/7M/FJ8tAGMqlS2KnZViYTYbX68lHHs+j8x4f/qliBkUtUBSF\nH+78kE//8bQqN+u4ODFR/f13FYwrIDBQiL2HOtE2iRmJbLe0HWcctny95/e4ONbz9Cz0+7SItDSx\nGD10qFUCEMqV4J+8cpJNFjRhbr6Kj/cGAzl+vJjKaJAmtwDPlBQ+7OPDkWfP3lVTNSFBJOa7cUOz\n4UvF6aun6TLXhek5NtzOXsCiReQTT9gk1vmjj0TemcC0dDp6evJcGXbg/3T8J3ZY1oEp2SoJFoV3\nwtnZvPxPRgkKEsV3N5mQFdMMEjIS2GZJG35ztPSpORbGxLCxt7f6/3+zs8lBg0QEkkbpUQooV4L/\n6qZXueDEAlUuDEmRmnfUKBFDq9Yd/Q4y8/P58YULrO/lZTykiyJFh40yBBfL8C3DLfoRqYbBIATf\nyoHvGzeK1CkFN+MVsbFs5+tr3iKelVgXvI6N5zfmldQrqvd97JhYPz9zxoJOzp4VYq9RKtH49Hi2\nWNTCbI1QFIXfRUXR1ceHl7SqAJ+bKzZw9uypWZg3WY4E/+L1i3SY7cA0vUobG9LSyOeeIwcM0Czk\n4sC1a3T18eHw0FAml7DoWNzOW1tScN2Npqu1JufOCdeOlVZNw8OFyPn7//c3RVE4PDSU79nIvWQM\nnxgfOro7Mjg+WLMx/vpLuHdKFc9w+bLIdKmqb6iIYVIu88EFD5pc40FRFE6JiGA7X1/GaR2HajAI\nf36PHppNMMuN4H+06yNOPTBVnasSGyuyXr7/vpjlq0ysXs/XQ0P5oLc3t5voJlIUsYRQFhcFx+8e\nzwm7J9jaDMGsWWIRV+M7Y06OyNZc1ANFWl4eW/j48C+NI7lMJSY1hg3nNeQ/Yf9oPtaMGWLTmVme\niYQEEQyxQMWn82IITw5ng3kNuDGk+CcJRVE4+cIFdvbzK3FCphqKInyE3bpp4sMtF4KfnJlMu5/s\neDVNheiZwEDywQc1CZfSGwycGx1Nh+PHOe3iRbMf+3/5xfS0ydakILFa5HXTCkBoSm6uuFlrHOM8\ndSr5wgvGvyKB6cKff16DnOfmkJmbyS7Lu1itnoGiiO/o66+b+PNJTRXrY19+qblttxMUH0TnOc7c\nc2FPkccVReGnERHs5OfHa9YO+VUUscns0UfFDmMVKReC/+3RbzlqmwphLH/+KVwC69WNL1cUhRsT\nEtjsxAm6BQWVetEnJUUs3paRiWMhvj7yNYdvGW5rMwT+/mIVUaMLdeSIqP9aUlTesitX2NHXl9k2\n8ucrisLXPV7n8C3DNa09eydZWWKWP3NmCQ1zcsSW3TFjbOKr9Ir2opO7E31ifAr9vcCN84gtxP4/\nI8hPPhE3QxVF/54X/IIkaaGJpasrSVJ88caPF5WqgtXzcSqKwn3XrrG7vz8f8fPjQRX+x733nnj4\nKGsUJFbzv+pfcmNrMGUK+dprqnd7/bpwNe8pemJYCEVR+EpICMeeP6+6HaYwx2sOuyzvYtX6swXE\nxYla30aDbRRF5AAfPNhqMehFsSt8F+vNqceziWdv/e2ryEh28PW1nhvHGIpCfvyxyNelkk//nhf8\n5aeW022tW+mvQHi4KJLt5qaaz0xRFO69KfStTp7k2vh41XLgnDolfkhlMAik7CRWI8U009WV3LZN\ntS4VRdxDJpixXHEjN5fNTpzgZo3TcNzJ4cjDrDenHi/duGTVcW/n9GnxwOxf1Bxg+nTxu7Oxy4sk\n1wSuYeP5jRmdEk33y5fZ6uRJJtgoKd9dKIooqtS9uyrRO/e04BsUA1ssasHDkYfN/+SKIlbcHByE\nc1yF+G29wcDf4+L4iJ8fW588yfXx8UUWEraUrl0LpQIvM+Tm57LFohaFa43akqNHyYYNVbuRr10r\nFs7NDZX2vZlvR6sSlHcSnRJNl7kuPHBRjcB4y/DwEE9EhTanL18ubsZWrudQHPO859Fl7Qd80NuL\nMVqFXpaWgjTsTzxhcZrSe1rwt4dtZ+flnc33T54/Tz7zjJhhqPC4HZmVxRmRkXTx8mL/wEDuTk7W\nNKvlH3+Q/ftr1r1FeIR6sOOyjrZPrFbAmDEivaOFxMaKEMzSpuBfEB3NrqdOWVQE3RT0eXp2W9mN\nPx3/SdNxzOHbb0XQSXY2yb17xS5aG6YILoo1cXGsdXg3u/wx0CYusBIxGMT3uE+fUm/OWr/+Hhf8\nnqt7cl3wOtM/cUoKOXmyeM6cO9eikMuUvDz+ERfH3gEBdPT05LjwcIZYaYelXi/WJMtYqDdJ4c7q\n/lt39dNblJbUVDHFtGAbqKKI7RglLkIW24cogj5JY6H7cOeHHLJxiFUXaUuiIHLn84EhVMpgcd/t\nSUms5+nJM+lpHLFlBAetH1So4leZIT9fbM569lmzc1Nv2CD2tFld8AE8ByAMwAUAnxtp88vN40EA\nOhlpw8bzG5uWRiEtTeStrVdPrHqWMnrjql7PlbGxHBAUxFrHjnFgcDA3JyRQb4Pt/NOnC9deWeTY\npWNssqAJs/PKyKPxzaLSpX0cXrlSnYSc13Jz2cTbm/9olKJjbfBauv7iqmraBLXIjEpgTNVm3DH0\nT1ubUohjN27Q0dOTJ2/6x3Pyc/jMmmf4/vb3y9RN8xZ5eeQrr4iYYBPXGTZtEmIfHGxlwQdQCUAE\ngKYAqgAIBND6jjbPA9h989+PAfAx0hfnec8r/pPGxoqAaQcHMcUICDDpAhUQp9dza2IiJ4SHs83J\nk7Q7fpyvhoRwQ0JC4RJmNiAmhrS313QXtkW8uOHFMuVW4MiRIrbZTKKixAOhRSkDbsMrJYXOnp68\nrLKf+GziWTq6OzIgzrzvuFXIziZ79GDKhC/ZoEHZWX8KTE+ns6cn999RayJNn8Yuy7vw6yNf28iy\nEsjNFYL/yisleim2bBHz3IKSEdYW/B4A9t72fiqAqXe0+RXA0NvehwGoV0RfTNUXoXZxcWIxtlcv\nUS5q3Lgi61/ejt5gYEhGBjclJPDrqCgOOXOGjb29aX/8OPsHBnLWpUv0TU3VZAHWEl57Taw3l0XO\nJ58vOykXSPLaNRE8f3uR4BJQFBEmrnbmZffLl9nd3181f35GTgbbLmlr1bJ+JqMowg3x6qukwUBv\nb3EDDdGw1LQpXMzKYgMvL240Ej0Vnx7P5gubc/mp0uXi1xy9Xrh2RowwGrK3fbtw/Z4+/d/fzBV8\nnTindOh0ulcAPEvy/Zvv3wDwGMnxt7XZAeBHkt433x+86frxv6MvXjl2DExNBaOioJw5A0NoKAzX\nryP36aeR8+yzyOneHZmVKiHdYECGwYCU/Hwk5eUhOS8Pibm5iMnJQXRODq7n5aFZtWpoXbMm2tSo\ngbY1a6JbrVpoXr06dDpdqT+v1nh6Au++C5w7B9x3n62tuZsJeyZAoYLFzy+2tSmCrVuBzz8HgoKA\n6tVLbL5iBfDbb4C3N1C5snpmKCQGh4SgRfXqmOfqanF/b297GwoV/PHiH2Xv+/rTT4CHB3DsGFCj\nBgBgzRrgm2+AkycBR0frm5SYm4snAgIwqVEjjG3Y0Gi7iOsR6Lm6J5a5LcPgVoOtaKGJZGUBbm5A\n8+biy3qbCOzbB4wcCezaBXTt+t8pOp0OJE3/kphzd7jzBeBlACtve/8GgEV3tNkB4Inb3h8E0LmI\nvlj9lVdY7ZVXeP8rr7DaV1+xyrp1rLtrF9ucPMlH/PzY7dQp9g0I4KDgYLbfsIGYOJF46y3ixReJ\n3r2JNm046Ycfipy5z5w5kwDues00snJnq/YuLleL3ARUFuxPykyio7sjw5LCyoQ9JMWe/8mTS2x/\n+XJhV47a9hTUw912059f6v4fAfERiCpl7/s5CGBqrVqiFOUdPP74cQJHCFSx6vfhi2+/5aOnTvHL\nyEiT2r8/8306uTvx+OXjNr+eRbZPTycff5wnu3a9rV0fAgkEevCtt97izJkzb71gZZdOdxR26XyB\nOxZuIVw6r9/23qhLRyIq/A0YYGsrjOPu6c5B6wfZ2oz/SEoSK1hexgvcFxTR+u47bU05kZJCJ09P\nRpYyzC40MZSO7o4MSbCxf6QogoLEHfPkySIP5+cLN/T771svq0KOwcD+gYF8LyzMrAXZfRH76DzH\nWdNMoxaRkiI253z8MY8fU+jkZDzJorUFvzKAixCLtlVR8qJtdxSzaCsR62EuLmSoBZkltCQ7L5tN\nf25aus1xWuHhQbZsaTSe2ZplcudHR7OLn5/Z+XYyczPZdklb/u/0/zSyzAISE0VU1LriQ6bT0sh2\n7ciff9beJIOi8I2zZ/lCcDDzSrF2si54HRvNb2TTncvFcv06M1p24sLqn3P/PuM3M6sKvhgPAwCc\nh4jW+eLm30YDGH1bm8U3jwehCHcOpeAX4ttvVdlbpBmbQzezw7IOzDeUoXwQQ4eK5FR3ULDBqiCq\nQWsK8u2MNnNTxTvb3uEbf79R9kIHc3LIp54ip00zqXlUlFhL37lTW7M+i4hgD39/iwrTLPRZyBaL\nWpSdQITbCAggWzomM7Vpe/Krr4y2s7rgq/WSgv8fSUkiRLMsZtEkhaj1Wt2Lv/r9amtT/iMpSSjN\nbc++iiJyeRXze9GE1Jv58/8wsR7uX0F/seWilmWjtOTtKIrw0QwebFaqkhMnxE02KEgbsxZER7PV\nyZOqJEObdnAauyzvUnSEoI0ICRFP+R4eFE9Xbdsa3SUoBb+cMGaM9YXKHALiAlhvTj1ez1I3v7dF\n7NhBNmsmfAsUuxHbtDF7E6MqnLmZPz+4hM1h55PP09HdkUHxGqmjJSxaJHw0aeZXnduwQSQFvKpC\nSYtC/SYksJG3t2qlCRVF4ZgdY9hrda8ykYIhLIxs0EBUG7tFQoIQ/Rkz7logkYJfTjh/XsySykDy\nQaOM3jGaH+/52NZmFObdd8n33ru1luvjU/IpWvFnXBxdfXx4w8hMNDsvmx2XdeQyv2VWtswEDhwQ\nO3zuiH4xh+++EzU/1MpQcvD6dTp7ejLIwoRjd2JQDBzmMYxua91M2+mvERERZKNGYs3pLhISxM33\nyy8Lib4U/HLEoEHksjKoBQUkZiTS0d2xUO5xm5OaSjZtyjm9dhTl0rc6E8LD+XxQUJHJ98buHMtX\nN71a9vz24eFih4+F9TcVhXz7bXLgQMsri/qnpdHJ05P/alAmkBSZYd3WunGYxzCbrE1FRYnCfL8W\n5yVNTBSiP3XqLdGXgl+O+Pdf8uGHVcn0rBkLTixg/z/7lynR8v7xKOMr1WdmlHVz1hdFrsHAXqdP\n3xUnvjl0Mx9a+FDZy5Nz44aIeFqxQpXucnNFJtjRo0sfrhmRlcX6Xl7conH65azcLPb5vQ9HbRtl\n1eywly4JT6RJu+yTkkTI2YQJpKJIwS9PKIpIQ+vhYWtLjJObn8u2S9pyc+hmW5tCUoQwN2pEXho2\nVRTBKQM3ooScHDbx9ubfNwUr4loEHd0d6XvF18aW3UFentiwYE5FGBNISxNliUtT2S0+J4fNT5zg\nsiI2e2lBRk4Gn1r1FD/Y/oFVJjGXLwuxNyuU9cYNskcP8r33pOCXN7ZtE2Uwy4BuGeXYpWNsNL8R\n0/TmL+6pzQcfiBdzcoQDedEiW5tEkvS7WTTFP+U6Oy/vzIU+C21t0t1MmCDqSmiQSPDqVeGyWL3a\n9HNu5Oayo68vv46KUt2e4kjTp7HHbz04btc4TUX/8mXyoYfIBQtKcXJ6OtmnjxT88obBINx2ptRc\ntSVvbn2Tn+771KY2HD4sZve3yoWGh6ubGtNC/oyLY+1DO/n8RusWITeJJUvIVq1UqyZWFGFhYiF9\n69aS22bm5/PJ06c5PjzcJtcqJTuFXVd05fjd4zVx71y6JMR+/nwLOsnKkoJfHlm7lnzySVtbUTzx\n6fF0cnfimQTbiGtGhvgB3ZWqd/Vqcce0UjnC4vAI9WCdDZPZ/ZSvTWouGGXfPhGRExGh+VD+/iL6\n7OBB421yDQa6BQVxRGioppXmSiIlO4U9fuvB97e/r6roX7ggnnYWqvCQJwW/HJKXRzZvLhZxyzJL\nfJfwqVVP2WRGNmkSOXx4EQcKyjONHm11m27nwrULdHJ3ok/MSb505gzfPneubMzyQ0OFAluxatW/\n/4ohi0rLk68oHB4aSregIOaWgZtiek46e//emyP/HqlK1axz58RT6HKVsjRLwS+nrFxZduveFpBv\nyGfXFV2tnsfdy0u4CowWnkpNJVu0uGM3i/XIys1ix2UdufjkYpJkRn4+H/Hz46xLl2xizy3i48WK\noTmOdZXYuVNEft5eU1hRFL4XFsbeAQHMsiBlgtpk5may/5/9+fLGly2q+hYUJDaD//GHerZJwS+n\n5OSImYFvGQvsuJOg+CA6ujsyNi3WKuNlZQkt31xSkFBBtkcbZKV79593+brH64Vm9Ff0ejbx9uYa\nE9MvqE5GhsjIOGOGbcanCEgoEH1FUTghPJzd/f2ZZuPqc0Whz9Pztc2vsdfqXryRbf46R8FTzcaN\n6tolBb8c88svItKwrPPV4a84eP1gq7gsPvtMFF8yif/9T+RasFJxepJcdXoVWy9uXWSenNCMjCJL\n8mlOfr7Y1ffmmzYP/9q2jXRyVjjK6yI7+fkZ3ZVcFjAoBn6852O2W9qOMakxJp+3bZsQ+wMH1LdJ\nCjrE/+gAABTmSURBVH45Rq8nGzcmvb1tbUnx6PP0bL24NTeFbNJ0HB8fsdZopKrd3SgK+dZbwtlv\nBaELjAuko7sjQxONP1Ucu3GDTp6eDChFvppSoSiiTGjfviYXzdbWHIWv7rnISmtOcqen7e0pCUVR\nOMdrDhvPb8zTV0+X2H7lSuFu9PPTxh4p+OWclSvJPn1sbUXJeEd702WuC5MzkzXpPztbTNbXrzfz\nxMxMsbHB3V0TuwpIzkxms5+bcf2Zkg3cnJDABl5eDLdG4qQffxRRSxqGX5qKoiicdvEi2/v6cv2e\nHDo5kbt22doq09gcupmO7o5cF1x0jQCDgfz8cxE5dv68dnZIwS/n5OaSrq7Fh7WVFSbumcihm4dq\n4tr57DPy5ZdLOVGPjharZxqpS54hj/3W9ONn+z8z+ZyVsbFsomIWyKIHWSkKmVhp12pxKIrCqTfF\nPvHmk8aJE+KJbc0aGxtnIoFxgWz2czNO2T+lUP6djAzypZdEGQGjgQQqIQW/ArB2LfnYYzZ3v5ZI\nZm4mWy9uzb+C1I2OOX5cPCZblFrF01M4Vs+dU82uAqbsn8J+a/qZHca3MCaGzU+c4BUt8jlv2SJu\ncuHh6vdtJgZF4ccXLrCDry+T7nArnT0r0irPmlX2v9+keJLr+0df9v69N2NSY3j5sniAfOst66Tl\nloJfASjYfbt9u60tKRn/q/50cndSrZRcerrYk2DKbs0S+e03kZ0uWT2304YzG9j056aldmX9dPky\nW508yXg1/euHDokIJX9/9fosJXkGA98+d449/P153cgCbUyMyIrx+utlOz14AfmGfH7373e0+96Z\ndXp4cM4c692szBX8+yC557jvPuC774Dp0wGDwdbWFE/n+p0xucdkvLXtLRgUy439/HPgiSeAF19U\nwbh33wUGDwYGDQKysy3uzi/WD+P2jMPWoVvhUMOhVH183qQJhjs7o2dAAGL0eottwvHjwNChwObN\nQOfOlvdnAXqDAa+ePYu4nBwc6NgR9lWqFNmuUSPg2DGgcmXgySeB6GgrG2ourIScA1+iypZ/UHPw\nFJx9+B3c0F+3tVVFY87dQcsX5AzfLBRFpFtQKYutpuQb8tlzdU/+dPwni/rZt0/sRVB1vdFgIEeM\nEGX8LIj/jk6JZsN5Dbnt3DZVzJoXHc2mJ07wgiVTXC8v7eIBzeRabi57BwTwtZAQ5pi4g1ZRyHnz\nhF9flSc6DYiIIHv2FIEU8fFkqj6VY3eOZb059bg6YLXmocmQLp2Kw6lTwpedUsZSqhfFpRuX6DzH\nmZ6XPUt1fmKiKP126JDKhpEiPLFfv1InbU/PSecjvz7COV5zVDVreWwsG3p58UxpKjydPCnEvgxk\n3QvPzGQLHx9OvnCB+aW4vl5eItrlvfeES68sYDCQixeTDg7ipnTnxmC/WD8+uuJRPrnqSfrEaFd2\nTQp+BWPUKBGxci+wK3wXG85ryLh083aXKorYJ6Tp50xNFYUlpk83S/TzDfkctH4Q3/3nXU1mc+vi\n4+lk7uas48eF2N+VSc76HLtxg/U8PflrrGU7r9PSxHe9eXPbR6gFB5O9epHdu4sMoMbIN+RzxakV\nbDy/MZ/76zmeiDmhui1S8CsYV6+SdeuKDHz3AjMOz2Cv1b3MimD59VcR+aD5PqGCEnJ31A01hqIo\nHL1jNPv+0Zc5+doZ9+9N0VxqSjjl7t1C7Pfv18weU1AUhYtiYujk6cl9Ku4k3rZNRJYOGWJRud1S\nkZhIjhkjLu8vv9w9qzeGPk/PX/1+ZZMFTdhrdS+uCVzDzFx1VqOl4FdAZs0iX3zR1laYRr4hn8/+\n+azJufPPnhUBJhpETxaNGaI/4/AMdl7e2SqFXy5kZrKljw8nhIcbzyK5YYNITmPjrdhpeXkcGhLC\njr6+lq1BGCErSxRIr1uXnDJFTHq0JCGB/Oor8T38+GPy+vXS9ZOTn8PNoZs54K8BtP/Jnu/+8y49\nQj1KFdGVZ8ij/1V/swVfRyG2Nken07Gs2HKvodcDrVsDK1YAzzxja2tK5lrWNTy68lF83+d7jOgw\nwmi77Gyge3fgo4+ADz6wooFJSUDfvsALLwA//ADodHc1Weq3FAt8FsDrHS8413S2ilk38vIw4tw5\n3MjPx7rWrdGsenVxgAQWLgTmzAH27AE6dLCKPUURlJGBoaGheMrODr+4uqJ6pUqajXXlCjB7NrB2\nLfDSS8Cnn4rfgVqcPSsu66ZNItBp8mTg4YfV6Ts2LRYbQzfiYORBeEZ7onnd5ujk0gkP2T+E5vbN\n0aBWA1S+rzIq3SeuX3xGPKJToxGTGoOghCD4XPFB4zqNcfajsyB59xfUCFLwywl79ghhPHMGqFnT\n1taUTEhiCPr+0RcbX9mIPs36FNnmvfeArCzxgy5Cc7UlKQlwcwPatBF30qpVbx1af2Y9PjvwGY6P\nOo5m9s2sapZCYsGVK5gdHY0lDz+MV+3sgHHjAB8fYMcO4MEHrWpPAbmKgh+jo7EkNhbzmjfHSBcX\nq4197RqwdCmwZIkI6Xz5ZfFq0cK8fkggLAzw8BBRrMnJ4js4bhzgrOE9Pc+QB7+rfghNDEXkjUhE\npkTiavpVGBQDFCogiHo166Fx7cZoUqcJWju1xhONn4BDDQfodDop+BWVESOA+vWBuXNtbYlpHIk6\ngqEeQ3H4rcNo59yu0LE//gB+/BHw8wNq1bKRgZmZ4qKmpwNbtgB2dlh3Zh0+3f8p9o/cf5fN1uRU\nWhpeDwnBY6dOYZ63N1xWrrTZhQpIT8eosDA0vP9+LG/RAo2qVbOJHfn5YtuBhwfw999AlSpAp05i\n+0GbNoCdHVCnjrhMGRniRpGcDFy4IL5nfn5AtWriaeG114DHHxd7XsoyUvArMElJQPv2wM6dwKOP\n2toa01gbvBbTDk/DiXdPoEGtBgCAkBCgTx/gyBGgne00VWAwAJ98Ahw8iO0/vYMx4fNwYOQBtHVu\na1u7PD2R+c47+Pbzz7GqZUvMbNoUHzZsiEpWfBSKy8nBjEuX8E9yMuY0b44369WDzuqPYkWjKEBU\nFHD6NBAQAJw/D6SkAKmp4v79wAOAoyPg4AA0awZ07SpeDRva2nLzMFfwbb5YW/CCXLRVhT//JDt0\nEEnW7hV+PP4jWy9uzbj0OKamki1bkr//bmurCnNi+ttMrnkfryz+0baG5OWRX38tdiPdzK0RkpHB\nnqdPs4OvLz0SEzWvA3sjN5dfR0Wx7vHjnBIRUaZz2Jd3IBdtKzYkMGCA2JL+5Ze2tsZ0vvv3O6w9\nsxYNDx5Ci/oNsWyZrS0SkMQc7zlY5LsIRx9ZiOZjp4up4JIl1nehnDsHvP++8DusWQM0aFDIzu3X\nruGHy5eRaTBg2oMP4hUnJ9yvok8iPCsLv1y5gnWJiRjo4IBvmzZF04KFY4lNkC4dCaKjhSZt2wb0\n6GFra0yn57TZOGX4DSFTDuMhh8a2Ngf5Sj4m7JkAz2hP7P5/e/ceHFV9BXD8e5pAKiSI4RGCpCJa\nHmIRgRJG6oAIrWMVGZEqvmiZqWhtDXWmCGVU8IEMHdROMk4By0MNSovW1yAQKVEcnCAaHoIZTC0G\ngklMwSQbmgSS0z9+twzoJpDsJjfZez4zmbl3c/fuucnO2d/+7u93fndsoF+3fq5fPyMDNm92N0um\nTWv9O8pVVa540qpV8Mgj7u58I4lcVck5dozFRUXsDoWY2quXq83TvTvfa0GcRTU1vFFezj/Ky/m0\nupp7UlP5zYUX0jchIdKrMlFgCd8ALtnPnu36Ly+4wO9ozm71anjiCZjxl6dZtS+Lt6a/5Ws/eVVt\nFbe/djs1J2tYP20953///DMPeO89eOABSE52Y/daYyhkbS28+CIsXOiGiS5ZAikp5/z0opoaXikr\nI7u0lKLaWtKTkhjTrRujkpJITUigd6dO9OrUiQagqr6eqpMnOVJXx+5QiN2hEDurqjhcW8sNPXow\npWdPfpac3KrDLE3zWcI3p2RkwKFDboBJO7mXFtYHH8DNN7scOmQIvLTnJR7c9CDLb1zOlMHRKIvZ\nPPlf5XPr+lu5pv81ZF2fRae48FUdOXnSDdlcsADS090n7IQJkf+xKyth2TJ49ln3QfLww27ISARK\n6+rIq6zkw8pKdoVClNbVUVZXR9mJE8SJkBQXR1JcHCmdOzOsa1euSExkeGIiP05KIr69D1UJMLtp\na06pqXElCbKy/I6kcbt2ucmhmzad+fiOwzs07ek0fXTro1rfcG7VFSPV0NCgmXmZ2mtJr3NamvCU\n6mpXtnToUDdL96mnVPPzm1eIraJCde1at4xXt26q06e7P44xTcBu2prTFRa6xuGrr8LVV/sdzZkK\nC2HcOHjmGTfu+dtKQiVM+/s04iSOZTcsY1DPQa0Xy9FCMjZmUBIqYd0t67g0+dLmn0TVjSV9/XU3\nEy4Uchc4YIAb+5eW5lr/tbXu59Ah2LPHzZb7/HP3D5o61dXo79GyevomWKxLx3xHTg7cdRfk5sLg\nwX5H4xQXu/w2d27TZRPqG+rJ2pHF4+8/TkZ6Bg/95CE6x3Vu/AnNVF1XzaJti1j28TLmjJ3D7DGz\no3f+wkI3A/bgQTcovKjIJfyEBPeTmuq6bIYNcxMOOsIUadOuWMI3Ya1eDY89Btu3QxvOeg+ruNjV\n/Ln7bpfwz0VRRRH3b7ifgvICt4LWFTM4r1PLhwQe/e9Rnv/keTJ3ZDLuonEsmbTk1MQvYzoKS/im\nUQsXunIrublupqEfCgtdsr/vPpgzp3nPVVXe//J9ln64lLziPO4deS+3XHYLl/e+/JxmeNbV15F3\nOI/svdms27eOGwfeSEZ6BiP7jmzh1RjjL0v4plGqMGuWqwL49tuutkhb2rPHTQpbsMDNH4pEQXkB\nWTuyeKfwHY6fOM7EARMZ0WcEKYkp9O7am24J3Sg/Xk5pqJQjVUfYfng7277cxsAeA5k8aDL3jLyH\nPok+f9UxJkKW8E2TGhpcaZitW2HjRteN3BY2b3b3ETIzw9+gjcQXx74g51857P96P2XHyygNlVJZ\nW0nPLj1JSUyhT9c+jOo7igkXT2jx4uLGtEeW8M1ZqcKiRbBypUvEl1zSeq/V0OAmiS5fDmvXukEr\nxpjoaG7Cj2/NYEz7JALz57tqgWPHurIwU6dG/3XKy+HOO91CJjt3tt23CWNMeDaFLsBmzYI33oB5\n89yImYqK6Jy3vt5NFB06FIYPhy1bLNkb0x5Ywg+49HRXb6drVzccfPVqqKtr+fm2bXOF27KzYdMm\nWLwY4u17pDHtgvXhm1Nyc10BswMH3PqdM2eeWwXgigrXP79ihVtk4skn4bbb2nf9HmNigd20NRH7\n6CO3OPTGja5bZtw4t5h4UpJb2jU+3k0ezc+HXbtc//ykSW7G7LXXtv9l4YyJFZbwTdTU1EBenqti\n+fHH7uZrba3r8klLc+uFDh8Oo0db6Rdj/GAJ3xhjAqK5Cd++fBtjTEBYwjfGmICwhG+MMQFhCd8Y\nYwLCEr4xxgSEJXxjjAmIFk96F5FkYB1wEXAQ+IWqfhPmuINAJVAPnFDV0S19TWOMMS0XSQt/LpCj\nqgOBLd5+OAqMV9UrLdkbY4x/Ikn4k4E13vYaYEoTx1pVFWOM8VkkCT9FVUu97VIgpZHjFHhXRHaK\nSIQL2xljjGmpJvvwRSQHCLfw5/zTd1RVRaSxughjVfUrEekF5IhIgapuC3fgggULTm2PHz+e8ePH\nNxWeMcYESm5uLrm5uS1+fotr6YhIAa5vvkREUoGtqjr4LM95FAip6tIwv7NaOsYY0wxtWUvnTWCG\ntz0DeD1MMF1EJMnb7gr8FNgbwWsaY4xpoUha+MnA34AfcNqwTBHpC6xQ1Z+LyADgNe8p8UC2qj7V\nyPmshW+MMc1g5ZGNMSYgrDyyMcaYsCzhG2NMQFjCN8aYgLCEb4wxAWEJvw1EMlGiI7Dr69js+oLD\nEn4biPU3nF1fx2bXFxyW8I0xJiAs4RtjTEC0q4lXfsdgjDEdTYecaWuMMaZ1WZeOMcYEhCV8Y4wJ\nCN8TvohcJyIFIvK5iDzkdzzRJCJpIrJVRPaJyKci8oDfMbUGEYkTkXwRecvvWKJJRLqLyHoR+UxE\n9ovIGL9jiiYRmee9N/eKyFoRSfA7pkiIyEoRKRWRvac9liwiOSJyQEQ2i0h3P2OMRCPX9yfv/blb\nRF4TkfObOoevCV9E4oAs4DrgMmC6iAzxM6YoOwH8XlWHAmOA+2Ps+v4vA9iPW84ylvwZ2KCqQ4Bh\nwGc+xxM1ItIf+DUwQlV/BMQBt/kZUxSswuWS080FclR1ILDF2++owl3fZmCoql4BHADmNXUCv1v4\no4FCVT2oqieAV4CbfI4palS1RFV3edshXMLo629U0SUi/YDrgeeJocXqvZbS1aq6EkBVT6pqhc9h\nRVMlrkHSRUTigS5Asb8hRcZbOvXYtx6eDKzxttcAU9o0qCgKd32qmqOqDd5uHtCvqXP4nfAvBA6d\ntn/YeyzmeC2qK3H/lFjyDPAHoOFsB3YwFwNfi8gqEflERFaISBe/g4oWVT0KLAWKgCPAN6r6rr9R\ntYoUVS31tkuBFD+DaWUzgQ1NHeB3wo+1LoCwRCQRWA9keC39mCAiNwBlqppPDLXuPfHACOA5VR0B\nVNOxuwPOICKXALOB/rhvnYkicoevQbUyb4WlmMw5IjIfqFPVtU0d53fCLwbSTttPw7XyY4aIdAJe\nBV5S1e+s+9vBXQVMFpF/Ay8DE0TkBZ9jipbDwGFV/cjbX4/7AIgVo4DtqvofVT2JW4r0Kp9jag2l\nItIHQERSgTKf44k6Efklrlv1rB/Yfif8ncAPRaS/iHQGbsUtjh4TRESAvwL7VfVZv+OJNlX9o6qm\nqerFuBt+/1TVu/2OKxpUtQQ4JCIDvYcmAvt8DCnaCoAxInKe9z6diLvxHmveBGZ42zOAmGp0ich1\nuC7Vm1S15mzH+5rwvZbFb4FNuDfbOlWNmZEQwFjgTuAab9hivvcPilWx9nX5d0C2iOzGjdJZ5HM8\nUaOqu4EXcI2uPd7Dy/2LKHIi8jKwHRgkIodE5FfAYmCSiBwAJnj7HVKY65sJZAKJQI6XX55r8hxW\nWsEYY4LB7y4dY4wxbcQSvjHGBIQlfGOMCQhL+MYYExCW8I0xJiAs4RtjTEBYwjfGmICwhG+MMQHx\nP9YqKxhgYoJdAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 57 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and check for their zeros..." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# zeros of Bessel functions\n", "n = 0 # order\n", "m = 4 # number of roots to compute\n", "jn_zeros(n, m)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 58, "text": [ "array([ 2.40482556, 5.52007811, 8.65372791, 11.79153444])" ] } ], "prompt_number": 58 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integration\n", "\n", "The numerical evaluation of the definite integral of a function is called numerical quadrature, or simply quadature. \n", "SciPy provides a series of functions for different kind of quadrature, for example the quad, dblquad and tplquad for single, double and triple integrals, respectively." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.integrate import quad" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 59 }, { "cell_type": "code", "collapsed": false, "input": [ "xl = 0 # the lower limit of x\n", "xu = 10 # the upper limit of x\n", "\n", "val, abserr = quad(lambda x: jn(0,x), xl, xu)\n", "\n", "print \"integral value =\", val, \", absolute error =\", abserr " ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "integral value = 1.06701130396 , absolute error = 7.43478946065e-14\n" ] } ], "prompt_number": 60 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can pass extra arguments to integrand function, by using the args keyword argument.\n", "For example we can use args=(3,) for the third order Bessel function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Exercise 11\n", "\n", "Create a function that returns the Bessel function of first kinf for different orders and then find the integral." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# type your solution here" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 61 }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Solution" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#%load solutions/bessel_integral.py" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 62 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interp1d function accepts given arrays describing points X and Y data.\n", "It returns a function that can be evaluated everywhere in the domain for an arbitrary value of $x \\in X$ and the image it is corresponding interpolated value. Options for interpoland kinds are linear, quadratic and cubic spline interpolation." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.interpolate import interp1d" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 63 }, { "cell_type": "code", "collapsed": false, "input": [ "x = np.linspace(0, 20, 20)\n", "y = lambda x: jn(5,x) # sampled on given values x\n", "f1 = interp1d(x, y(x)) # linear interpolation\n", "f2 = interp1d(x, y(x), kind='quadratic')\n", "f3 = interp1d(x, y(x), kind='cubic')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 64 }, { "cell_type": "code", "collapsed": false, "input": [ "# new points where we want to evaluate the interpoland\n", "xnew = np.linspace(0, 20, 100)\n", "\n", "fig, ax = plt.subplots(figsize=(10,4))\n", "plt.plot(x,y(x),'o',xnew,f1(xnew),'-', xnew, f2(xnew),'--',xnew, f3(xnew),'--',xnew,y(xnew),'k')\n", "plt.legend(['data', 'linear', 'quadratic','cubic','true'], loc='best')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAEACAYAAABmuKihAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdcVuX/x/HXYQgoiiAORIZ7pSXuCU5cuc2RW3Ml0rfS\nVByYmqVmirtylJblIBeJGoo4UDAtzS0uFLeIILKv3x+aP3cg49zA5/l43I+4ua9znfftCe4P51zn\nujSlFEIIIYQQIvMY6R1ACCGEECKnk4JLCCGEECKTScElhBBCCJHJpOASQgghhMhkUnAJIYQQQmQy\nKbiEEEIIITJZugsuTdNaapp2StO0s5qmffaadjU1TUvSNK1TevcphBBCCJGdpKvg0jTNGJgPtAQq\nAT00Tav4inZfAf6Alp59CiGEEEJkN+k9w1ULOKeUuqiUSgR+Adq/pJ0HsA64lc79CSGEEEJkO+kt\nuOyB8KeeX3n8vSc0TbPnURG26PG3ZGp7IYQQQuQq6S24UlM8zQHGqEdrCGnIJUUhhBBC5DIm6dz+\nKuDw1HMHHp3lelp14BdN0wBsgVaapiUqpTY93UjTNDnzJYQQQohsQymV6pNIWnoWr9Y0zQQ4DTQF\nIoAQoIdS6uQr2i8HNiulfF/ympKFtEVqeHt74+3trXcMkU3I/y8iteT/FZEWmqalqeBK1xkupVSS\npmkjgG2AMbBUKXVS07Qhj19fkp7+hRBCCCFygvReUkQptRXY+tz3XlpoKaX6p3d/QgghhBDZjcw0\nL7IdNzc3vSOIbET+fxGpJf+viMyUrjFcGUnGcInsyM8vCB+f7cTHm2BmlsTIkS1o06aR3rGEEEJk\nsiwdwyVEbubnF4Sn5zbCwqY9+V5YmBdAjii6Ht9ZLAyE/EEqRPYmZ7iEeEMtWnixY8dQ7AoFY2Vz\nHBPTK2B0C2USQdmitgzu148WXbtibGysd9Q38vivN71jCORYCGGI0nqGSwouIdLo/v37/PTTT4wa\nNZkHD1IoYG2BKlIQU7MCGJvYYJk/iTyWGg9PneL+nTu4urrSpEkTmjVrRqVKlfSOn2ryIW845FgI\nYXjSWnDJoHkhUilg714GDx6Mk5MTf/zxBxaNLWHCHe53sifapi93L63i1uGNlDN34dSWLVw6d46T\nJ0/SrVs3jh07RvPmzenerRtXf/1V77cihBAii0nBJcR/8N++Havy5WndoQNOTk4cPHKQBx0f4NSk\nEM7rRsLuiVD4BAypjtmI4th1ucXNBzcBKFasGD169OC7777j7NmzFC1UCKcPPmBokyYkX7ig8zsT\nQgiRVaTgEuIVkhMT6d6vH627dsW9Zk2iz5+n94jedNnahZIFS3LAcx/zZ3bEvfQ+XKNK0PzoYCbU\n+xSjwolUW1KNnRd2PtNf3rx5mbtwIYu3bOGHW7co3LAhWz/9FBITdXqHOUe/fv2YMGGC3jGEEOKV\nZAyXEC9xZPduWgwbxsOEBNbNmUPLtm35M+JP2v3Sjk/qfsL/6vzvtXfx7QjbQZ8NfRhWYxheDb0w\nNnp24HxsUhI9pk1j85w5uJcqhV9AAEYFC2b220qT7DRuqH///jg4OPD555+/tp2bmxu9e/dm4MCB\nWZQsY2SnYyFEbiFjuIRIp9WrV9OkXTsqV6nCzWPHaNm2LRtPbaTlTy2Z32o+H9f9+D+nTGheujl/\nDv6TgAsBtPyp5ZNLjP/Ka2LCxkmT2Pf33+yNi6Pfxx+TkpKSmW8rw/j5BeHuPh43N2/c3cfj5xek\nSx/PS01BIlNdCCF0o5QyiMejKELoJykpSQ0bNkxVqFBBHTp0SCmlVEpKivom+BtV/OviKuRKSJr7\nTExOVF4BXsr+a3sVeCHwpW2ioqJUw4YNVf/+/VVSUlK63kNGetnP5JYtu1Xp0uMUqCeP0qXHqS1b\ndqe634zo4/Dhw6patWoqf/78qlu3bqp79+5q/PjxKjIyUrVp00YVLlxYWVtbq7Zt26orV64opZQa\nN26cMjY2Vubm5srS0lJ5eHgopZQaOXKkcnBwUAUKFFDVq1dXe/bsSXWOrCK/H4UwPI9/LlNf56Sl\ncWY+5BeK0FNCQoLq0aOHcnNzU/fv31dKPSqWRviNUJUXVFYXIy+mq3//s/6q6MyiauruqSo5JfmF\n12NiYpSbm5vq3bu3wRRdL/uZbNHC65lC6d+Hu/v4VPeb3j7i4+OVo6OjmjNnjkpKSlLr1q1Tpqam\nasKECerOnTvK19dXPXz4UEVHR6uuXbuqDh06PNnWzc1NLV269Jn+Vq1ape7evauSk5PV119/rYoV\nK6bi4+NT/X6ygvx+FMLwpLXgkkuKIteLjYykU6dOREdH8/vvv5M/f36i46Np/0t7Tt85zb4B+3Aq\n6JSufbiXcefQ4EP4h/nT+qfW3Hpw65nX8+XLh5+fHxEREfTu3ZukpKR07S+zxMe/fHGKuLjUT+6a\n3j4OHDhAUlISnp6eGBsb07lzZ2rWrAmAjY0NHTt2xNzcHEtLS8aNG8fu3buf2V49d+nx/fffx9ra\nGiMjIz7++GPi4+M5ffp0qt+PEEKkhhRcIle7sXs3jrVqEQP4+vpiYWHB1ftXabSiEcUti+PX0w8r\nc6sM2VeJAiXY1XcX1YpVw+VbF/Zc2vPM63nz5mXz5s3cuXmT5pUqEX/1aobsNyOZmb28EDQ3T86y\nPiIiIrC3t3/me05Ojwrihw8fMmTIEJydnbGyssLV1ZWoqKhniqznx3HNmjWLSpUqUbBgQaytrYmK\niuL27dupfj9CCJEaUnCJXOviunWU7dcP67Jl2erri6mpKX9d/4s6S+vQvXJ3vn33W0yNTTN0nyZG\nJkxvNp0lbZfQZW0Xvtz7JSnq/wfLW1hYsH7TJo7ky0fDzp3BwAbSjxzZgtKlvZ75XunS4/DwaJ5l\nfdjZ2XH1uWL00qVLKKWYNWsWZ86cISQkhKioKHbv3v30sIUXiq09e/Ywc+ZM1q5dy71794iMjMTK\nykruCBRCZDhZvFrkShfWraPyRx9RrkED/vz5Z4yNjPA740e/jf1Y2HohXSt3zdT9ty7bmtAPQum+\nrjt7Lu/hhw4/YJvXFgDLvHk56OdHZRcXJg4dyufffpupWdLi30W5582bQFycMebmyXh4tEzTYt3p\n7aNevXqYmJjg4+PDsGHD2Lx5M6GhoTRp0oSYmBgsLCywsrLi7t27TJ48+ZltixYtSlhY2JPn0dHR\nmJiYYGtrS0JCAl9++SX3799P9XsRQohUS8uAr8x8IINCRRaJOnBAWZUsqd7u1k0lJz8awD7/4HxV\nbFYxtf/y/izNkpCUoD7d9qlymO2g9l7a+8xrC9auVUYFCqj9a9ZkaaZ/GfLP5KFDh164S3HChAkq\nIiJCubm5KUtLS1W+fHm1ZMkSZWRk9OQ4BwcHq3Llyilra2vl6empkpOT1YABA1SBAgWUnZ2dmjFj\nhipZsqQKCAjQ+R0+y5CPhRC5FWkcNC8Tn4pcJSUlhbZNm3IqJYXTO3dipMGn2z/FP8wfv55+lLIu\npUuuzac3M2jzIEbVG8UndT95cumr1ZAh7Nu2jeuhoeQtXDhLM8lkm4ZDjoUQhietE59KwSVylc8+\n+4x9+/bxxx9/kGyUTE/fnkTHR7P+vfVYW1jrmu3SvUt0W9eNIvmKsKLDCmwsbEhKSaFMw4Y0c3Hh\n+3nzsjSPfMgbDjkWQhgemWleiFdYvHgxv/32Gxs3biQyMRLXFa7YWNjg38tf92ILwKmgE0H9gyhr\nUxaXJS4cvHIQEyMjDm/ezI5Nm9i4caPeEYUQQrwhOcMlcgU/Pz8GDRrE3r17ibWMpe3qtnzg8gFe\nDb0McrmXjac2MnjLYMbUH8NHdT4iODiYDh06EBoa+mQKhMwmZ1UMhxwLIQyPXFIU4jm7N22iy8CB\nbN68mSjbKHr/1pu5LefSo0oPvaO91oXIC3Rb143i+YuzvP1yvpv3HRs2bGDPnj0YG6d+otE3JR/y\nhkOOhRCGRy4pCvGUMxs20HTAAMbNnMnfJn/Td0NffLv5GnyxBVDSuiR7B+zFycoJl29dcO3hiqmp\nKd8uWqR3NCGEEGmU7jNcmqa1BOYAxsD3Sqmvnnu9PfA5kPL4MUoptfMl/cgZLpGhEm/epHijRji/\n8w6NBzqy8fRG/Hr6UcamjN7R0mz9ifUM8xvGBza9+cZjEY3r9udBcmHMzJIYObJFmubBSi05q2I4\n5FgIYXiy9JKipmnGwGmgGXAVCAV6KKVOPtUmn1LqweOvqwC/KaVe+MSTgktkKKVo2qwZodev03RC\nJe48vMFv3X6jUN5Ceid7Y2F3w6i/pAHGvo4Ql4eIfx4tDVS6tBdz57pneNElH/KGQ46FEIYnqy8p\n1gLOKaUuKqUSgV+A9k83+LfYeswSkEXKRKabMWYMgUeOUGFgEcxMjNnRe0e2LrYAStuUxv4vN+65\nRHHtyj+UclgGQFjYNObN26FzOiGEEK+T3qV97IHwp55fAWo/30jTtA7AdMAOaJHOfQrxAj+/IHx8\nthMfb4JR8m32HltFuZ4VaVC5GrNazMJIyxnDFfPfKE9sMpSqEsW18GlAb8CUuLjMH0RvKJydnVm6\ndCkhISGcP3+e7777Tu9IQgjxn9JbcKXqHLdSagOwQdO0hsBKoPzL2nl7ez/52s3NDTc3t3TGE7mB\nn18Qnp7bCAubBiSCaW2Ma8dTv7oLs91n6x0vQ5mZJcG2bzg/5C3Mw60o85Yn5/5ZiLl5st7Rsoym\naWiaxtixY/WOIoTIRQIDAwkMDHzj7dNbcF0FHJ567sCjs1wvpZTao2maiaZphZRSd55//emCS4jU\n8vHZ/rjYAgr0gkLHSY5cxpW1p2Cgvtky2siRLQgLm0dY4DQKuazi/NY1ODub4+HRQe9ouUZKSgpG\nRjnjjKkQIvWePxE0efLkNG2f3t8ah4CymqY5a5qWB+gGbHq6gaZppbXHM0tqmuYC8LJiS4g3FR//\n+O+Gop8D6yD5Fzj+fo68zNamTSPmznWnhW040U4XKWCfl6pVwzPlLkVDppTC29ub3r17A3Dx4kWM\njIz48ccfcXJyonDhwnzxxRfPtP/yyy8pU6YMtra2dOvWjcjIyCevd+3aFTs7OwoWLIirqysnTpx4\n8lq/fv0YNmwYrVu3xtLSMl1/4Qohcq90FVxKqSRgBLANOAH8qpQ6qWnaEE3Thjxu1hk4pmnaEWAu\n0D09+xTieWZmSWB2CqIng+XXcLkjQI69zNamTSO2+U9l35itGL8fy779uzh27JjesbLcy1YI2Ldv\nH2fOnCEgIIDPP/+c06dPA+Dj48OmTZsICgri2rVrWFtb8+GHHz7Zrk2bNpw7d45bt27h4uLC+++/\n/0y/q1evZsKECcTExFC/fv3MfWNCiBxJZpoX2d7mLYH0GtmL+0ZFIexPAEqXHsfcuS1z/JmfcQHj\n2L56O4XCC7Ft27YM7fu/piLQJmfMkkhqUtp+7kuWLMn333/P3r17OXfuHCtXruTixYuUKlWKK1eu\nULx4cQBq167NJ598wnvvvUfFihVZsGABTZo0AeDatWs4OTkRFxf3wuXBe/fuYWNjQ1RUFPnz56df\nv34ArFixIt3v9U3JtBBCGJ60TguR3jFcQujOZ+N3JKYY07JMfR6W8MbcPBkPj5xfbAGMbzSeX/7+\nhXD/cPbv30+9evWybN9pLZSyQrFixZ58nTdvXmJiYgC4dOkSHTt2fKa4MjEx4caNGxQpUgQvLy/W\nrVvHrVu3nrS5ffs2+fPnR9M0SpQokbVvRAiR40jBJbK1XwK2ErBmC/MneTL848/1jpPl8prmZU6b\nOYw6MJzRgwax96mxR+L/OTo6snz5curWrfvCaytXrmTTpk0EBATg5OT05AyXnFESQmQkudVGZFsP\n4uPpN/hDmtSpkiuLrX+1LdcWoxrmHLh9m++//VbvOFkmLQXR0KFDGTduHJcvXwbg1q1bbNr06P6e\nmJgYzMzMsLGx4cGDB4wbN+6N9yOEEK8iBZfItuq83xVLU1M2/LxO7yi6MtKMGNJoBDUavsVYHx+9\n42SJf+fienrg/MsG0f/L09OTdu3a0aJFCwoUKEDdunUJCQkBoE+fPjg5OWFvb89bb71F3bp1X+j3\ndX0LIURqyKB5kS35BfjRqUM3fpk1no5DxugdR3f34u5RdkZpIucbs2TKFAYOGfLfG/0HGahtOORY\nCGF4snTx6owkBZdIrejoaOzK2tF8SDN+m7xB7zgGY+iWoZxcFc7JYxe5efx4uvuTD3nDIcdCCMOT\n1YtXC5HlugzugomTCasn/KJ3FIMyotYIzlc9yt1bt1i2erXecYQQQjxFCi6Rraz/fT1/bP2D9cvW\nY25irnccg/JWkbcoY1+Gth90ZNXy5XrHEUII8RQpuES2ER0dTZ/+fegyqgtNKzfVO45B8qjlwU2n\nY5w/c4b9+/frHUcIIcRjMoZLZBuVWjXmzq1/uHzgCmYmZnrHMUhJKUmUmluKvsl9CdkRkq7Z52Xc\nkOGQYyGE4ZExXCJHGrd8CWdCjvLrRxOl2HoNEyMThtUYRrhzOKdPn+bAgQN6RxJCCIGc4RLZQHhk\nJM4VK/BB7Yos3hiodxyDd+vBLcrNL8coo1H8FfoXa9aseaN+5KyK4ZBjIYThkTNcIsep36MzzkVs\nmbdqo95RsoXC+QrTvnx7kt9Jxn/LFg5t3653JCGEyPWk4BIGbcyyhUSE/sUP40Zimt9K7zjZxoha\nI1h6fCnOrq54zpihd5xsy9nZmYCAgAzrb8+ePVSoUCHD+hNCZB9ScAmDFRMTw7xRYxjk5kKD7umf\nOT03qVG8BsUsi/Fe/xYEh4RwOyxM70jZUnqX9TEyMuL8+fNPnjds2JBTp05lRDQhRDYjBZcwWO0/\naI9F2TzMX+uvd5RsyaOWB4FxfhSuUIH/eXnpHSfHSU5OTlU7GXslhAApuISBWrt1Lbv8drHph02Y\nGJnoHSdb6lKpC//c/If+fd9jTWAgydHRekfKMEeOHMHFxYUCBQrQvXt3unfvzoQJE1ixYgUNGzZ8\npu3TZ5n8/PyoVq0aVlZWODo6Mnny5Gfarly5EicnJ2xtbfniiy+eec3b25suXbrQu3dvrKys+OGH\nHwgNDaVu3bpYW1tTvHhxPDw8SExMBKBRo0YAvP322+TPn5+1a9cSGBiIg4PDkz7Dw8Pp1KkTRYoU\nwdbWFg8Pjwz/txJCGAYpuITBiX0YS98Bfen1WS/qla+nd5xsy8zEjMHVBxPldB4jS0tmLF6sd6QM\nkZCQQIcOHejbty+RkZF07doVX1/fVF3+s7S0ZNWqVURFReHn58eiRYvYuPHRzRgnTpxg+PDh/PTT\nT0RERHDnzh2uXLnyzPabNm2ia9euREVF0bNnT4yNjZk7dy537twhODiYgIAAFi5cCEBQUBAAR48e\nJTo6mq5duz7TV3JyMm3btqVkyZJcunSJq1ev0r1794z6ZxJCGBgpuITBaT+0PRb2FiwfK8vTpNeQ\n6kP45fhqxo76Hzv8M/jSrLc3aNqLD2/v1Ld/VdvXOHDgAElJSXh6emJsbEznzp2pWbNmqi7dubq6\nUrlyZQCqVKlC9+7d2b17NwDr1q3j3XffpUGDBuTJk4cpU6ZgZPTsr8h69erRrl07AMzNzXFxcaFW\nrVoYGRnh5OTE4MGDn/T3X0JCQrh27RozZ87EwsICMzMz6tevn5Z/CiFENiIFlzAoX/38PbvXB7Hp\nx00YGxnrHSfbsy9gT4vSLbCsHM+pkyc5duxYxnXu7Q1Kvfh4XcGV2ravERERgb29/TPfc3JyStW2\nBw8epHHjxhQpUoSCBQuyZMkS7ty586TfEiVKPGmbN29eChUq9Mz2T78OcObMGdq2bYudnR1WVlZ4\neXk96e+/hIeH4+Tk9EJRJ4TImeQnXRiMuzHRTBg7hY4tXKlfSf7Szygjao5gyV9LGD58OHPmzNE7\nTrrZ2dlx9erVZ7536dIlAPLly0dsbOyT71+/fv2Zdj179qRDhw5cuXKFe/fuMXTo0CdnxooXL054\nePiTtrGxsS8UT89fshw2bBiVKlXi3LlzREVFMW3aNFJSUlL1PhwcHLh8+XKqB9/nBklJSQQEBODl\n5cWsAQNYVaECOxo3JqJbN1I+/RR++QXk30tkU1JwCYPRqHdXCuY15+efZYLTjNTAsQEWJhaUa1EO\nX19fbt68qXekdKlXrx4mJib4+PiQmJiIr68voaGhaJrG22+/zfHjx/n777+Ji4vD+7kzaDExMVhb\nW5MnTx5CQkL4+eefn7zWuXNntmzZwr59+0hISGDixIn/WTzFxMSQP39+8ubNy6lTp1i0aNEzrxct\nWpSwV0zJUatWLezs7BgzZgyxsbHExcXligXH/fyCcHcfj5ubN+7u49m4ehOLf/uNQYMGPfn3MDIy\n4qqZGd8WK0bLa9cosWMHJnPnUmLMGLaMHKn3WxDijUjBJQzCnF+XcWJnMD9OGYexuYXecXIUTdMY\nUWsEK8NW0qVLF5YsWaJ3pHQxNTXF19eXFStWUKhQIdasWUOnTp1QSlG2bFkmTpxIs2bNKF++PA0b\nNnzmrNTChQuZOHEiBQoUYMqUKXTr1u3Ja5UrV2bBggX07NmT4sWLY2Nj88wdhS8blD9r1ix+/vln\nChQowODBg+nevfszbby9venbty/W1tasW7fumT6MjY3ZvHkz586dw9HREQcHhzdehim78PMLwtNz\nG9u3T+Xg3hGcvR5Mx8Hv8/H4CVSsWJHQ0FBCQ0OZMmUK3yxaRFBgIEknT3Lj+nUCr16l8ocf0n7d\nOuq0acOZM2f0fjtCpEm611LUNK0lMAcwBr5XSn313OvvA6MBDYgGhimljr6kH1lLMZe6//ABhSu/\nhfs75djku03vODlSbGIsjt84srL+Snq26ML148cxc3Z+7TbZaf2+/v37U6JECaZMmaJ3lEyRnY7F\n67i7j2f79qnY2ezkXv5B5MlrTfzt73F18cXfP3XH7vjdu6xcsIDv586lR48eTJw4kcKFC2dyciFe\nlKVrKWqaZgzMB1oClYAemqZVfK7ZeaCRUqoqMAX4Nj37FDlP5yEdKKil4LtaLiVmlrymeen/Tn/+\neBBAkrMz06ZO1TtShsoJxUhuEB9vQln7aVxP7kTRQrWIanqPuF7tOWcbSkxCTKr6qGxjw5cTJnDy\n5Ek0TaNixYpsnDULHo/jE8JQpfeSYi3gnFLqolIqEfgFaP90A6VUsFIq6vHTg0AJhHhs095NBKwP\nYP1vP2FiZq53nBxtWM1h/PDXCjp26cy3e/ZAUpLekTJMepfgEZkvKSmJiIv+REVMpnTR8VysfQKO\n9odf1xNjfRnnOc54BXhxI+ZGqvorXLgwPj4+bNu2jSFTprC6Rg24cCGT34UQby69BZc9EP7U8yuP\nv/cqA4Hf07lPkUPEJcTRs09Penj0oEHVBnrHyfFKWZeinkM93m5RnJvXrxOyYoXekTLM8uXL+fzz\nz/WOIV7hzp07NGvWDEubZPKWHMi5huvhfDMI8qK0xW8sbbmYA4MOEBkXSYUFFRiyeQhn7qRujFb1\n6tXZsXcv/VNSGN2uHTy1dqUQhiRdY7g0TesMtFRKffD4eS+gtlLqhfUpNE1rDCwA6iulIl/yupo0\nadKT525ubri5ub1xNmH4Wn7YktBdodw8dhNjY5lzKytsD9vOp9s/JcW/JEUiI9n5eDb0l8kp44Zy\ngux8LOLj42nWrBkuLi5MnzGdBgtduXkukdLH22NhnoKHR3PatGn0pP3NBzeZHzKfRYcW0cipEaPr\njaZ2idr/uZ+VoaH0a9OGwcWKscjXF8qUycy3JXKhwMBAAgMDnzyfPHlymsZwpbfgqgN4K6VaPn4+\nFkh5ycD5qoAvj4qzc6/oSwbN5yIb92+kY/OOBO0LosE7cnYrqyilqLigIt3tP2XqwNE8OH8es4IF\nX9o2O3/I5zTZ9VikpKTQp29f4h4+ZPUvq+n5W09SVAq/dvn1P9dIfZDwgGVHlvF18Nc4FXRiVL1R\ntC7bGiPt1Rdmfvv7b7q0bEmPokVZFRwMFnLHs8g8aR00n96CywQ4DTQFIoAQoIdS6uRTbRyBnUAv\npdSB1/QlBVcu8eBhDA6VyuLe0Y3Vs1frHSfXmXdwHkGX93B6djjeo0fTqWPHl7bLrh/yOVF2PRZt\nhg3jz337CAsOxjPAk4v3LuLX0w8zE7NU95GUksTa42uZsX8GCckJjKo3ip5VepLHOM9L2+84dYpW\nzZvTtmdPNnz11UvbCJERsrTgerzDVvz/tBBLlVLTNU0bAqCUWqJp2vdAR+Dy400SlVK1XtKPFFy5\nRP2OLfj7WBh3j58gj1nqf/GKjHE//j7Oc5yZZD2JbRu28fvvLx9WmV0/5HOi7Hgsxnz2GbOWLWP/\nzp0citvDsiPLCOwXiGUeyzfqTynFH+f/YOb+mZy4dQLP2p4Mrj4YK3OrF9oGnz9PK1dXls6ZQ+fO\nndP7VoR4qSwvuDKKFFy5w7er5jF0uBe//Pgt73XornecXGvE7yOw1Cz5vvf3/Pnnny9dizA7fsjn\nVNntWKz67jv6jBrFivnzadGhGVUWVWF3v91UKlwpQ/o/cu0IM/fPZFvYNgZVG4RnHU+K5y/+TJuQ\nkBDatm3LgQMHKFWqVIbsV4inScElDFZ09D2KudSgVrWK7FqzWe84udqp26dwW+FGl0tdsLG2eekd\nftntQ/6/uLm50bt3bwYOHPjCa5cvX6Zy5crcv3/fIKeXyE7H4tCuXdTp3JmPPvqIWRMn0nN9T5wL\nOvNF0y8yfF8X711kdvBsVh1dRccKHfm03qdULPz/U0HOmTOHn376ib1792ImZ9NFBsvSiU+FSItm\nPTpjbKSx/SdfvaPkehVsK1ClaBVKuJVg6dKlJOWgOble5XVzdTk6OhIdHW2QxVZ2Enf3Lq26daNJ\nhw7MmjiRHWE7CL4SzPhG4zNlf84FnfFp5cNZj7M4F3TG7Qc32q1ux97Le1FK4enpib29PaP79IFU\nLiouRGaRgktkiZ//+Ikju/bz87zZmJqa6h1HAB61PNgQuQF7a2uWDx+udxyRA0ybOJH6Njb4L11K\nXFIcw39m6DSrAAAgAElEQVQfzvxW88lrmjdT91sobyEmuE7goudFWpVpRb8N/WiwvAF7Lu9h+bJl\nbNq4kTkt2jyzaLaf36unRBEiM0jBJTLd/dj7DBgwkP6j+9C2xbt6xxGPtSnbhusx1ynargWTgoIg\nLk7vSKkWHh5Op06dKFKkCLa2tnh4eODt7U3v3r2ftLl48SJGRkakPHVm49y5c9SuXRsrKys6dOhA\nZGTkS9vevXuX/v37Y29vj42NDR1fcSen+H9Hjx5lyZo1LNy5EyNN48u9X1K1aFXalGuTZRksTC0Y\nVnMYp0ecZliNYfTd0Jc+2/rSa+QHjDq4lzOH67F7tzfbt0/F03ObFF0iS0nBJTJdyyEtsbG1YfHE\nxXpHEU8xNjJmeM3hmL99l+sRERxeskTvSKmSnJxM27ZtKVmyJJcuXSIiIoLu3bv/5+VApRQ//vgj\ny5cv59q1a5iYmDBy5MiXtu3duzdxcXGcOHGCmzdv8vHHH2fGW8kxkpOTGTRoEF988QXFixfnzJ0z\nzA+Zz9yWc3XJY2xkTK+qvTj14SmaODfhqzzLKO5Sixv5RmBi/GilubCwacybt0OXfCJ3koJLZKql\nfks5uOEg29dul/ExBmhgtYH8cXEjZV0b8dn69Wna1vvCBbTAwBce3q9Yz+5l7V/V9nVCQkK4du0a\nM2fOxMLCgjx58lC/fv3/HFSuaRp9+vShUqVK5M2blylTprBmzZoXtrt27Rr+/v4sXrwYKysrTExM\naNiwYZpz5iY+Pj5YWloycOBAlFIM9xuOV0MvShTQd+lcMxMz/lf3f9Q64MFl4xqYW5pT9J3/X+43\nLk5WuBBZ5/VT/QrxBvz8gvDx2U5MXALBx7+mfe8uvFX6Lb1jiZcolLcQnSp04r6VDRsnrCD5/HmM\nU3kLvXfJkniXLJnqfaW1/auEh4fj5OSEkVHa/150cHB48rWjoyOJiYncvn37hf5tbGywsnpxfifx\nosNnzjBl2jRCDhxA0zRWH1vN7djbeNR+YYU33eQzNoGAqZg4vUNE+FCMHXaRHN4Yc/NkvaOJXETO\ncIkM5ecXhKfnNrZvn8Kl+zsxs3Dg6LbSMlbCgI2oNYLg2F/Q8uXj22xwWdHBwYHLly+TnPzsh6Wl\npSWxsbFPnl+/fv2FbS9fvvzM16amptja2r7Q/927d4mKisrg5DmPiomhTefOvD1gAGXKlOFe3D0+\n2f4Ji9su/s+le7LSyJEtKF3ai1uXeqAa1CfZwoOSpcbi4dFc72giF5GCS2QoH5/thIVNo0z93lwL\nv4jJg/WcP/+FjJUwYNXsquFc0ImGHdz4IyJC7zj/qXbt2tjZ2TFmzBhiY2OJi4tj//79vPPOOwQF\nBREeHk5UVBTTp09/ZjulFKtWreLkyZPExsYyceJEunbt+sKlbjs7O1q1asXw4cO5d+8eiYmJBL1m\nke/c7NM+fbgbG8uGKVMA8Arwol35dtQpUUfnZM9q06YRc+e64+4+gWqaAzy4yLsd1DOLZguRWn5+\nQbi7p32qEym4RIaKjzfB6q0lnD/5OyWLfkhMZHVAxkoYOo9aHsSWPMPOLVuIiYnRO85rGRkZsXnz\nZs6dO4ejoyMODg6sWbOGZs2a0a1bN6pWrUrNmjV59913nymm/h3D1a9fP+zs7EhISMDHx+eZ1/+1\ncuVKTE1NqVChAkWLFn2mnXjkvK8vcwICmO3jg5WZGSFXQ/A95cv0ptP/e2MdtGnTCH//KRzesoSG\nQ1xY8eMiOYsp0uz/r+JMTfO2MtO8yFB12g/k1PlgrEyLcfnIziffd3efgL//FB2TiddJTE7Eea4z\nZbaWYWDvgfTp0ydbzW6e0xncsUhOpoqLCynlynF87VqSUpKo9V0tPq77Mb2q9tI73X86fO0wDTs2\npH+t/sz3ma93HJGNuLuPZ/uOiRQo3YT75/bJTPNCH9djrnMn+Sjx9x9w7e//n02+dOlxMlbCwJka\nmzKk+hDMapixYsUKveMIA3fgq684fvkyvgsWALAgZAEFzQvyfpX3dU6WOi52LjTo34AfV/3A38HB\nescR2Uh8vAn2NduSknw1zdvKGS6RIR4mPqTWF7UImxnGlDHe7NgbTVycMebmyXh4NJexEtnA9Zjr\nVJhbAeNvjDn852GcnZ0N66xKLmZoZ7jaly9PrRYt8Jo3j6v3r/L24rfZN2Af5W3L6x0t1UKvhjKz\neSOuRBZk79Wrb3TXq8h9XJp35thfQVgXnM6tcx/I4tUiayml6Lq6K3+M/YOvJ3790sWBRfbQc31P\nrq4Iw9XUnim//WZQH/K5mSEVXPv376dHjx6cPn0ac3Nz3lv7HhVsK/B54xcXQDd07y1sSuj4wzT3\nGMm3kyfrHUcYuCObN9O0Qwcc3Svw99bjgCxeLbKY104v9i3bR9OaTRkwYIDecUQ6eNTy4Hz5BL7a\nt0/vKMIAKaUYM2YM3t7emJubs/XsVg5fO8zYBmP1jvZG/vfuVN5pbMdSHx8u3LihdxxhwG5cuECH\nzp1xdTHBc9Ak3N0npLkPKbhEuiwKXcSKX1dgcsaE77/7XmaTz+bqlKiDbVmNlDx59I4iDJC/vz+3\nb9+md+/ePEx8yIitI1jQegEWphZ6R3sjdR3q8vDd4pQsVZJeH32kdxxhoBISEujcqxeVqtryzoyx\n9O/03hvdBCYFl3hjm05vYqKvN/G/xbN2zVqsra31jiTSSdM0PGuPpEStMnpHEQYmJSWFsWPHMm3a\nNExMTJi2ZxrV7arjXsZd72jpMqHJZCq9HU2wnx//XLyodxxhYJRSjBgxAvJpnOhpyugGn71xX1Jw\niTdy4MoBBqwbSszP5owY+AF16hjWRIfizXWr3I17b4XrHUMYmHGzZ3MT6NChAydvnWTJn0uY03KO\n3rHSrb5jfR64OlCyUR16T5qkdxxhYEbOns2+/fu42+ous1vOTtfZXCm4RJqdvXOW9r90wNzfmeK2\ntkz6PPsNlhWvZmFqwdCGXfWOIQxIXGgo38yYwRAvLwCG+Q1jYqOJFM9fXOdkGWOSmzcJDa5xYuNG\nLl26pHccYSBWHD7MwilTaDOhM8ULFadTxU7p6k/uUhRpcvPBTeotrYf90RqE+O8lfPdubEuX1juW\nyGCX7l3C2VqmhTAUet+lOKBJE9ZHRXHv0CFWHl2Jz0EfDg46iLFRzllBwm2FGwX2FqAIRfj+++/1\njiN0di8hgWJ169KlVXO2Wy1jV99dVC5S+Zk2j38u5S5FkfEeJDzg3dXv4vKwCXvWbmXLjBlSbOVQ\nTgWd9I7wWs7OzuzcufO/G4p0u7dtGz/+9RdfzZhBZFwko3eMZnHbxTmq2AKY6DqRE+VOsHHjRs6c\nOaN3HKGzdqNHk//WLfLWuEOPt3q8UGy9CSm4RKokpSTRfX13nPOUYsPXGxjTsydNe/bUO5bIpV53\nxicpKSmL0+Rsw2fOxKZiRYY2bcrYP8bSpVIXahSvoXesDNfYuTF2he1o0rMJk2QsV662+uBB9i1b\nxkKPvmw6t5nJjTNmjjYpuMR/Ukrxod+HxCXEEfVzJB+0bcsXCxfqHUvkUr179+by5cu8++675M+f\nn5kzZ2JkZMSyZctwcnKiWbNm7N69GwcHh2e2c3Z2JiAgAHj0//SXX35JmTJlsLW1pVu3bkRGRurx\ndgxa4t9/43/gAFMmTGB/+H42n9nMtCbT9I6VKTRNY2Kjifzt9DeBfn78LWdQcyWlFJ8MHMj7NWsy\np9BOpjaZSkHzghnSd7oLLk3TWmqadkrTtLOapr1wv6SmaRU0TQvWNC1O07RP0rs/kfWm751OSEQI\nZULLkJyczJwlS0Dm2xI6WblyJY6OjmzZsoXo6Gjee+89AIKCgjh16hT+/v4vPfuladqTeeJ8fHzY\ntGkTQUFBXLt2DWtraz788MMsfR/Zwep9+6hWvjwDmjdl6JahzHafjZW5ld6xMk2zUs2wtrJmYLkS\njPrgA73jCB0s/+gjip09i/vk94hPiqf/O/0zrO90FVyaphkD84GWQCWgh6ZpFZ9rdgfwAGalZ19C\nHz/+/SPf/vktXaK7sHvnbtauXYupqanesYQB+LeASe8jvf4trry9vbGwsMDc3Pw/t1myZAlTp06l\nePHimJqaMmnSJNatW0dKSkq68+QUKSkpzFi4kM+mT8fnoA/FLIvRrXI3vWNlKk3TmOQ6iV0dk9l5\n/z5LV6/WO5LIQlfPneOz+fNZOHc2n4VMwaeVT4aOVUzvGa5awDml1EWlVCLwC9D+6QZKqVtKqUNA\nYjr3JbLYjrAdjNoxijFFxzD/6/ls2bKFggUz5tSqyP6UUhnyyCjPX0J8nYsXL9KxY0esra2xtram\nUqVKmJiYcEOWd3ni999/J0+ePFSoVYHpe6ezoPWCXLGShHtpdxKLWNG6ST1GT0n7bOIie1JKMfR/\n/2P40KFsLHOFJiWbUM+hXobuI70Flz3w9AyJVx5/T2Rzf13/i/d932e47VhGDf2M9evXU6pUKb1j\nCQHw0g/+p7+XL18+YmNjnzxPTk7m1q1bT547Ojri7+9PZGTkk0dsbCx2dnaZGzwb+eqrr/jss8/w\n3ObJiFojKFuorN6RsoSmaUx0nci12heIun2bZatW6R1JZAFfX18uXLjAe+OH893h7/iy2ZcZvg+T\ndG6foRPDeHt7P/nazc0NNze3jOxepNKle5do+3NbPq0wkXEDP+fjbt2oVy9jK30h0qNo0aKEhYXR\npEmTl75erlw54uLi+P3332nevDlffPEF8fHxT14fOnQo48aN44cffsDR0ZFbt24RHBxMu3btsuot\nGLSdQUFERERg9pYZx3ceZ3Xn3HVprU3ZNkyynERTt7p4zZ3LgF699I4kMtGl+/f5+NNPWb50KZ/t\n+ozR9Ue/dFLfwMBAAgMD33xH6bwUUAfwf+r5WOCzV7SdBHzymr6U0N/d2Luq4vyK6vPtXyrzMmVU\nK3d3pVJS9I4ldGDIP5MbN25Ujo6OytraWs2aNUsZGRmp5OTkZ9qsWLFC2dnZqSJFiqhZs2apkiVL\nqoCAAKWUUikpKWr27NmqfPnyKn/+/Kp06dLKy8tLj7eSKll6LBISVPG6dVWv6VOV0zdOakfYjqzb\ntwHZcHKDquJTQxkVKqQ2HDigdxyRiap//LEq26yZ2nJ6iyo3r5yKT4pP1XaPfy5TXTOla6Z5TdNM\ngNNAUyACCAF6KKVOvqStNxCtlPr6FX2p9GQR6ReXFIf7KneqFKrOau/d2AIng4MxypNH72hCB3rP\nbi7+X1Yei99nzqTt9Ol8uHIQd+Ov8lOnn7Jkv4ZGKUW1JdUofNqFfBEP2PDrr3pHEplg+9mztKxZ\nk+C9u+m9qytzW86lVdlWqdo2rTPNp3tpH03TWgFzAGNgqVJquqZpQwCUUks0TSsGhAIFgBQgGqik\nlIp5rh8puHSUolLosb4HKckphH4TTnJMDOf27MHM2lrvaEInUnAZjiw7FkpRrnp1bN+qwNkqOzg2\n7BjFLItl/n4NlO9JX6bumEr4lHAOHDhAaVlZI0dRSlG8SROq2NrS7KOaBF0KYkvPLanePssLrowi\nBZe+Pt3+KQfDD+K825mLZ8+xdeNGLIsU0TuW0JEUXIYjq45F6MqV1P7wQ1ymV2dAra4Mrzk80/dp\nyFJUCm8vfpsqx6uQPzE/S5Ys0TuSyECzFy9m9GefcSR0F403tODAoAOUsSmT6u2l4BJpNvfAXBYf\nWkzdo3U5e/Is/v7+5MuXT+9YQmdScBmOrDoWNRo04L6VBQU7RBE8MDjHrZf4JtYcX8NX27/iwpQL\nHDt2DHt7uRE/J0iOi6NA9eoMdHPjXrNo7PPbM73Z9DT1IYtXizRZf2I9M/bNoNHpRhw7cgw/Pz8p\ntoTIhe4fPszJ0FAia1zMkYtTv6kulboQZxZHq+oVmN2mjd5xRAbZPGQIJa5cocsn77Hzwk68Gnll\n+j6l4MrF9lzawzC/YbS73o7g3cFs27aNAgUK6B1LCKGDFQEBlKpkS886rXGxc9E7jsEw0owY33A8\nVxvG8v0//3Dh0CG9I4l0ij9/nk9XrWLuN7P53x8f81Wzr7DMY5np+5VLirnUyVsnabTclaKB73D3\n6AmO/PknRYsW1TuWMCBySdFwZPaxSElJoVSZUtxrdY9LX1/K0eslvonklGQqL6xM3NoilLbIS4C/\nv96RRDrMmjSJ3Vu20O7boaz4ewV7++99o1UU0npJMb0Tn4ps6Fr0NVr/2AbzTaW5HHGOv/z9pdgS\nL5UblnLJzfz8gvDx2U5ExAWuxtzgvbd7SLH1EsZGxoxvNB6fez+wa2YwN48epUjVqnrHEm8gMjKS\nrxYuxO8PP97d8S5b39+aZb/n5AxXLhMdH43bwsZc/CEB48Rkjm/bRuESJfSOJbKRlAcPyFOyJP/7\nYiQzB43XO454Q35+QXh6biMsbBqYukKzP3G+MIT5s9rTpk0jveMZnKSUJCouqEj0zzbUK1wU302b\n9I4k3oCXlxc3btwgX+d8xCXFseTdN7/zVAbNi1dKTE6ky+JOnF58HytTMy4GB0uxJdLMKF8+qlYu\nz9pf/PSOItLBx2f7o2IrzwkwCoGYj7h46mvmzduhdzSDZGJkgldDL0o0K8im4GDiEhL0jiTSKPji\nRb5esIBuw7ux+p/VTGs6LUv3LwVXDufnF4S7+3gaNh1H8UGV2D/9AG85OHJ2/37yygB58YbGdelK\n+F9nuRh5We8o4g3Fxz8aUWJRaSyaW3kI/QSAuDi5O/FV3q/yPpGFz2Hu6MDYZcv0jiPSqP/EibzT\nvj1f/vMlE10nYpvXNkv3LwVXDvbvJYPtxzqx12YVt3+9hmlyNcZ/NAFjU1O944lsrMvgYeSxsODT\nBVn7F6LIOGZmSZhZXCXu3E6KxNeFuEerSpibJ+uczHCZGpsyrsE4HJrmJ2DJErmpJBsJ8PXlzMaN\n9OrfkJsPbjK0xtAszyAFVw4212cbYUUKQmVX+C0KHv5E5O0g5s//Q+9oIrszNaV+lUrs2n5IPnSy\nqZEjW+BYtR82pRy5cfBR4Vy69Dg8PJrrnMyw9X67N7HFLhIdc49du3bpHUekRnIyH8yaRcPmzZn1\n11R8WvpgYpT19wxKwZVDXYu+xt8lf4WbM2BnCUj+E2gPyCUDkTFmzJnHvcNHCTofpHcU8QbcW9bn\nwuXjlCthhWstH9zdJzB3bksZMP8f8hjnYVyjcVi6WTJr1iy944hU2DJ/Phf/+QeXXmWoZV+LxiUb\n65JD7lLMYRKTE1l1dBUzVkznvG8slnHluXt7M5D3SRt39wn4+0/RL6TIMUq9U4oybcqwfdp2vaOI\nNBo9cxrz5i3m0rEQiljZ6R0nW4lPiqfM7DI8/Pohu7dup7KLTBRrsBISaFC+PBb163C46nYODz6M\nU0GnDOla7lLMpR4mPmRByAJKzS3Foi+XcfqHmzSrVRer/LV4utiSSwYiIw3pP4TALYE8SHigdxSR\nRlvmz6NlVWcptt6AmYkZYxuPpVyVfMxu3hzkZIHBOjF5MmciIsjbLBbP2p4ZVmy9CTnDlc3FJMSw\n+NBiZgfP5u3Etwn96QoPoqL47vPP6dWvH35+Qcybt4O4OGPMzZPx8GgulwxEhrlz5w7FHIox/4/5\nDKk3RO84IpWCDgfRvJ4rYX+GUqJyDb3jZEtxSXGU/KIk0VPvcOCHH3irRw+9I4nnpaTQxcqKYp2b\n41ftCCeGn8DC1CLDuk/rGS4puLKpyIeRzAuZx/yQ+bg5ulHqXCm+mT4f1+bN2fDDDzLlg8gyNRrX\n4EHZB5z89qTeUUQqVe9UHSMTI0LXhOodJVvzOejDtLGrqZBowu49e/SOI55z+PBh2rRujfVnNkxr\nOY2OFTtmaP9ySTGHu/ngJmP+GEOZeWW4cO8Cy+st58rcK+zdupe/goPZ/ttvUmyJLOXZrgeR/pc4\nH3le7ygiFe7cv8ORbUeYNVYGfKfXBy4fkFI/mj1Hj3JHCi6DM3nyZOr1rI+9rT0dKnTQO44UXNnF\nlftX8NzqSYX5FYiOj8a/nT95t+elb/u+9OrVi6CgICpWqaJ3TJELtWnoyo27RswI+FbvKCIVRs0Z\nhW1ZW1yrueodJduzMLVgXOuB5C9bkhHTp+sdRzxl/6FDHAw5yG7b3cxtOdcg1oWVgsvAhd0NY/Dm\nwVRdVBUTIxM2tNvBwe+u0qJOc8zNzTl58iTDhw/HyEgOpdCHTY0alCxbls2rA0lRKXrHEa+hlOKX\nFWsYPnS43lFyjCE1hmDa0Ih1hw+TmJiodxzxWKcxYyjw7jv0culFpcKV9I4DSMFlsI7fPE4v317U\n/r42xSyLEdDzICErb+BWuzHx16+ze9Mmvv76awoXLqx3VCEYVK8ed47eYtcFmQjSkH2+dglx0WaM\nGzhO7yg5Rl7TvIzu3ANlU4BJP/2kdxwBrA4O5taRw0SVOoa3m7fecZ6QQfMG5vC1w0zbM429l/fy\nUe2PKJdcgwkz53AqKIhSZcqw/KOPaNi3r94xhXhG3Jkz5KtVi6bTWrH9w9V6xxGvUNytDk6aCcG7\n9uodJUeJSYjBbmhlSv1jw98hR/SOk7slJ2Nfty6JZcz4clR/BlQbkGm7kkHz2dTey3tp9VMr2q1u\nRz27eswoNgP/8f4Me+997B484M9Fizj3559SbAmDZF6uHG9VqMARv7NExUXpHUe8xOHwU1w79A+z\nPvhA7yg5jmUeS8b0GUzYpbMcOHBA7zi5mu/8+Vw/c4bidePo904/veM8Q85wZRE/vyB8fLYTH2+C\nmVkSI0e2oHXrhuw4v4Npe6Zx+d41OuRpS/zxOHzX+1KpUiWGDRtGhw4dMJWFpkU2sGrFCj7w+oQ5\nG79gSA2Zk8vQ1BnQiUshp7l27B8wgAHEOU10fDR2Xe1oaNSQrRu26h0nd1IKxxo1SCkG65YsoE6J\nOpm6O5mHywD5+QXh6bmNsLBHC8Ricw6b+kPIX+cGRucLEncqmetHj1PFwYFu779Ply5dKFeunL6h\nhUgjpRT2zvYU6lOIY1OO6R1HPCU+MZ78ZcsyunULpi78Xu84OdaErROY8d4MTh84jHPlynrHyXVO\nffstNTw9ab20K2t6/pjp+8vygkvTtJbAHMAY+F4p9dVL2vgArYBYoJ9S6oWL3Dm54Grh7sWOo12h\n1CoosJYikRbcv51E3NWr2Do60qJyZUZ16sQ7XbuCmZnecYV4YxMmTuCbwG8IXR9KxcIV9Y4jHpu6\nZDwzxiwg8swZjOVGm0wTFRdFveqFaRpTBJ9LV/SOk7soxfuFC3HK/gFb9l7ELn/mL1mVpQWXpmnG\nwGmgGXAVCAV6KKVOPtWmNTBCKdVa07TawFyl1Avn+TRNUy1aeDFyZItMX3rmZZf33mSficmJ3Iq9\nxfWY608ed25eIv70Ka4ej+BYzF1ORd3g7sVI1FUjiDEGk3ewK1iBInlLUcb2OusOLJDT+yLHOHPm\nDC51XRj20zBmtpypdxzxWMmahXm3YCl8dhzUO0qON2rVSBYNWMyZffsoXrOm3nFyjXP791OzUQM+\n2TyJ8a0mZck+s7rgqgtMUkq1fPx8DIBS6sun2iwGdimlfn38/BTgqpS68VxfChSlS3sxd657phVd\nL1zeg2f2maJSiHwY+UwRdePBDW5FXuX2vQjOP7xJxM0Ibty+QdTtKCzjLCl51Yb46Lzci00gKiaG\nh3fuYFyoEKXti+JarxZ7AsM5dewroDKPTgQ+4u4+AX//KZnyPoXQS1WXqlytcZUbi29gYmSid5xc\nb9/pfTR8pyHhF8KxL2avd5wc717cPZwbNKKBnQNbNvvpHSfXaN65OUceHiFicwR5jPNkyT7TWnCl\n97ehPRD+1PMrQO1UtCkB3HiuHc4N3yfRxIiPv/2KL39fho2poohtHhSPikKlFEopIu8lciXeBJRC\nQ4ECVAqFjBOxsTF+0k4pRUpKCrF347kTDVpKCndvR5FSNC8lizQlPjaGiMQYwowe0NF7Dnm+1IiN\njaXIQxtQNiQnJZOUmERiQgIJD+NQDx6gAVbWVtha21KtRDUcSjiQUsWCKEtLSpcpQ9XKlWlWtSol\nnlpe51GR9+tzRd44PDxapusfXwhDNKh3f2atmsK2c9toU66N3nFyvdFfj+at+m9JsZVFCpoXpPl7\ndfht5m/E376Nma2t3pFyrH+vVkVG3SP07wDG+EzKsmLrTaS34Ert6bHnK8CXbnfr1A40pTAihRjL\nfNzOb8I/NtrjDv6/C7NoRUrcv2eKNB7VlxqxJslcsDVG0zQ0tEf/NdIwilHEJZmBkTHxiUmkqATQ\njFHxZnDHAZIKUKToCQZ1bkmRgkWIeRjP5QcPyF+gAPmtrLApUgQnZ2cqFCuGc8GCaV4i4N+zdfPm\nTSAuzhhz82Q8PFpm+qVTIfTQqWRJPj4Vz6IDy6Xg0llUXBQHNhxg/ar1ekfJVZZ4fInf9zv539hR\nLPxuud5xcqRnrlaVagCly7D2qwQaFA/KtM/WwMBAAgMD37yDp88GpfUB1AH8n3o+FvjsuTaLge5P\nPT8FFH1JXwqUAqXc3cerzNKihdeT/Tz9yMx9CpGrJCcr2woVlF2PaurWg1t6p8nVRizwUAVKFFAp\nKSl6R8l1mnq+p8zKlJJ/+0zy5LPcPkRhaamw35/ln+WPSqjU10zpnfj0EFBW0zRnTdPyAN2ATc+1\n2QT0AdA0rQ5wTz03futpjy61NU9nrFcbObIFpUt7Zek+hchVjIx4r1YtYo8n8POxn/VOk2sppVjy\n826a9OlrEAv35jY/fj6X+JgYZmz6Ve8oOVJ8/OMLdCVHQ7V6cLUuAHFxxq/ZSl/puqSolErSNG0E\nsI1Ho8GXKqVOapo25PHrS5RSv2ua1lrTtHPAA6D/q/pzd5+Q6Zfa5PKeEJlv4siRLHR1ZcmenxhZ\ne6TecXKlZfvXkPhXGN/4yCUtPRQvUIyKjcvg/b9hbP3mVLruiBcvMjNLwqrUWqL+OoRmvu/JOCVz\n8+/WHNEAACAASURBVGRdc72OTHwqhMgUJV1cSHE2ZcOCxVSzq6Z3nFynfKfGJFy6x4U/ZW0/Pfj5\nBTH801+5fGkRmP8OkS0z/S783GTLlt0Mmf4xRnFmXDm8H3h0tWru3Kw7gSJrKQohDMLA9u25e/w2\ny/+SMyxZ7ULkRcIOnGJ0+/Z6R8m1fHy2c/nUArBpAraPlroKC5vGvHk7dE6WM5zPc4JrR89QzaoS\nrq7euLtPyNJi603IGS4hRKZITEzEzt6OlAEpXJt6DTMTWUUhq7Qd34c/Vuwk9uxZjCws9I6TK7m5\nebN7tzcYXQLTUuQpspKE8J64unoTGOitd7xsLTYxlmZv21IwuTi/nz6nWw45wyWEMAimpqb07N4T\nqzNWbD6zWe84uUZcUhxH1+7i3Vo1pNjSkZlZ0qMvUpyoRQlqpgwGLdmgxxhlF1O3fs6pUw+ZP+UL\nvaOkiRRcQohM06tXL+IOx7HsyDK9o+QaS/cs5c6FCL7z/lzvKLna03fER8XP5sjtJIr8X3v3HV/T\n/cdx/PXNtDcl9qjZGFW0KIlfBTFaVVSrQbX2rFq1R41WtUajVWpUSVrUSJoQ46IqihCCSIxYIbEj\n1UTG9/eHUFpU5rn35vN8PDzcM3LO+7bXyed+v9/zPS3ayB3x6XT+1nkWz/qKVuXLUqFjR6PjpIo8\nd0MIkWnq1atHHvs87Pp9F5HtInHK62R0JKs38+uZvNLWlQI1axodJVv75x3x16Irkssumldfq21w\nMss2zGcYsYdtGb93k8U9h1hauIQQmUYpRbdu3Sh1uhQ/BP9gdByrt/f8Xi7tuMTkj6V1yxy0bt0E\nf/8pmEwTmfpeV87vu8gUk2V1g5mTPef3sGnVJtxbuFOlShWj46SaDJoXQmSqs4GBVGvWjJJTyxE2\n9KhMwpmJ3Ca5ceSHI0SGR8p/ZzOjExLIX6UK+RsWZNf8NZQrUM7oSBYlWSfzomdDTow/wR+mnTg7\nOxsdSQbNCyHMS9kGDchRpgx2QfkIvBBodByrdfXOVUxrTAwdMFSKLTOk7O3p36ED0ftiGBEw0ug4\nFmfF4RVE7HWk8Ev1zaLYSgspuIQQmUspOjVpQtSRP2VOrkz06ab5JJ2xp1f3XkZHEU8wefp09J04\nNu89we/nfzc6jsWIvRvLcL9xxPgeYc6UKUbHSTMpuIQQmW7S6NHcPHWGjXs2cSfhjtFxrE5SchK+\nC36l/IsvUqBAAaPjiCewt7Nj+siRVNhvw2D/wSTrZKMjWYSZv80k/8HncSpTlg716hkdJ82k4BJC\nZLrnypenorMzhY+UYe3xtUbHsTrrj67ndPBZxnl0NzqK+A/9evbk4vGLJFxOYNmhZUbHMXtnb55l\nfuBCTm0+wJw337S4OxMfJgWXECJL9H/zTcKPnJM5uTLBnNnzcMiZE4/33zc6ivgPOXPmZMCAAZQ6\nUoox28YQEx9jdCSzNmLLCCqFvEiJ0qXpMNKyx75JwSWEyBIDPvqIoloTFBRExM0Io+NYjbBrYRzZ\nH8VbrVrKYHkLMXDgQAK37qFxrleYssNyxyRltl1nd7H75G7Or9nG6vbtwcHB6EjpIgWXECJL2Nna\n8n637pQ9W1a6UjKQ56Kx3Io4zxeT5Re3pShQoAB9ChYk109XWHJoCWHXwoyOZHaSdTJDNg3B/XAN\nXrW15eUxY4yOlG5ScAkhsoyHhwcXdl9gSdASGTCcAWLvxuLvtZ73natTtGhRo+OIVBgyeTLegfvo\nVGEAH236yOg4ZmfpoaXYJdqxYfMBxn/2mcW3boEUXEKILFSpUiWqV60O4bAjYofRcSze8qDlREQo\nhi6S6TYsTZEuXahdsyamFYcIuxaGX7if0ZHMxu3424zdNpZXLr9CYxcXnAcMMDpShpCCSwiRpbp1\n60a+4/lkTq500lozc/FMKleuTPXq1Y2OI1JLKRb37k3o1q30rTmOoZuGcjfprtGpzMK0XdNwLemK\n10Ivxo8fb3ScDCMFlxAiS3Xs2JHTgWH8tnuz3KGVDr+d+42rO64yeshoo6OINKrevTt1qlRhmacP\n5QuWZ/4f842OZLjTN07zXdB3XD9UgnL16lHTih7CLgWXECJL5c+fnyK1a1PkxPP8dPQno+NYrGnr\npmF7zZYOHToYHUWklY0Ni4YN47CfH31rj2XarmlExUYZncpQwwOG0716fzZ9u5gxEycaHSdDScEl\nhMhyvTt3JuTwOZYGyZxcaRF5O5KAXyNo+153HKxgMHF2VqdLF5xbtWLx4jV41PJg7LaxRkcyjCnC\nRNClIHZ7n6JS06a0rVvX6EgZSgouIUSWGz5wIHf/+ou7u25w4uoJo+NYnO/XTCd5z0X69ulndBSR\nAdZ8+im7ly9ncO3BbAzbyIHIA0ZHynJJyUkM8R9Cv6oj2Lt+A8s7djQ6UoZTWmujMwCglNLmkkUI\nkfmav/02pyJO0nlqc6a/Nt3oOBbjbtJdXFo1IPyK5srBQ0bHERnEw8ODKlWqUKxlMZYFL2NXj13Z\naiLbhQcWsuLwCq6tykHO+Hj27zD/u5iVUmitn/l/krRwCSEMMXPsWCJCQln1+0qSkpOMjmMxNgR5\nERx+g34ffGh0FJGBRo8ezZw5c3ir4lvcSbiDV4iX0ZGyzK24W4zfPp5hxXtxcvduVk2xzkl809zC\npZQqBHgDZYEIoJPW+uZj9vseaA1Ea62dn3I8aeESIpt5uVUrLhUN55tx82j1fCuj41iEzl3rs2Zz\nBHciI3GwszM6jshAHp06Ub5IEV4b2YV31r5DaP9QcjvkNjpWpvt488fcjLvJ7RG/UStvXj7Zt8/o\nSM8kK1u4RgEBWuvKwNaU5cdZArRMx3mEEFZqwqBBsD9B5uR6RocvB3PgUAzN2rSRYssKTW7cmPkL\nF1LRrgKNyzRm5u6ZRkfKdOHXwll6aCmdbFuyKzycwV7W27KXnhauUKCp1jpKKVUcMGmtqz5h33LA\nRmnhEkI8LCkpifIVynOj3Q3OfnaWQjkLGR3JrPX5uiMrh6zl+MVLlCxWzOg4IqMlJ9PO2ZnQ8uXZ\n+uMCan9bmwO9DlCuQDmjk2Wa171ep1HpRmyZtIX2zZvTd/hwoyM9s6xs4XpOa31/wpAo4Ll0HEsI\nkQ3Z2trS68NeFA8rzsojK42OY9Zuxt1k2Xof2r77thRb1srGhrkzZnBy504OhEYyuMFghgdYTgGS\nWltObyEkOgTnP505c+YMHwwZYnSkTPXUNmmlVABQ/DGbHnlst9ZaK6XS3Tw18aFJzlxcXHBxcUnv\nIYUQZq5Hjx7MnDWTRYGLGFDfOp6ZlhkW7VuE2q8YteNJozeENSjXti1tX3qJPsOGcWZ7ANW+rsb2\nM9txLe9qdLQMlZicyBD/Icxs9hmjuo1l6tSp2NvbGx3rqUwmEyaTKc0/n94uRRet9WWlVAlgu3Qp\nCiHSol3TplzXoXj+HEDN56znUR4ZJVkn49TNiVKnSrF/936j44hMdjM4mCKurnz1ww8UL/8Xk3dM\nJqh3EHY21jNuz3OfJ6uPrabCjXZ4L1rErcOHsbGxrIkTsrJLcQPQLeV1N2BdOo4lhMjG6r75JofO\n52LZ/sVGRzFLAacCuL3zNuNGjDM6isgCBWrVoufQoXwycSJvVGlPoZyF+O7Ad0bHyjA3/rrBpB2T\nGPHSRJZMncr8+fMtrthKi/ROC/ETUIaHpoVQSjkB32mtW6fstwpoChQGooHxWut/3ZIkLVxCZF/x\nCQnkdnLiFbeCbF0egoOtPK7mYbU/bcupBX9w82wktra2RscRWSAhIQGnqlX5fv58ytRxwm2FG8f7\nH7eKG0uG+A8hLjGOwNlhxDs6ctzPz+hIaZJlLVxa6+ta69e01pW11m735+DSWkfeL7ZSlrtorZ20\n1o5a69KPK7aEENmbo709rzZvTmR4bnzCfIyOY1Yirp8myvs0Ll3ek2IrG7G3t+fbzz9n/Cef4FzM\nmTervslE00SjY6Vb6NVQVhxewcuXnDkcFMS6OXOMjpRlrL8NTwhhEWaMGcOZ4yfx3jDP6ChmZfEX\nI4g6c475I2WwfHbTvn17HB0dWbVqFVOaTcErxIuQ6BCjY6XLsM3D+KTBcIbO/JLOXbpQpXJloyNl\nGSm4hBBmoUGNGhStXJlTW29y5sYZo+OYhbjEONZtPUnVVxtTtkgRo+OILKaUYtasWYwcORL7u3aM\nazKOIf5DsNThN/4n/Qm/Fk5hzwOUun6dZXPnGh0pS0nBJYQwG2OGDyc8LJKFBxYaHcUsrP75c44e\nP8Psydb5bDnx3xo3bkyrcuUY89pr9HmpD5diL7H+xHqjY6VaQlICQzcN5fOyA/nk59V8t3AhDmY+\nDURGS/Og+Ywmg+aFEFprqr5QleiG0UR9E5WtB89rrWnlWos/btlw/eAho+MIA90IDaVyzZpM+eYb\nKrmUobdPb472O0oOuxxGR3tmc/fOxSfMhxp7KhFz6hSLN20yOlK6ZeW0EEIIkaGUUowcNhK1V7E+\n1PK+xWekLacCCNt7lMUff2x0FGGwglWr0rNvXwZNnEi9og1wLubMl3u+NDrWM7t25xpTd06lf9n+\n/Lh6NTNWrDA6kiGkhUsIYVbi4uIoXqo4VYdVJXB0oNFxDKG15oWPXiDGJ4ZzYedQ6pm/RAsrpZOT\nKensTOUXXmDxgmk0WNSAw30P45TXyeho/2nArwNISkomePohunbtSr9+/YyOlCGkhUsIYdFy5MjB\nwP4DCV4XzMnrJ42OY4gtp7dwZuMZpk+YLsWWAEDZ2LDO05Odfn6c+P04ver2YtQW879z9Wj0UbyP\nenN0iw1RNjb06dPH6EiGkRYuIYTZiY6OpnyZUvSb8x6f985es89rrXEe7cyVH69w8cxF7Oys53Eu\nIv069O3L1qAgIkybqbGgOms6reHlUi8bHeuxtNa0WNGCwnfq4P3JIvbu20e9ChWMjpVhpIVLCGHx\nihUrRql69djgdZD4xHij42Spzac2E+ETwaQxk6TYEv+yau5ciInhu5WrmfG/GQzyG0SyTjY61mP5\nhPlwN/QsayYsY9isWVZVbKWFFFxCCLM0Y/hITgZH4BX4g9FRsozWmv4rZ5Ac7UiPHj2MjiPMkIO9\nPQHLl/PFmDG4lXDD1saW5cHLjY71L3eT7jLS7yMu+OemnLMzn8vnWboUhRDmq0D1apQqnZsi8W1w\ndExk0CA3WrduYnSsTOHru5Nxyz05FxpC/oJlmTt8pNW+V5F+o0aNIigoiEmLJtHh5w6EDggln2M+\no2M9MHvPbHZOWcWvhy5y9uhRShQsaHSkDCddikIIq+Dru5Nyuapz8ng0O450ZPPmqQwevAlf351G\nR8twvr47GTTYnyvJEVyPOM/5wO+t9r2KjDF16lSSk5NZP3cNLSq68enOT42O9MDFmItsXDCR3zcF\n4f3tt1ZZbKWFFFxCCLM0d+5mgg94k5yURNnaEwE4depT5s0LMDZYJpg7dzOnbRpx90wcZcq7khD/\nnNW+V5Ex7Ozs8PLywuubb3D1imXxwcWEXws3OhaBFwJp4lmfRO+/GNKlC+3btjU6ktmQgksIYZbi\n4+0AO+xzv8vZiwGQ7wIAcXG2xgbLBHHxtuSvP4UrJ05x69S0v9db4XsVGadIkSKs9fam3zofPP5s\nx0ebPzIsi6/vTmq815pXPZtxe2lu/ixZjpHLlhmWxxxJwSWEMEuOjokA3AmfBKeToEUzcIwhR44k\ng5NlvJjixyl8+DplKzTgZmz1B+ut8b2KjPVi69a8N2QI8xf6EXU0DP+T/lmeYd3GLbyzcjDH8p8m\ncU4/rpxw4IZui7//7izPYs6k4BJCmKVBg9yoWHEMkBMSJsHORHJ2r0Wf/i5GR8tQF2IucL6iicsh\np7A9PPbB+ooVP2HgwOYGJhOWYsGMGVSvVYsz6+0ZuLY3u87uyrJzR96OpPuO7sQklYct3eGvNcAm\nIiJmS5f4P8gkL0IIs3T/Dr1588bx55/J7Nt3lYp5KvJL0nJe182sYgb2W3G3cP/RnZrnX8CxRV60\n3kqpOBM5ciQxcGBLuUtRPLPf162jRO3a2G9/HrebbchzqxRVz7syqnenNH+OkpKTiP4zmsjbkVy8\nfZGLMRe5ePviI8vnbp2j8JV63DraGqJHQtHtEF0SkC7xf5JpIYQQFsHT05P1G9dzq8MtmldozpRm\nU4yOlC7xifG0+rEVZW3KsmHIBoKCgihbtqzRsYQFm7/0ZwaOGgyF3oIiReGVryhwugoL3xtHx3at\nHtn3dvztfxVSF2MuEhkb+WA5KjaKgjkLUjJvSUrmK4lTHidK5iv593JeJ0rnK81r7T8iaM86yLsR\nLjZ+cI4WLcbh72/Z/06fJrXTQkgLlxDCIvTs2ZNpY8cy/u3P+CxkBmXyl+HDuh8aHStNknUyPdb3\noECOAtgG2NKzZ08ptkS6bVwVDH/9AeGvw/E6cPAgN/83ga57O7Eh6Q0u3b70oLBK0kmPFE4l85ak\nUqFKNC3X9MFyibwlcLB1eOo5v/fz42DgLzQv04iAkL+LrXtd4i0z+y1bFCm4hBAWwdHRkZatWjFi\n1hz2mfxosrQJJfOVxP15d6OjpdqoLaM4d+scs1/8ijYD3QkLCzM6krAC8fF2EFMK2AF0gti+sP4n\nXkh0oNlbLz/SOpXfMX+6u+XnL1/OoMGDGdemDfW79MJm3jji4mylS/wJpOASQliM+V98wfJatVi7\n4Td+6fwL7Va1w+9dP+o61TU62jObu3cuG05s4Lcev1GtZXs6Dh5MgQIFjI4lrMD9O3shD7Ae6IUd\njSkc+yo96mTco3Xi4uIY0aMHG1avZkXnzryzbBnY2kqB9R/kLkUhhMXIUbw4fdzcmPLZZ9R3qs/C\ntgtp59WOMzfOGB3tmaw5toaZu2fi39WfRRsCuHXiBJ+9+67RsYSV+PvOXgB74HtaOEZxMvA7/Pfu\nzZBzHDt2jAZ163L5l1846OnJOytWgK0Mjn8WMmheCGFREi5fJt+LL/LxxIlM6dWL+X/MZ/4f89n9\n/m4K5ypsdLwn2nV2F2/+9Cabu26mRlFn8lerxoCaNfl8zRqjowkr4uu7k3nzAh507Q3q68reLycz\nJSiIkg0bMmvaNDq/+GKqj5ucnMyiRYsYM2YM06dPp2f79qjC5vvvLSukdtB8mgsupVQhwBsoC0QA\nnbTWN/+xT2lgOVAM0MBCrfXcJxxPCi4hxDMZ7eHBnMBAYkNDsbGxYUTACH4//ztbPLaQwy6H0fH+\n5fiV47gsc+GH9j/gVtGND774gp/mzuXmjh3YlCtndDxh7bTm8qpVfDB7Nn7h4RR5+WWmTp7MB/Xr\nP3Uc152kJFaGhPCttzchGzdSycaGn729qVq1ahaGN19ZWXB9BlzVWn+mlBoJFNRaj/rHPsWB4lrr\nQ0qpPMAB4A2t9fHHHE8KLiHEM0m6e5cX6tZl5qef0q5dO5J1Mu+ufZfE5ES83/LGRpnPaInI25E0\nXNyQya6T8ajlwdU//6R4xYp4tmxJr6VLjY4nspPERK4tWECfvXv5xdeXl11d6eTiQtGiRR/543fp\nEnO8vDjm54c6d466rVrRr0sX3nZ3x8Hh6XctZidZWXCFAk211lEphZVJa/3UslcptQ6Yp7Xe+pht\nUnAJIZ7Zr7/+ypAhQzh8+DA5cuQgPjEetxVu1C1Rl9ktZhsdD4CY+BiaLGlCpxqd+OTVTwCYMno0\na5ct42BwMBQtanBCkV3dunWLJUuXcurkSa5cuXLvz7FjRMfGEmtjQ10nJ/rHxPBGu3Y4zJsHdnKP\n3T9lZcF1Q2tdMOW1Aq7fX37C/uW4d69qDa117GO2S8ElhEiVt956i2rVqjFlyr3JFW/8dYNG3zei\nd93eDH55sKHZ7ibdxf1Hd54v9DyerT1RSnH9+nWqVqnCjmnTqPahZc4hJqzYokUQFgaJifD669C4\nsQyIf4oMLbiUUgFA8cdsGgMse7jAUkpd11oXesJx8gAmYKrWet0T9tETJkx4sOzi4oKLi8szvAUh\nRHYVGRlJrVq1MJlM1KhRA4CzN8/S6PtGzGk5hw7VOxiSS2uNxzoPYuJjWNtpLbY2935pde/enbx5\n8zJv3jxDcgkh0s5kMmEymR4sT5o0KUu7FF201peVUiWA7Y/rUlRK2QM+gJ/W+qunHE9auIQQqfbN\nN9/ww5Il7Pz9d2xTvo0fvHQQtxVurOu8jkZlGmV5ptFbRrM9Yjvbum0jl30uADZs2MDQoUMJDg4m\nT548WZ5JCJGxUtvClZ6RpRuAbimvuwH/arlK6WpcDBx7WrElhBBp1evdd7kUEcEbY8c+WFenRB1W\ntF9Bh586cOLqiSzN47nPkzXH1+Dzjs+DYuvatWv07duXJUuWSLElRDaVnoJrBtBcKRUGNEtZRinl\npJTyTdmnEdAVcFVKHUz5Iw9XEkJkGJu8eVnYqxe+np7sCA9/sL5FpRZM/9903Fe6ExUblSVZ1oWu\nY+rOqfh39adIriIAJGtNAw8P2rz1Fk2ayEzcQmRXMvGpEMLyJSXRsmFD9tvbE71rFzYPzS000TQR\n33BfTN1M5HbInWkR9pzfQzuve48aesnppQfre3z9Nd6ffkr0/v3kcXLKtPMLIbJWVnYpCiGEebC1\nZc28efwZHk7/r79+ZNOEphNwLuZM59WdSUxOfMIB0ifsWhjtvduz/I3ljxRbv50+zbJx41hev74U\nW0Jkc9LCJYSwGkt69uSDjRsJCQmhWrFiD9YnJCXQZlUbyhcoz4LWC546u3ZqXY69TMPFDRnbZCzv\n13n/oXMm8ZyLC/Vu3mSTyQTZ/DEoQlgbaeESQmRbPebPp+VrrzHs/fdJSkp6sN7e1p6fO/5M4IVA\npv82PcPOF3s3ljYr29CtVrdHii2Ad2bPJuHUKdZPmSLFlhBCWriEENYlISGBli1b4uzszFdfPXpz\n9P3H7ExtNpWuNbum7zxJCbTzakepvKVY2HbhI61mwcHBvNqoEV5t2uDu5ZWu8wghzJO0cAkhsjV7\ne3tWr17Npk2b+Pof47mc8jrh+44vwzYPY9uZbWk+h9aa3j69sVW2LGjzaBflyZMncXd35/sPP8R9\n4cI0n0MIYV2khUsIYZVOnz5No0aNWLx4Me7u7o9sM0WY6PRzJ7Z6bMX5OedUH3vC9gn4nfRje7ft\nj9z5GBkZSePGjRk1ahS9evVK93sQQpgvaeESQgigQoUKrFm9mm5vv82PW7c+ss2lnAtzWs6h9crW\nXIi5kKrjLjywkB+P/IjPOz6PFFs3btygRYsWfPjhh1JsCSH+RQouIYTVatiwIYP/9z88unRhe2jo\nI9u6OHdhQP0BuP/ozq24W890PJ8wHyaYJuDf1Z9iuf++C/L09eu4tGqFm5sbo0aNytD3IISwDtKl\nKISwbomJdHRxYf2lS/y0ejVv1KnzYJPWmoF+Awm9Gsqv7/6Kg63DEw+z7cw23l79Nj7v+FC/ZP0H\n66/cuUNFFxeeL1KEfT4+2NjI91ghsgPpUhRCiIfZ2fGTvz9dSpbkTVdXus2YQXLKlzulFHNaziGP\nQx4+2PABT/rSt/LISrqs6cLPHX9+pNjaEBJC+caNKXLjBoHTp0uxJYR4Irk6CCGsnsqTh2U7drB2\n8GC8v/ySZu+9x507dwCwtbFlZYeVhF0LY9z2cY/8nNaaz3d/zqgto9jmsY2m5ZoC8FdiIi2GD+eN\nRo1ooxSh3t7Y16qV5e9LCGE5pEtRCJGtXLtxg/79+hFy5AheXl688MILAFz58woNv29Iy7yvE7Yq\nB3HxNkRU9UNVuMZvfXZSKl8pAMLDw+ns7s7Z+HjW9ejBq+PHg62tkW9JCGEA6VIUQoinKFywIKtW\nruTjjz/G1dWVkSNHYjKZyGeThxFOk1hw/Bs2n6/DzueOcS4uD2pJe4J3nebKlSvMmjWLV155BQ93\nd6IDA3l10iQptoQQz0RauIQQ2VZ4eDhLly5ly+bNHD50CIfnq1AwT2WoehibqErcvlKYqydvkSNh\nF46OChcXF7744gsqVqxodHQhhMFS28IlBZcQQgDnjx/nnZ6jOH37Jlevn8YxV17y5yzHhcudqV/+\nAHv2fCWD4oUQD6S24LLLzDBCCGEpSlerRq68zkTumQrAXeB2yraCL56UYksIkS5yBRFCiBSDBrlR\nseKYR9ZVrPgJAwc2NyiREMJaSJeiEEI8xNd3J/PmBRAXZ0uOHEkMHNic1q2bGB1LCGFmZAyXEEII\nIUQmk2khhBBCCCHMjBRcQgghhBCZTAouIYQQQohMluaCSylVSCkVoJQKU0ptVkoVeMw+OZRSe5VS\nh5RSx5RS09MXVwghhBDC8qSnhWsUEKC1rgxsTVl+hNY6DnDVWtcGagKuSqnG6TinEJhMJqMjCAsi\nnxfxrOSzIjJTegqudsCylNfLgDcet5PW+k7KSwfAFriejnMKIRdFkSryeRHPSj4rIjOlp+B6Tmsd\nlfI6CnjucTsppWyUUodS9tmutT6WjnMKIYQQQlicpz7aRykVABR/zKZHpmLWWmul1GMn0dJaJwO1\nlVL5gU1KKRettSmNeYUQQgghLE6aJz5VSoUCLlrry0qpEtxrvar6Hz8zDvhLaz3rMdtk1lMhhBBC\nWIysenj1BqAbMDPl73X/3EEpVQRI1FrfVErlBJoDkx53sNSEFkIIIYSwJOlp4SoE/ASUASKATimF\nlRPwnda6tVKqJrCUe2PFbIAftNafZ0RwIYQQQghLYTbPUhRCCCGEsFaGzzSvlGqplApVSoUrpUYa\nnUeYN6VUhFLqsFLqoFLqD6PzCPOhlPpeKRWllDry0Lr/nKBZZE9P+LxMVEpdSLm+HFRKtTQyozAP\nSqnSSqntSqmjSqkQpdSglPWpur4YWnAppWyB+UBLoDrQRSlVzchMwuxp7t2sUUdrXd/oMMKsLOHe\nteRh/zlBs8i2Hvd50cDslOtLHa21vwG5hPlJAIZqrWsALwP9U2qVVF1fjG7hqg+c1FpHaK0TUehE\n1QAAAexJREFUAC/gdYMzCfMnN1iIf9Fa7wJu/GP1M03QLLKfJ3xeQK4v4h+01pe11odSXscCx4GS\npPL6YnTBVRI4/9DyhZR1QjyJBrYopfYrpT40Oowwe880QbMQDxmolApWSi2WLmjxT0qpckAdYC+p\nvL4YXXDJiH2RWo201nWAVtxr1n3V6EDCMuh7dwjJNUc8zQKgPFAbuAR8YWwcYU6UUnmANcBgrfXt\nh7c9y/XF6ILrIlD6oeXS3GvlEuKxtNaXUv6+AvzCvW5pIZ4kSilVHCBlguZog/MIM6a1jtYpgEXI\n9UWkUErZc6/Y+kFrfX/e0VRdX4wuuPYDzyulyimlHIDO3JtQVYh/UUrlUkrlTXmdG3ADjjz9p0Q2\nd3+CZnjCBM1C3JfyS/O+9sj1RQBKKQUsBo5prb96aFOqri+Gz8OllGoFfAXYAou11tMNDSTMllKq\nPPdateDeUxJ+lM+LuE8ptQpoChTh3niK8cB6HjNBs1EZhfl4zOdlAuDCve5EDZwBej80RkdkU0qp\nxsBO4DB/dxuOBv4gFdcXwwsuIYQQQghrZ3SXohBCCCGE1ZOCSwghhBAik0nBJYQQQgiRyaTgEkII\nIYTIZFJwCSGEEEJkMim4hBBCCCEymRRcQgghhBCZTAouIYQQQohM9n/XjwQf8XFARAAAAABJRU5E\nrkJggg==\n", "text": [ "" ] } ], "prompt_number": 65 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sparse matrices\n", "\n", "Sparse matrices are often useful in numerical analysis dealing with large systems.\n", "Whenever the problem involves vectors or matrices with mostly zero values then one should take advantage of the sparsity.\n", "SciPy has an extensive library for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create a sparse matrix we have to choose the format that will be stored in:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.sparse import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 66 }, { "cell_type": "code", "collapsed": false, "input": [ "# dense matrix\n", "M = array([[10,0,0,0], [0,.3,0,0], [0,2,1.5,0], [-3,0,0,100]])\n", "print M" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 10. 0. 0. 0. ]\n", " [ 0. 0.3 0. 0. ]\n", " [ 0. 2. 1.5 0. ]\n", " [ -3. 0. 0. 100. ]]\n" ] } ], "prompt_number": 67 }, { "cell_type": "code", "collapsed": false, "input": [ "# convert from dense to sparse\n", "M_sparse = csr_matrix(M); # or csc_matrix(M)\n", "print M_sparse " ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " (0, 0)\t10.0\n", " (1, 1)\t0.3\n", " (2, 1)\t2.0\n", " (2, 2)\t1.5\n", " (3, 0)\t-3.0\n", " (3, 3)\t100.0\n" ] } ], "prompt_number": 75 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can go back to the dense form:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "print M_sparse\n", "print M_sparse.todense()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " (0, 0)\t10.0\n", " (1, 1)\t0.3\n", " (2, 1)\t2.0\n", " (2, 2)\t1.5\n", " (3, 0)\t-3.0\n", " (3, 3)\t100.0\n", "[[ 10. 0. 0. 0. ]\n", " [ 0. 0.3 0. 0. ]\n", " [ 0. 2. 1.5 0. ]\n", " [ -3. 0. 0. 100. ]]\n" ] } ], "prompt_number": 69 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The easiest way to create a sparse matrix is using the lil_matrix command. It creates an empty matrix and then we allogating non-zero entries by using matrix indexing. This way we avoid creating large dense matrices." ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = lil_matrix((4,4)) # empty 4x4 sparse matrix using a linked list format\n", "A[0,0] = 1\n", "A[1,1] = 3\n", "A[2,2] = A[2,1] = 1\n", "A[3,3] = A[3,0] = 1\n", "A" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 70, "text": [ "<4x4 sparse matrix of type ''\n", "\twith 6 stored elements in LInked List format>" ] } ], "prompt_number": 70 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also change the format based on the sparsity:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "A = csr_matrix(A); A # if the rows are sparse" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 71, "text": [ "<4x4 sparse matrix of type ''\n", "\twith 6 stored elements in Compressed Sparse Row format>" ] } ], "prompt_number": 71 }, { "cell_type": "code", "collapsed": false, "input": [ "A = csc_matrix(A); A # if the columns are sparse" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 72, "text": [ "<4x4 sparse matrix of type ''\n", "\twith 6 stored elements in Compressed Sparse Column format>" ] } ], "prompt_number": 72 }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Exercise 12\n", "\n", "Construct a $10000 \\times 10000$ lil_matrix $A$ and add values to it by setting all rows and columns $2000$ to $3000$ into random values. Also set random diagonal entries (you can do this without a loop!). Then use the spy function from matplotlib.pyplot to \"spy\" the matrix.\n", "\n", "Using a vector $b$ with $10000$ random values solve the system $A x = b$. \n", "To do so use the direcet solver linsolve from scipy.sparse.linalg.dsolve.\n", "Change to different sparse formats (CSR and CSC) and measure the time needed to solve the system. You can use the time package." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# type your solution here" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 73 }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Solution" ] }, { "cell_type": "code", "collapsed": false, "input": [ "#%load solutions/sparse_matrices.py" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 74 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References:\n", "\n", "#### NumPy:\n", "\n", "- [NumPy tutorial](http://www.scipy.org/Tentative_NumPy_Tutorial) \n", "- [NumPy tutorial video](http://showmedo.com/videotutorials/series?name=i9KuJuNcG) \n", "- [NumPy: lock 'n load](http://mentat.za.net/numpy/intro/intro.html)\n", "- [NumPy for Matlab users](http://www.scipy.org/NumPy_for_Matlab_Users)\n", "\n", "\n", "#### SciPy\n", "- [SciPy getting started](http://www.scipy.org/getting-started.html)\n", "- [SciPy tutorial](https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/index.html)\n", "\n", "#### More general\n", "\n", "- [General Python tutorial](http://docs.python.org/tutorial/)\n", "- [Dive Into Python](http://diveintopython.org/)\n", "- [Hans Petter Langtangen book](http://link.springer.com/book/10.1007/978-3-642-02475-7/page/1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other packages\n", "Besides NumPy and Scipy, there are many other useful Python packages for scientific computing. Here is a short list: \n", "\n", "- [matplotlib](http://www.matplotlib.org/) - plotting and visualization\n", "- [sympy](http://sympy.org/) - symbolic computation\n", "- [mpi4py](http://mpi4py.scipy.org/) - parallel computing\n", "- [petsc4py](http://code.google.com/p/petsc4py/), [pytrilinos](http://trilinos.sandia.gov/packages/pytrilinos/) - Python bindings for the \"big 2\" parallel scientific libraries\n", "- [pyCUDA](http://mathema.tician.de/software/pycuda), [pyOpenCL](http://mathema.tician.de/software/pyopencl) - GPGPU computing\n", "- [FENiCS](http://fenicsproject.org/), [FiPy](http://www.ctcms.nist.gov/fipy/), [PyClaw](http://numerics.kaust.edu.sa/pyclaw/) - solve complicated PDEs with very sophisticated numerical methods\n", "- [networkX](http://networkx.github.com/), [pygraphviz](http://networkx.lanl.gov/pygraphviz/) - graphs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Copyright 2015, Maruan Al-Shedivat, Yiannis Hadjimichael, ACM Student Chapter.*" ] } ], "metadata": {} } ] }