Skip to content

Instantly share code, notes, and snippets.

@jhale

jhale/03_numpy.ipynb

Created May 7, 2019 14:29
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jhale/007de2b2b37ac3de56e2b8030d9fa9c8 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jhale/007de2b2b37ac3de56e2b8030d9fa9c8 to your computer and use it in GitHub Desktop.
Download ZIP
Draft notebooks from Scientific Python UL Course
Display the source blob
Display the rendered blob
Raw
03_numpy.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Lecture 3\n",
"\n",
"### Numpy"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Covered topics:\n",
"\n",
"* `numpy` - a package for storing and working with N-dimensional arrays of data.\n",
"* Why Python's `list` object is not adequate.\n",
" * Dealing with high-dimensional data tricky.\n",
" * Does not have 'maths' notation, e.g. $x + y$\n",
" * Speed.\n",
"* Slicing.\n",
"* Examples:\n",
" * Mandelbrot set.\n",
" * Markov chains.\n",
" * Calculating $\\pi$ via random sampling."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "notes"
}
},
"source": [
"Working with numbers is at the core of most activities in scientific and engineering computations. The topic is so important that there entire libraries dedicated to working with numbers. `numpy` is the most important and widely used package for dealing with numerical data in Python."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"We have already seen Python's `list` data structure. We used it for storing a list (sometimes called an array) of data of arbitrary type. We also know how to *slice* the list to get out particular elements we are interested in, and *iterate* over the list.\n",
"\n",
"*Reminder:*"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0, 3, 4]\n",
"test\n",
"0\n",
"3\n",
"4\n",
"5\n",
"test\n"
]
}
],
"source": [
"mylist = [0, 3, 4, 5, 'test']\n",
"print(mylist[0:3])\n",
"print(mylist[-1])\n",
"for l in mylist:\n",
" print(l)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### 'Problems' with Python's `list`.\n",
"\n",
"1. Dealing with high-dimension data is tricky.\n",
"\n",
"\\begin{equation}\n",
"A = \\begin{bmatrix} \n",
"0 & 1 & 2 \\\\\n",
"3 & 4 & 5 \\\\\n",
"6 & 7 & 8\n",
"\\end{bmatrix}\n",
"\\end{equation}\n",
"\n",
"Idea: A list of lists! Let's put each row in its own `list` and then put the rows in another list."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0, 1, 2], [3, 4, 5], [6, 7, 8]]\n"
]
}
],
"source": [
"A = [[0, 1, 2], [3, 4, 5], [6, 7, 8]]\n",
"print(A)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Exercises*: \n",
"1. Select the first row.\n",
"2. Select the first number of the second row.\n",
"3. Print the first column. (Not so easy)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0, 1, 2]\n",
"3\n",
"[0, 1, 2]\n",
"0\n",
"[3, 4, 5]\n",
"3\n",
"[6, 7, 8]\n",
"6\n"
]
}
],
"source": [
"print(A[0])\n",
"print(A[1][0])\n",
"for row in A:\n",
" print(row)\n",
" print(row[0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"2. `list` does not conform to 'standard' mathematical conventions.\n",
"\n",
"$$ x = [1, 2, 3]$$\n",
"$$ y = [4, 5, 6]$$\n",
"$$ x + y = \\; ? $$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[1, 2, 3, 4, 5, 6]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = [1, 2, 3]\n",
"y = [4, 5, 6] \n",
"x + y"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Something better: `numpy`"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Search results for 'create array'\n",
"---------------------------------\n",
"numpy.array\n",
" Create an array.\n",
"numpy.memmap\n",
" Create a memory-map to an array stored in a *binary* file on disk.\n",
"numpy.diagflat\n",
" Create a two-dimensional array with the flattened input as a diagonal.\n",
"numpy.fromiter\n",
" Create a new 1-dimensional array from an iterable object.\n",
"numpy.partition\n",
" Return a partitioned copy of an array.\n",
"numpy.ctypeslib.as_array\n",
" Create a numpy array from a ctypes array or POINTER.\n",
"numpy.ma.diagflat\n",
" Create a two-dimensional array with the flattened input as a diagonal.\n",
"numpy.ma.make_mask\n",
" Create a boolean mask from an array.\n",
"numpy.ctypeslib.as_ctypes\n",
" Create and return a ctypes object from a numpy array. Actually\n",
"numpy.asarray\n",
" Convert the input to an array.\n",
"numpy.ndarray\n",
" ndarray(shape, dtype=float, buffer=None, offset=0,\n",
"numpy.recarray\n",
" Construct an ndarray that allows field access using attributes.\n",
"numpy.chararray\n",
" chararray(shape, itemsize=1, unicode=False, buffer=None, offset=0,\n",
"numpy.pad\n",
" Pads an array.\n",
"numpy.asanyarray\n",
" Convert the input to an ndarray, but pass ndarray subclasses through.\n",
"numpy.copy\n",
" Return an array copy of the given object.\n",
"numpy.diag\n",
" Extract a diagonal or construct a diagonal array.\n",
"numpy.load\n",
" Load arrays or pickled objects from ``.npy``, ``.npz`` or pickled files.\n",
"numpy.sort\n",
" Return a sorted copy of an array.\n",
"numpy.array_equiv\n",
" Returns True if input arrays are shape consistent and all elements equal.\n",
"numpy.dtype\n",
" Create a data type object.\n",
"numpy.ufunc\n",
" Functions that operate element by element on whole arrays.\n",
"numpy.choose\n",
" Construct an array from an index array and a set of arrays to choose from.\n",
"numpy.nditer\n",
" Efficient multi-dimensional iterator object to iterate over arrays.\n",
"numpy.swapaxes\n",
" Interchange two axes of an array.\n",
"numpy.full_like\n",
" Return a full array with the same shape and type as a given array.\n",
"numpy.ones_like\n",
" Return an array of ones with the same shape and type as a given array.\n",
"numpy.empty_like\n",
" Return a new array with the same shape and type as a given array.\n",
"numpy.zeros_like\n",
" Return an array of zeros with the same shape and type as a given array.\n",
"numpy.asarray_chkfinite\n",
" Convert the input to an array, checking for NaNs or Infs.\n",
"numpy.diag_indices\n",
" Return the indices to access the main diagonal of an array.\n",
"numpy.nested_iters\n",
" Create nditers for use in nested loops\n",
"numpy.chararray.tolist\n",
" a.tolist()\n",
"numpy.put_along_axis\n",
" Put values into the destination array by matching 1d index and data slices.\n",
"numpy.ma.choose\n",
" Use an index array to construct a new array from a set of choices.\n",
"numpy.savez_compressed\n",
" Save several arrays into a single file in compressed ``.npz`` format.\n",
"numpy.matlib.rand\n",
" Return a matrix of random values with given shape.\n",
"numpy.datetime_as_string\n",
" Convert an array of datetimes into an array of strings.\n",
"numpy.ma.empty_like\n",
" Return a new array with the same shape and type as a given array.\n",
"numpy.ma.make_mask_none\n",
" Return a boolean mask of the given shape, filled with False.\n",
"numpy.around\n",
" Evenly round to the given number of decimals.\n",
"numpy.source\n",
" Print or write to a file the source code for a NumPy object.\n",
"numpy.diagonal\n",
" Return specified diagonals.\n",
"numpy.nan_to_num\n",
" Replace NaN with zero and infinity with large finite numbers.\n",
"numpy.einsum_path\n",
" Evaluates the lowest cost contraction order for an einsum expression by\n",
"numpy.histogram2d\n",
" Compute the bi-dimensional histogram of two data samples.\n",
"numpy.busdaycalendar\n",
" A business day calendar object that efficiently stores information"
]
}
],
"source": [
"import numpy\n",
"import numpy as np # np is now an alias (shorthand) for numpy\n",
"np.lookfor('create array')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x: [1 2 3]\n",
"y: [4 5 6]\n"
]
}
],
"source": [
"x = np.array([1, 2, 3])\n",
"y = np.array([4, 5, 6])\n",
"print(\"x: {}\".format(x))\n",
"print(\"y: {}\".format(y))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([5, 7, 9])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x + y"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"3. Speed"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Exercise*: Write a function `numpy_add` that takes two numpy arrays, loops over the entries, and computes the sum entry-by-entry and places the result into a third numpy array.\n",
"\n",
"Tips:\n",
"\n",
"* `def` - define a function\n",
"* `return x` - return `x` from the function\n",
"* `len(a)` - get the length of an array.\n",
"* `np.zeros(n)` - make a numpy array of length `n`\n",
"* `for i in range(0, len(a)): print(i)` - counting indices $0, ..., \\mathrm{len}(a) - 1$ "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"def list_add(a, b):\n",
" assert(len(a) == len(b))\n",
" c = np.zeros(len(a))\n",
" for i in range(0, len(a)):\n",
" c[i] = a[i] + b[i]\n",
" return c\n",
"\n",
"\n",
" \n",
"a = np.ones(1000000)\n",
"b = np.ones(1000000)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 680 ms, sys: 0 ns, total: 680 ms\n",
"Wall time: 676 ms\n"
]
}
],
"source": [
"%time c = list_add(a, b)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 10 ms, sys: 0 ns, total: 10 ms\n",
"Wall time: 8.94 ms\n"
]
}
],
"source": [
"%time c = a + b"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Key points:\n",
"* Avoid writing 'tight' loops in Python.\n",
" * Tight loop - looping over thousands of entries of data and performing a single operation on each piece of data.\n",
"* Think about how algorithms can be written without tight loops.\n",
"* Always prefer built-in functionality that can act on many pieces of data in one operation."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Example*: Plot the function $\\sin(x)$."
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {
"scrolled": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x: [0. 0.12822827 0.25645654 0.38468481 0.51291309 0.64114136\n",
" 0.76936963 0.8975979 1.02582617 1.15405444 1.28228272 1.41051099\n",
" 1.53873926 1.66696753 1.7951958 1.92342407 2.05165235 2.17988062\n",
" 2.30810889 2.43633716 2.56456543 2.6927937 2.82102197 2.94925025\n",
" 3.07747852 3.20570679 3.33393506 3.46216333 3.5903916 3.71861988\n",
" 3.84684815 3.97507642 4.10330469 4.23153296 4.35976123 4.48798951\n",
" 4.61621778 4.74444605 4.87267432 5.00090259 5.12913086 5.25735913\n",
" 5.38558741 5.51381568 5.64204395 5.77027222 5.89850049 6.02672876\n",
" 6.15495704 6.28318531]\n",
"y: [ 0.00000000e+00 1.27877162e-01 2.53654584e-01 3.75267005e-01\n",
" 4.90717552e-01 5.98110530e-01 6.95682551e-01 7.81831482e-01\n",
" 8.55142763e-01 9.14412623e-01 9.58667853e-01 9.87181783e-01\n",
" 9.99486216e-01 9.95379113e-01 9.74927912e-01 9.38468422e-01\n",
" 8.86599306e-01 8.20172255e-01 7.40277997e-01 6.48228395e-01\n",
" 5.45534901e-01 4.33883739e-01 3.15108218e-01 1.91158629e-01\n",
" 6.40702200e-02 -6.40702200e-02 -1.91158629e-01 -3.15108218e-01\n",
" -4.33883739e-01 -5.45534901e-01 -6.48228395e-01 -7.40277997e-01\n",
" -8.20172255e-01 -8.86599306e-01 -9.38468422e-01 -9.74927912e-01\n",
" -9.95379113e-01 -9.99486216e-01 -9.87181783e-01 -9.58667853e-01\n",
" -9.14412623e-01 -8.55142763e-01 -7.81831482e-01 -6.95682551e-01\n",
" -5.98110530e-01 -4.90717552e-01 -3.75267005e-01 -2.53654584e-01\n",
" -1.27877162e-01 -2.44929360e-16]\n"
]
}
],
"source": [
"x = np.linspace(0.0, 2.0*np.pi, num=50)\n",
"y = np.sin(x) # Apply sin function to all of the entries of x\n",
"print(\"x: {}\".format(x))\n",
"print(\"y: {}\".format(y))"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fa44f781e48>]"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEKCAYAAADenhiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8lPW17/HPyg3CNSAXIYCAIFqlgka6Fa0IWmztEWqt1Z621mo5dte22yq7cOruRXuhpW6tbY/VVndxt1vUShHrBUXEtlqVICA3uQgUQkBQDNeE3Nb5YyY4SWYgl8k8z8x8369XXpl55pnJGg2z8vv91m895u6IiIgkU07QAYiISOZRchERkaRTchERkaRTchERkaRTchERkaRTchERkaRTchERkaRTchERkaRTchERkaTLCzqAoPTp08eHDh0adBgiImll2bJl77p73+Odl7XJZejQoZSWlgYdhohIWjGzf7bkPE2LiYhI0im5iIhI0im5iIhI0im5iIhI0im5iIhI0oUmuZjZg2a228xWJ3jczOweM9tkZm+a2Vkxj11rZhujX9emLmpJZP7yHYyftZhhM55i/KzFzF++I+iQRCSFwlSK/HvgV8BDCR7/ODAy+vUR4F7gI2bWG/geUAI4sMzMFrj7+x0ecZabv3wHsxeup7yikoFFhUyfPIqpY4uZv3wHM+etorKmDoAdFZXMnLfq6PPiPUdEMktokou7/9XMhh7jlCnAQx65LvOrZlZkZgOACcDz7r4XwMyeBy4FHu7YiLNbogTi7sx69q2jxxtU1tRx2/xV1NQ5R2rrGz0HUIIRyTChSS4tUAxsj7lfFj2W6HgzZjYNmAYwZMiQjokyS8xeuD5uArnlsZXUe/znHDxS1+xYZU0dP1v41tHX1IhGJDOkU3KxOMf8GMebH3S/H7gfoKSkJMFHoDQVb/qrvKIy7rn1Dj0L89lXWdPi1y+vqOKWx1ZSF81KGtGIpL/QLOi3QBkwOOb+IKD8GMclCRqmv3ZUVOJEPvhveWxl/OwNFBcV8oPLT6cwP7fR8cL8XHp1yY/7HIOjiaVBZU0dsxeub/8bEJFApNPIZQFwk5nNJbKgv8/dd5rZQuDHZtYret7HgJlBBZlp4k1/1dU7BbmGmR1dP4FIAomdzmo62gEardM0PKfp6zcor6hMWDQgIuEWmuRiZg8TWZzvY2ZlRCrA8gHc/TfA08AngE3AYeC66GN7zewOYGn0pW5vWNyX9ks0/VVT59z12TMTfvBPHVucMAk0fc7shevZEefnOGi6TCRNWaT4KvuUlJS4uiIf2/a9h5l050tU19U3e6y4qJCXZ0xMys9pWnkG0Dk/B3cajYw64meLSOuY2TJ3LzneeaEZuUjwYqegirrkc/hILQbk5xo1dR/8EdIw/ZUsiabRbn5kRdzzE42mRCQ8lFwEaD56eP9wDTkG37nsNE7o2qnD1z3iTaMlmi7r36NzUn+2iCSfkosA8Rfu6x0e/PtWXp4xMZA1jumTRzWbLgM4XF3LL1/YyNyl27XQLxJS6VSKLB0o0VRTkFNQU8cW85MrRlNcVIgRWWu5dfIp5Jhx5/MbGpVHz5y3Sv3LREJEIxdh9Y59mEG82o6BRYWpDyhGvOmyP7y6DZps0mzYF6PRi0g4aOSS5Uq37uWa375Kj875dMpr/OuQ7IX7ZHlnX1Xc41roFwkPjVyyUENV2I6KSgw4oVsBT9x0Pku37E2LDYsDiwrjLvQHPcoSkQ8ouWSZplVhDhysqmXplr3H3PgYJokW+s8Z2ivBM0Qk1TQtlmXiVYVV1danVR+vpgv9A4s68+FBPZi/opy5r28LOjwRQSOXrBPGqrC2aDrKqq6t5ysPlTLzz6tYu3M/L6zbHfrpPZFMppFLFtlz4Ag5OfGuUJD+6xUFeTn85vNnM+yErjz0j3+qTFkkYEouWaKqpo6vPFSK4WlTFdZahQXxOyyrfb9I6mlaLAvU1zu3PLqSlWUV3Pu/z6Yq+mGbidNGu1SmLBIKSi4ZLLbkGODyMwdw6RknApnbsl5lyiLhoGmxDBV7BckGz699J+PXHqZPHtXsKpi5OZYR034i6SQ0ycXMLjWz9Wa2ycxmxHn8LjNbEf3aYGYVMY/VxTy2ILWRh1O8kuPKmvQqOW6LpmXK3TrlUlfv1GfpdYtEghKKaTEzywV+DVwClAFLzWyBu69tOMfdb445/+vA2JiXqHT3MamKNx1kSslxW8SWKdfW1fP5B15j5rxVnNK/O2cU9ww4OpHsEIrkAowDNrn7ZgAzmwtMAdYmOP8aIpdBljjcnU75OVTVNL+KY7atPeTl5vCrz53FJ+/5O1988DU65eWya19VxhUyiIRNWKbFioHtMffLoseaMbOTgGHA4pjDnc2s1MxeNbOpHRdmevif17dRVVNPXpM9LZlSctxafbp14upxg9l7qIad+6q0/0UkBcKSXOLt7Es0SX418Cd3j11QGBK9pvPngLvN7OS4P8RsWjQJle7Zs6d9EYfU+l0HuP3JtVwwsg+zP/3hRtdC+ckVo7P2L/XHSsuaHdP+F5GOE5ZpsTJgcMz9QUB5gnOvBr4We8Ddy6PfN5vZEiLrMW83faK73w/cD1BSUpJxK7xVNXV8/eE36N45jzuvOpN+3TvzqbMHBR1WKGTzGpRIEMKSXJYCI81sGLCDSAL5XNOTzGwU0Av4R8yxXsBhdz9iZn2A8cDPUhJ1CDTsZSmvqKSwIJfD1XXM+fI4+nXXdeZjaf+LSGqFYlrM3WuBm4CFwDrgUXdfY2a3m9nlMadeA8x1b1RXehpQamYrgReBWbFVZpksdi+LA4er68jLMd4/VB10aKETb/9Lnva/iHSYsIxccPengaebHPtuk/vfj/O8V4DRHRpcSMXby1Jb77rcbxwN/z0aRnldCnI5VF3HgJ4a4Yl0hNAkF2k9rSO0Tuz+l0NHavn4L/7GrX9ayTPf/CjdOumfgkgyhWJaTNom0XqB1hGOr2unSNFD2fuV/OipdUGHI5JxlFzS2I0Thjc7lq17WdrinKG9mXbBcB5+fRtn3/E8w2Y8xfhZi7X3RSQJlFzS2Otb3ifXoF/3TtrL0kYj+nXDgPcOVWtzpUgSaaI5Tb2w7h2eXFnOty45hW9MGhl0OGnr7kUbm+3WbdhcqSQt0nYauaShA1U13DZ/NaP6d+fGC+M2I5AWUlGESMdQcklDP3t2Pbv2VzHr06MpyNP/wvZQUYRIx9C0WJqI3YnvwIUj+zB2SK+gw0p70yePYua8VY32C2lzpUj76c/eNNB0Jz7Aa1v3atE5CZpeXKwwPwd35ywlbpF2UXJJA/F24ldlwVUlU2Xq2GJenjGRLbMu48VbL6Jzfi7fW7Aa19UrRdpMySUNaNE5dU7s2ZmbLzmFF9fvYeGaXUGHI5K2lFzSgBadU+tL5w3l1BO784Mn13LoSG3Q4YikJSWXNPC/zhzQ7Jh24necvNwcfvSpM9i5r4qP/PgF7dwXaQNVi4VcTV09z699hz5d8ynIy2Wnrv+eEtv3VpKbYxyMjlwadu4D+u8u0gJKLiE355WtvL3nEA9cW8Kk0/oHHU7WmL1wPXX1jRf0tXNfpOU0LRZiew4c4ReLNjJhVF8mntov6HCyioooRNonNMnFzC41s/VmtsnMZsR5/EtmtsfMVkS/boh57Foz2xj9uja1kXecn0dLkP/jkx/CzIIOJ6uoiEKkfUKRXMwsF/g18HHgQ8A1ZvahOKc+4u5jol+/iz63N/A94CPAOOB7Zpb2O+BWbq/g0WXb+fL5wzi5b7egw8k6uiyySPuEZc1lHLDJ3TcDmNlcYAqwtgXPnQw87+57o899HrgUeLiDYu1Q85fv4GcL36K8ooocg+F9ugQdUlZqelnkTnmRnfvnj+wTcGQi6SEUIxegGNgec78seqypT5vZm2b2JzMb3Mrnhl5Dm5fyiioA6h1+8OQ6lcAGJHbn/lPfvIA6h7ue3xB0WCJpISzJJd6CQtPeG08CQ939w8AiYE4rnhs50WyamZWaWemePXvaHGxHidfmpaFCSYJ1ct9ufP5fTuLh17ex4Z0DQYcjEnphSS5lwOCY+4OA8tgT3P09dz8Svftb4OyWPjfmNe539xJ3L+nbt29SAk8mVSiF2zcnjaRbpzx+9NS6oEMRCb2wJJelwEgzG2ZmBcDVwILYE8wsdpv65UDDv/CFwMfMrFd0If9j0WNpp3/PznGPq0IpHHp1LeAbk0by0oY9LFm/O+hwREItFAv67l5rZjcRSQq5wIPuvsbMbgdK3X0B8A0zuxyoBfYCX4o+d6+Z3UEkQQHc3rC4n25G9uvGrn1VjY6pzUu4fOHck7h3ySZumFNKXb2rW4JIAqFILgDu/jTwdJNj3425PROYmeC5DwIPdmiAHWzru4f4x9vvce7Jvdn2XiXlFZX64AqhZ1bt4kBVLbXR3ftqCyMSX2iSS7ab/dx68nNz+MVnx9KvR/zpMQne7IXrqa5TWxiR4wnLmktWW7m9gqfe3MlXLhimxBJyKroQaRkll4C5Oz95Zh0ndC1g2oUnBx2OHIfawoi0jJJLwJas38Orm/fyzYsjZa4SbvHawuSqLYxIM/o0C0hsm5fcHKNrQe7xnySBa9oWprAglyM1dZQMTft2diJJpeQSgIY2Lw278evqndvmryE3J0eLwmlg6tjio/+fdu6r5MLZS7h70UZ+/pkzA45MJDw0LRYAtXnJHAN6FvLFfzmJeW+UsVFtYUSOUnIJgCqOMsu/XjSCLgV5/Pw5/XEg0kDJJQAnqs1LRundtYCvXDCchWveYcX2iqDDEQkFJZcAnDWkqNkxtXlJb9dfMIwTuhbws2ffCjoUkVDQgn6K7aus4W8b3+X0gd2pOFyrNi8ZolunPL520Qhu/8tazr7jefYeqtb/V8lqSi4p9sDfNrO/qpafXXkmpw/sGXQ4kkQ9OudhwHuHqgH1HZPspmmxFNp7qJoH/r6Fy0YPUGLJQHct2tjsKnWqApRspeSSQr956W0qa+q4+ZKRQYciHUBVgCIfUHJJkd37q5jzylamjilmRL/uQYcjHUB9x0Q+oOSSIr9+cRN19c43L9aoJVPF6ztWmJ+jKkDJSqFZ0DezS4FfELkS5e/cfVaTx78F3EDkSpR7gC+7+z+jj9UBq6KnbnP3y1MW+DHMX77jaA8qB84d3puTTugadFjSQWL7ju2IToVdefYgLeZLVgpFcjGzXODXwCVAGbDUzBa4+9qY05YDJe5+2My+CvwM+Gz0sUp3H5PSoI+jaf8wgDe2VTB/+Q592GSwhr5j9fXOJ+75Gy9veo/aunrycjVJINklLL/x44BN7r7Z3auBucCU2BPc/UV3Pxy9+yowKMUxtkq8/mFHautVOZQlcnKMf7v4FDa/e4gFK8uDDkck5cKSXIqB7TH3y6LHErkeeCbmfmczKzWzV81saqInmdm06Hmle/bsaV/Ex6HKIZl8en8+NKAHv3hhI7V19UGHI5JSYUkuFudY0y0DkRPNPg+UALNjDg9x9xLgc8DdZhb3ko7ufr+7l7h7Sd++fdsb8zGpckjMjJsvOYV/vneYect3BB2OSEqFJbmUAYNj7g8Cms0lmNnFwHeAy939SMNxdy+Pft8MLAHGdmSwLTF98ihyrXHOVP+w7HPxaf0YXdyTXy7eSI1GL5JFwpJclgIjzWyYmRUAVwMLYk8ws7HAfUQSy+6Y473MrFP0dh9gPBBbCBCIMYOLqHena6dcDCguKuQnV4zWYn6WMTO+dckpbN9byTk/XMSwGU8xftZi5mskIxkuFNVi7l5rZjcBC4mUIj/o7mvM7Hag1N0XEJkG6wY8ZpERQUPJ8WnAfWZWTyRZzmpSZRaIXy7eREFeDi/eOoF+3eO32JfsUHG4GjOoqKwB1HNMskMokguAuz8NPN3k2Hdjbl+c4HmvAKM7NrrW2fruIeav2MGXzhuqxCL8/LkNeJMVxIaeY0oukqnCMi2WUX714ibycoz/c+HwoEOREFDloGQjJZck2/ruIf68fAef/5eTNGoRQJWDkp2UXJJMoxZpKl7Psc7qOSYZLjRrLumsaQ+xC0/po1GLHBWv59iUMQO13iIZTcmlneL1EHtt8171EJNGGnqOuTtTfv0yr7z9HjV19eSr55hkKP1mt1O8HmJV6iEmCZgZ35g4ku17K7XXRTJaq5OLmXWNdjEWVAkkrTfptH6cPrAHv35xk3qOScY6bnIxsxwz+5yZPWVmu4G3gJ1mtsbMZptZVl/9SpVA0lpmxjcmjWTre4fVMVkyVktGLi8CJwMzgRPdfbC79wMuINL6fla0mWRWmj55FLk56iEmrfOxD/XntAE9+NXiyBVKRTJNS5LLxe5+h7u/6e5Hx/DuvtfdH3f3TwOPdFyI4XbOsN64O10L1ENMWi6y9jKCze8e4i9vavQimee41WLuXgNgZncDN7s3bWTxwTnZ6DdL3iY3x3juWxdSrKkwaYXJp5/IiT06ccujK/m3uSsYWFTI9Mmj9IeJZITWLOgfBBaYWVcAM/uYmb3cMWGlh137qnhk6XauPHuwEou02oKV5bx3qJraesf5oKGlqsgkE7Q4ubj7bcDDwBIz+ztwCzCjowJLB/f99W3q3PnXCXGvTSZyTLMXrqemrvFEQENDS5F01+JNlGY2CfgKcAgYAFzv7ln7r2D3gSr+57VtXDG2mMG9uwQdjqQhlbFLJmvNtNh3gP9w9wnAlcAjZjaxQ6JKA7/962Zq6ur52kUjgg5F0pTK2CWTtWZabKK7/z16exXwceCHyQrEzC41s/VmtsnMmk23mVknM3sk+vhrZjY05rGZ0ePrzWxysmJK5L2DR/jDq9uYMqaYoX26dvSPkwwVr6FlQa4aWkpmOO60mJlZggqxndGpsoTntFR0x/+vgUuAMmCpmS1ockXJ64H33X2EmV0N/BT4rJl9iMhlkU8HBgKLzOwUd2/ckyUJGhpUNjQfPPXE7sn+EZJFYhtalldUkpNj9O/RiSljBgYcmWSi2Aa7qahMbNEmSjP7upkNiT0Yvdb9uWY2B7i2nXGMAza5+2Z3rwbmAlOanDMFmBO9/SdgkkWudzwFmOvuR9x9C7Ap+npJ1dCgckfMfPjdizaqskfaZerYYl6eMZEtsy7jx586g+3vV7Jkw56gw5IME/v5larKxJYkl0uBOuBhMys3s7VmthnYCFwD3OXuv29nHMXA9pj7ZdFjcc9x91pgH3BCC5/bbvEaVKqyR5LpU2MHUVxUyD0vbKQdEwEizQTx+XXc5OLuVe7+/9x9PHASMAk4y91PcvevuPuKJMRhcY41/deV6JyWPDfyAmbTzKzUzEr37GndX4eq7JGOVpCXw1cnnMzybRW8vOm9oMORDBLE51eruiK7e42773T3iiTHUQYMjrk/CGjaE+PoOWaWB/QE9rbwuQC4+/3uXuLuJX379m1VgKrskVT4TMkgTuzRmXsWbww6FMkgQXx+tTi5RKu1Pmdm/9fMvtvwlaQ4lgIjzWxYdC3namBBk3MW8MHazpXA4mgRwQLg6mh8w4CRwOtJiuuoeJU9alApydYpL5f/c+FwXt+yl1c3a/QiyTF98ihyLbUNdltzJconiKxzLAOOJDMId681s5uAhUAu8KC7rzGz24FSd18APAD8t5ltIjJiuTr63DVm9iiwFqgFvtYRlWJNK3vUB0o6yjXjhnDnc+u59sHXqa6t1++atNtZQ3pR707XTrkcPlKXkt8pa+nCoZmtdvczOiySFCspKfHS0tKgwxBpZv7yHdz62EpqY1rxF+bnqtu2tNmMx99k3vId/O3fL6J/j87tei0zW+buJcc7rzVrLq+Y2eh2xCQiLTB74fpGiQVUmShtV/b+YR5/o4yrzxnc7sTSGq2ZFjsfuC5ahnyESJWWu/uHOyQykSylykRJpt+89DYAN16Y2ga7rUkulxJNKB0Ui4gQqeDZESeRqDJRWmvnvkoeXVrGlWcPTvnvz3GnxaLt9QHWAKuA1dGvNdHvIpJE8SoTO+Wp55i03n0vbaY+oMuCtORKlOdHv6uRlkgKNK1MBBjRr5sW86VVdu+v4uHXt3HFWcFcFqQ113P5DPCsux8ws9uAs4A73H15h0UnkqWmji0+mkx+sWgjdy3awOod+zijuGfAkUnYNW2we0r/YMYFrakW+49oYjkfmEykieRvOiYsEWnwpfFD6d45j19q174cR7wGu3c+tyGQBrutSS4NGxMvA+519yeAguSHJCKxehbmc934YSxc8w7rdu4POhwJsTA12G1NctlhZvcBVwFPm1mnVj5fRNro+vHD6NZJoxc5tjCVsbcmOVxFpD3LpdHGlb2B6R0SlYg00rNLPl86byhPr9rF+l0Hgg5HQipMDXZbc5njw+4+z903Ru/vdPfnOi40EYl1/fnD6FqQq9GLJPTVCcObHQuqwW5rNlGKSIB6dS3gi+cN5d4lb/P6lkXsOXBETS2lkbL3qwDo171T4L8fSi4iaaQ4Or2x+0CkMXnD5WoBJZgst/dQNQ/9YyuXnzmQe64ZG3Q4WpAXSSf3Lnm72TE1tRSA3/5tM5U1dXxj0oigQwGUXETSSpiqgSQ89h6qZs4rW/lfHx7IiH7haKai5CKSRsJUDSThEbZRC4QguZhZbzN73sw2Rr/3inPOGDP7h5mtMbM3zeyzMY/93sy2mNmK6NeY1L4DkdTR5balqTCOWiAcC/ozgBfcfZaZzYje/3aTcw4DX3T3jWY2EFhmZguj+20Aprv7n1IYs0ggYptaNrT4+PqkEVrMz0JNe4idPjA8iQXCkVymABOit+cAS2iSXNx9Q8ztcjPbDfQFKhDJMg1NLSsOV3P+T19kVdm+oEOSFGvoIRbb6uXuRZvo36MwNH9oBD4tBvR3950Q2ZgJ9DvWyWY2jkhPs9iymR9Fp8vuiralEcl4RV0K+PL4oTyzehdry9VzLJuEqYdYIilJLma2yMxWx/ma0srXGQD8N3Cdu9dHD88ETgXOIdKSpumUWuzzp5lZqZmV7tmzp43vRiQ8rj9/ON0753H3og3HP1kyRjpUDaYkubj7xe5+RpyvJ4B3okmjIXnsjvcaZtYDeAq4zd1fjXntnR5xBPgvYNwx4rjf3UvcvaRv377JfIsigejZJZ8bzh/Oc2vfYfUOTY9li3SoGgzDtNgC4Nro7WuBJ5qeYGYFwJ+Bh9z9sSaPNSQmA6aiSy9Llrnu/KH00Oglq3zlo8OaHQtb1WAYksss4BIz2whcEr2PmZWY2e+i51wFfBT4UpyS4z+a2SpgFdAH+GFqwxcJVo/O+Uz76HAWrdvNOT9cxLAZTzF+1uJALhAlqbFp90FyDPr36IQRaQv0kytGh2YxH0JQLebu7wGT4hwvBW6I3v4D8IcEz5/YoQGKpIE+3SJ1LHsOqudYptu+9zCPLN3O5z4yhB9OHR10OAmFYeQiIu30y8Wbmh0LW/WQJMc9L2wkx4ybLhoZdCjHpOQikgHSoXpI2u/tPQd5/I0yvvAvJ3Fiz85Bh3NMSi4iGSAdqoek/e5etJHO+bncOOHkoEM5LiUXkQwQv+dYTqiqh6R91u3cz5Mry7lu/NCja2xhFviCvoi0X7yeY1efM1iL+RkgtodYQ2VYOlByEckQDT3H6uqdS+/+Ky9tfJfaunrycjVBka6a9hBz4I6/rKNLQV7o/3DQb51IhsnNMW752Cg27znEvDe01yWdpUMPsUSUXEQy0OTT+3Pm4CLuWrSBqiYfTpI+0rkKUMlFJAOZGd+ePIqd+6r4w6v/DDocaaMBRfHLjdOhClDJRSRDnTeiD+eP6MP/W/I2B6pqgg5H2mDSaf2bHQtbD7FElFxEMtj0yaPYe6ia82YtVs+xNFNdW89L6/cwoEcnBhZ1Dm0PsURULSaSwba8e4gcgwNVtYB6jqWTP772T7btPcycL4/jwlPS7xIhGrmIZLDZC9dT742PpUu1UTbbX1XDPS9sZPyIE/joyD5Bh9MmSi4iGSydq42y2X0vvc37h2uYcelpRC5VlX6UXEQymHqOpZ9d+6p44O9buPzMgYwe1DPocNpMyUUkg8XrOVaQp55jYXbX8xuoq/e0/38U+IK+mfUGHgGGAluBq9z9/Tjn1RG52iTANne/PHp8GDAX6A28AXzB3as7PnKR8IvtOVZeUUmOQb9unbj8zIEBRyZNzV++gx8/vY7dB47QtVMuy/75PoN7dwk6rDYLw8hlBvCCu48EXojej6fS3cdEvy6POf5T4K7o898Hru/YcEXSy9Sxxbw8YyJbZl3GnVeNoayikvkrVI4cJg09xHYfiFxJ9NCROmbOW5XWZeNhSC5TgDnR23OAqS19okVWuiYCf2rL80WyzeVnDuTMQT352bPrqaxWW5iwSOceYomEIbn0d/edANHv/RKc19nMSs3sVTNrSCAnABXuXhu9XwaoeF8kgZwc47ZPfohd+6v47d82Bx2ORGViVV9K1lzMbBFwYpyHvtOKlxni7uVmNhxYbGargP1xzvM4xxrimAZMAxgyZEgrfrRI5jhnaG8+MfpE7l3yNp89ZzD9e4T7crnZoEdhPvsqm7foSeeqvpQkF3e/ONFjZvaOmQ1w951mNgDYneA1yqPfN5vZEmAs8DhQZGZ50dHLIKD8GHHcD9wPUFJSkjAJiWS6b196KgtX7+Kiny+hsrqOgUWFTJ88Srv2A7D3UDXVtXXkGI02vKZLD7FEwjAttgC4Nnr7WuCJpieYWS8z6xS93QcYD6x1dwdeBK481vNFpLHl2yowMw5X1+F80BYmnReQ09Wdz62nui5SelxcVJh2PcQSCbwUGZgFPGpm1wPbgM8AmFkJcKO73wCcBtxnZvVEEuIsd18bff63gblm9kNgOfBAqt+ASLqZvXA9tU36wjQsIKfzB1q6WVO+j4df38YXzx3KVyeM4KsTRgQdUtIEnlzc/T1gUpzjpcAN0duvAKMTPH8zMK4jYxTJNJm4gJxu3J0fLFhLUZcCbr74lKDDSbowTIuJSIqpLUzwnnxzJ69v3cv0yaPo2SU/6HCSLvCRi4ik3vTJo5g5b1WjvRW5OZbWC8jpYP7yHUe7JZjBoKLOXFUyOOiwOoRGLiJZaOrYYn5fDDwBAAANUUlEQVRyxeijC8hdC3Kpq3dG9OsWdGgZq2EX/o6KSpxIZdjug9U8uTJhgWtas0jBVfYpKSnx0tLSoMMQCYV9lTVMuvMlinsVMu+r55Gbk55t3sNs/KzF7IizplVcVMjLMyYGEFHbmNkydy853nkauYgIPQvzue2y01i5vYK5S7cFHU5GyrYiCiUXEQFgypiBnDv8BH76zFu8e/BI0OFknGwrolByEREAzIw7pp7OwSO1XPDTFxk24ynGz1qsjZVJcuOE4c2Opfsu/GNRchGRo1bv2E+OGZU12rmfbK9veZ9cg37dO2XMLvxjUSmyiBylnfsd44V17/DkynK+dckpfGPSyKDDSQmNXETkqGxbdE6FA1U13DZ/NaP6d+fGC08OOpyUUXIRkaOybdE5FX767Fvs2l/FrE+PpiAvez5yNS0mIkfF27lvwL9dnB1TOckSuxPfgQtH9mHskF5Bh5VS2ZNGReS4mu7c79OtAAfW7TwQdGhpo+lOfIDXtu7NuqIIjVxEpJGpY4sbLd5/74nVPPjyFi75UH/OPfmEACNLD7MXrm808gOoqqnPuqIIjVxE5Ji+/fFTGXpCF259bCUHqppfilcaU1FEhEYuInJMXQryuPOqMXzmN6/wlTmlbH+/kvKKSl0aOYH+PTuza19Vs+PZVhQR+MjFzHqb2fNmtjH6vdmql5ldZGYrYr6qzGxq9LHfm9mWmMfGpP5diGS2s0/qxUWn9uPVLXuPriVog2Vz9fVOr8Lm12bJ5J34iQSeXIAZwAvuPhJ4IXq/EXd/0d3HuPsYYCJwGHgu5pTpDY+7+4qURC2SZdaV7292rGGDpUT81ytbWbfrAFecVXy0KCLTd+InEoZpsSnAhOjtOcAS4NvHOP9K4Bl3P9yxYYlIrJ1xpnog+9YSElm9Yx+znlnHxaf1587PnIlZdl+2IAzJpb+77wRw951m1u84518N/GeTYz8ys+8SHfm4e9yWrmY2DZgGMGTIkPZFLZJlBhYVxr0eSbatJcSK3c+Sk2N0Lchl9pUfzvrEAimaFjOzRWa2Os7XlFa+zgBgNLAw5vBM4FTgHKA3xxj1uPv97l7i7iV9+/ZtwzsRyV7TJ4+iMD+30bGC3JysW0to0HQ/S129U1Vbz0sb9gQdWiikZOTi7hcneszM3jGzAdFRywBg9zFe6irgz+5+tB6yYdQDHDGz/wJuTUrQItJIw5pBw1/qebmG4Zw5uCjgyIIRbz9LdW327WdJJAwL+guAa6O3rwWeOMa51wAPxx6IJiQsMg6dCqzugBhFhEiCeXnGRLbMuowXb51Al0553PjfyzhcXRt0aCmn/SzHFoY1l1nAo2Z2PbAN+AyAmZUAN7r7DdH7Q4HBwEtNnv9HM+tLpAXSCuDG1IQtkt0G9erCPdeM5doHX+d///Y1dh+ooryiKmv2v/Tr0Yl39jdf3s3mNahYgScXd38PmBTneClwQ8z9rUCz31Z3n9iR8YlIYheM7MsnRp/IX97cdfRYw/4XIGMTzIGqmrjTPtm4nyWRMEyLiUgaW76totmxTN7/UltXz03/s5zdB6u58cLhWb+fJZHARy4ikt7KK7Jn/4u78/0n1/DShj385IrRXDNuCDM+flrQYYWSkouItEs27H9p2M/S8D4njurLNeO0V+5YNC0mIu0Sb/8LwHXjh6Y+mA4Qu5+lwT82v6eeaseh5CIi7dL0AmP9uneiMD+H37+ylZ370n9qLN5+lsro9VkkMU2LiUi7Nb3A2JtlFXzut68x5Vd/Jycnh3f2pW+JcrwpP8jMNaVk0shFRJLuw4OKuG78UHYfqGbXvqq0bdH/+LKyhI9l0ppSR1ByEZEOMe+N5kkknUqUH/rHVm55bCUj+3Wlc37jj0rtZzk+TYuJSIdIx/YoTavCTh/Yg8e/eh7Prt51tKdauk7vpZqSi4h0iEQlykVdml+pMQwiVWFvUllTf/TY23sO8uzqXc3WlOT4NC0mIh0iXolyjsH7h2u4Yc5Szpv1AsNmPMX4WYtDsQ4z65m3GiUWgCpVhbWZRi4i0iGatugfWFTIzReP5PE3yli07oMra4ShF9nza99h1/7s6TSQCkouItJh4k0n3bVoQ7PzGhb6U5FcYq8eeWLPzgw9oQv/2LyXvByjtt6bna+qsLZRchGRlDpWL7LYD/6OWDhv2G3fsCly574qdu6rYuKpffnEGQP4jyfWNNowqaqwttOai4ikVKKRgAO3Prby6GWDO2JfTLzd9gDrdx3kypLBjToNqMtx+wQ+cjGzzwDfB04DxkWv4xLvvEuBXwC5wO/cfVb0+DBgLtAbeAP4grtXpyB0EWmD6ZNHNRo9AHTOy6HOnZq6xtNS7ZkuazoKuvLs4uPutldVWPIEnlyIXJb4CuC+RCeYWS7wa+ASoAxYamYL3H0t8FPgLnefa2a/Aa4H7u34sEWkLeIt9E+fPIqbH1kR9/wdFZVU1dQl3GsSbyoNaJTAdlRU8osXNmFERkhNaV0l+cw93n/q1DOzJcCt8UYuZnYu8H13nxy9PzP60CxgD3Ciu9c2Pe9YSkpKvLQ07iBJRAIwftbihCOLglyjrt6JHdh0zs/h02OLmbe8vNEoKD/XyDHjSG19s9fpWZhHda03W1fR9FfLmdkydy853nnpsuZSDGyPuV8WPXYCUOHutU2Oi0iaibcvpnN+DjdeOJy83ByazJhRVVPPH1/f3mwNpabO4yYWgP2VtVpXSZGUTIuZ2SLgxDgPfcfdn2jJS8Q55sc4niiOacA0gCFDdKEfkTBJNF02dWwx9720OSk/Y2BRodZVUiQlycXdL27nS5QBg2PuDwLKgXeBIjPLi45eGo4niuN+4H6ITIu1MyYRSbJEH/yJWsnkmlEXZ2q/qDCfI7X1KisOULpMiy0FRprZMDMrAK4GFnhkwehF4MroedcCLRkJiUgaiTdlVpifyzUfGRz3+PcvP13TXwELvFrMzD4F/BLoCzxlZivcfbKZDSRScvyJ6GL9TcBCIqXID7r7muhLfBuYa2Y/BJYDDwTwNkSkAx1ryqzkpN4JN14qmQQnNNViqaZqMRGR1su0ajEREUkjSi4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0WVstZmZ7gH+28el9iGzgTFfpHj+k/3tI9/gh/d9DuscPwbyHk9y97/FOytrk0h5mVtqSUrywSvf4If3fQ7rHD+n/HtI9fgj3e9C0mIiIJJ2Si4iIJJ2SS9vcH3QA7ZTu8UP6v4d0jx/S/z2ke/wQ4vegNRcREUk6jVxERCTplFxaycwuNbP1ZrbJzGYEHU9rmNmDZrbbzFYHHUtbmNlgM3vRzNaZ2Roz+2bQMbWWmXU2s9fNbGX0Pfwg6JjawsxyzWy5mf0l6Fjawsy2mtkqM1thZmnXwdbMiszsT2b2VvTfw7lBx9SUpsVawcxygQ3AJUQuYLYUuMbd1wYaWAuZ2UeBg8BD7n5G0PG0lpkNAAa4+xtm1h1YBkxNl//+AGZmQFd3P2hm+cDfgW+6+6sBh9YqZvYtoATo4e6fDDqe1jKzrUCJu6flPhczmwP8zd1/F73GVRd3rwg6rlgaubTOOGCTu29292pgLjAl4JhazN3/CuwNOo62cved7v5G9PYBYB2QVhfs8IiD0bv50a+0+gvPzAYBlwG/CzqWbGRmPYCPEr12lbtXhy2xgJJLaxUD22Pul5FmH26ZwsyGAmOB14KNpPWiU0orgN3A8+6ebu/hbuDfgfqgA2kHB54zs2VmNi3oYFppOLAH+K/o1OTvzKxr0EE1peTSOhbnWFr91ZkJzKwb8Djwb+6+P+h4Wsvd69x9DDAIGGdmaTNFaWafBHa7+7KgY2mn8e5+FvBx4GvRKeN0kQecBdzr7mOBQ0Do1n+VXFqnDBgcc38QUB5QLFkpuk7xOPBHd58XdDztEZ3KWAJcGnAorTEeuDy6ZjEXmGhmfwg2pNZz9/Lo993An4lMeaeLMqAsZsT7JyLJJlSUXFpnKTDSzIZFF9GuBhYEHFPWiC6GPwCsc/f/DDqetjCzvmZWFL1dCFwMvBVsVC3n7jPdfZC7DyXy+7/Y3T8fcFitYmZdowUhRKeTPgakTQWlu+8CtpvZqOihSUDoilrygg4gnbh7rZndBCwEcoEH3X1NwGG1mJk9DEwA+phZGfA9d38g2KhaZTzwBWBVdM0C4P+6+9MBxtRaA4A50crDHOBRd0/Lct401h/4c+RvFfKA/3H3Z4MNqdW+Dvwx+kfuZuC6gONpRqXIIiKSdJoWExGRpFNyERGRpFNyERGRpFNyERGRpFNyERGRpFNyERGRpFNyERGRpFNyEQmJ6LVqLone/qGZ3RN0TCJtpR36IuHxPeB2M+tHpOPz5QHHI9Jm2qEvEiJm9hLQDZgQvWaNSFrStJhISJjZaCK9x44osUi6U3IRCYHoJZz/SOTKpofMbHLAIYm0i5KLSMDMrAswD7jF3dcBdwDfDzQokXbSmouIiCSdRi4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0Si4iIpJ0/x8JXCmB/OD4lwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"plt.xlabel(r\"$x$\")\n",
"plt.ylabel(r\"$\\sin(x)$\")\n",
"plt.plot(x, y, 'o-')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Slicing"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Slicing* allows us to select one or more than one elements of an existing array."
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0 1 2 3 4]\n"
]
}
],
"source": [
"import numpy as np\n",
"a = np.arange(0, 5)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Selecting single elements of an array"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select the first element"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[0]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select the third element (remember, index starts from zero!)"
]
},
{
"cell_type": "code",
"execution_count": 104,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"2"
]
},
"execution_count": 104,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select the fifth element (remember, index starts from zero!)"
]
},
{
"cell_type": "code",
"execution_count": 106,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"4"
]
},
"execution_count": 106,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[4]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select the last element, without knowing the length of the array."
]
},
{
"cell_type": "code",
"execution_count": 107,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"4"
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[-1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Selecting multiple elements of an array "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* Select the first and second elements of the array.\n",
"* syntax: `[start:end]` exclusive of `end`\n",
" * default start: `0`\n",
" * default end: `N`"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0 1]\n",
"[0 1]\n"
]
}
],
"source": [
"print(a[:2])\n",
"print(a[0:2])"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* First to end"
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 2, 3, 4])"
]
},
"execution_count": 116,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[0:5]\n",
"a[:]\n",
"a"
]
},
{
"cell_type": "code",
"execution_count": 114,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 2, 3, 4])"
]
},
"execution_count": 114,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[0:]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Exercise: Select the second, third and fourth elements"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 2, 3])"
]
},
"execution_count": 121,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[1:4]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Exercise: Select everything."
]
},
{
"cell_type": "code",
"execution_count": 309,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 2, 3, 4])"
]
},
"execution_count": 309,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[:]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Advanced slicing"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select every second element.\n",
"* Syntax: `[start:end:stride]`\n",
" * default start: `0`\n",
" * default end: `N`\n",
" * default stride: `1`"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 2, 4])"
]
},
"execution_count": 122,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[0::2]"
]
},
{
"cell_type": "code",
"execution_count": 123,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 2, 4])"
]
},
"execution_count": 123,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a[::2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Let's make a longer array."
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0 1 2 3 4 5 6 7 8 9]\n"
]
}
],
"source": [
"b = np.arange(0, 10)\n",
"print(b)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* Exercise: Select every third entry starting with the second entry."
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 4, 7])"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b[1::3]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* *Exercise*: Select every third entry starting with the last entry. (Hint: Try playing with a `stride` of `-1`)."
]
},
{
"cell_type": "code",
"execution_count": 128,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([9, 6, 3, 0])"
]
},
"execution_count": 128,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b[::-3]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* *Exercise*: Select every third entry starting with the second to last entry."
]
},
{
"cell_type": "code",
"execution_count": 129,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([8, 5, 2])"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b[-2::-3]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## N-dimensional indexing\n",
"* Slices in each axis are seperated by commas.\n",
"* `[start0::stop0::step0, start1::stop1::step1]`\n",
" * Default: `[,]` = `[::,::]`"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0 1 2 3]\n",
" [ 4 5 6 7]\n",
" [ 8 9 10 11]\n",
" [12 13 14 15]\n",
" [16 17 18 19]\n",
" [20 21 22 23]\n",
" [24 25 26 27]\n",
" [28 29 30 31]]\n"
]
}
],
"source": [
"c = np.arange(0, 32).reshape(8, 4)\n",
"print(c)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Select the element on the third row in the fourth column."
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"11"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[2, 3] # [axis0, axis1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* *Exercise*: From every second row, select every element."
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0, 1, 2, 3],\n",
" [ 8, 9, 10, 11],\n",
" [16, 17, 18, 19],\n",
" [24, 25, 26, 27]])"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[::2, ::]\n",
"c[::2]"
]
},
{
"cell_type": "code",
"execution_count": 140,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0, 1, 2, 3],\n",
" [ 8, 9, 10, 11],\n",
" [16, 17, 18, 19],\n",
" [24, 25, 26, 27]])"
]
},
"execution_count": 140,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[::2, ]"
]
},
{
"cell_type": "code",
"execution_count": 141,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0, 1, 2, 3],\n",
" [ 8, 9, 10, 11],\n",
" [16, 17, 18, 19],\n",
" [24, 25, 26, 27]])"
]
},
"execution_count": 141,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[::2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Exercise: Select the elements in every third row, starting with the second row, and in every second column."
]
},
{
"cell_type": "code",
"execution_count": 142,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 4, 6],\n",
" [16, 18],\n",
" [28, 30]])"
]
},
"execution_count": 142,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[1::3, ::2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Slicing based on integers\n",
"* Select the elements at (1, 1), (2, 2) and (3, 3)."
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 5, 10, 15])"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[(1, 2, 3), (1, 2, 3)]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Exercise: Select the elements at (5, 3) and (1, 2)."
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 2, 3],\n",
" [10, 11],\n",
" [18, 19],\n",
" [26, 27]])"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c[[5, 1], [3, 2]] # Can use lists, ndarrays or tuples.\n",
"c[0::2, (50000, 100000)]\n",
"c[0::2, 50000:100000:50000]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Slicing based on booleans"
]
},
{
"cell_type": "code",
"execution_count": 158,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ True, True, True, True],\n",
" [ True, False, False, False],\n",
" [False, False, False, False],\n",
" [False, False, False, False],\n",
" [False, False, False, False],\n",
" [False, False, False, False],\n",
" [False, False, False, False],\n",
" [False, False, False, False]])"
]
},
"execution_count": 158,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c < 5"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {
"scrolled": true,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ True True True True]\n",
" [ True False False False]\n",
" [False False False False]\n",
" [False False False False]\n",
" [False False False False]\n",
" [False False False False]\n",
" [False False False False]\n",
" [False False False False]]\n",
"[[False False False False]\n",
" [False False False False]\n",
" [False False True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]]\n",
"[[ True True True True]\n",
" [ True False False False]\n",
" [False False True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]\n",
" [ True True True True]]\n"
]
},
{
"data": {
"text/plain": [
"array([ 0, 1, 2, 3, 4, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,\n",
" 22, 23, 24, 25, 26, 27, 28, 29, 30, 31])"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(c < 5)\n",
"print(c >= 10)\n",
"print(np.bitwise_or(c < 5, c >= 10))\n",
"selection = np.bitwise_or(c < 5, c >= 10)\n",
"c[selection]\n",
"# Select all of the elements of c that are less than 5 or greater than or equal to 10"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Mini-project: Mandelbrot set\n",
"\n",
"#### Key Skills:\n",
"\n",
"* Complex numbers\n",
"* Creating a grid of points in the complex plane\n",
"* Iteration\n",
"* Avoiding tight loops with operations on `numpy` arrays\n",
"* Boolean masking\n",
"* Basic plot"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The Mandelbrot set is the set of complex numbers $c$ that remain bounded under repeated application of the function:\n",
"$$z_0 = c \\quad z_{n + 1} = z_{n}^2 + c$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### A reminder of complex numbers\n",
"\n",
"$$z = x + j y$$ \n",
"\n",
"where $j = \\sqrt{-1}$ and $x = \\Re(z)$ and $y = \\Im(z)$ are the real and imaginary parts of $z$."
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Real part of z: 3.0\n",
"Complex part of z: 4.0\n",
"z + z = (6+8j)\n",
"z**2 = (-7+24j)\n"
]
}
],
"source": [
"z = 3.0 + 1j*4.0\n",
"print(\"Real part of z: {}\".format(np.real(z)))\n",
"print(\"Complex part of z: {}\".format(np.imag(z)))\n",
"print(\"z + z = {}\".format(z + z))\n",
"print(\"z**2 = {}\".format(z**2))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 1: Create a grid of points $c$ in the complex plane\n",
"\n",
"1. Tool: `np.linspace`. Create 300 points between -2 and 1 and place the result in `x`.\n",
"2. Tool: `np.linspace`. Create 300 points between -1.5 and 1.5 and place the result in `y`.\n",
"3. Tool: `x[:, np.newaxis]`, `y[np.newaxis, :]`. What does this do?"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(300, 300)"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = np.linspace(-2, 1, 300)\n",
"y = np.linspace(-1.5, 1.5, 300)\n",
"c = x[:, np.newaxis] + 1j*y[np.newaxis, :]\n",
"c.shape"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 2: Apply the function \n",
"\n",
"The Mandelbrot set is the set of complex numbers $c$ that remain bounded under repeated application of the function:\n",
"$$z_0 = c \\quad z_{n + 1} = z_{n}^2 + c$$\n",
"\n",
"1. Tool: `for i in range(0, 50):`\n",
"2. `numpy` element-wise operations `z**2`"
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"num_iterations = 100\n",
"z = c\n",
"for i in range(0, num_iterations):\n",
" z = z**2 + c"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 3: Select the $c$ \"that remain bounded\"\n",
"\n",
"1. Tool: Boolean operation `np.abs(z) < bound`"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"bound = 100.0\n",
"mandelbrot_set = np.abs(z) < bound"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 4: Plot the points $c$ that are in the Mandelbrot set\n",
"\n",
"1. Tool: `plt.imshow`."
]
},
{
"cell_type": "code",
"execution_count": 104,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASAAAAENCAYAAACvsbhoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFLpJREFUeJzt3WuwJGV9x/Hvz10W1IhcFmFh1wDllgYTE3SzQU2lUEwChGK9kawahag5RVKaFymr3BSVmKKSStQXUUtKWCMl5AIoXjgKBkRAkhcoS7lclnXjsqVyaldWXEJEFFz458X0mHHsmdNz6X768vtUnTpzeU7Pf/pM/+Z5np7uUURgZpbCM1IXYGbd5QAys2QcQGaWjAPIzJJxAJlZMg4gM0umNgEk6XJJ+yXdN+L+0yU9Kml79vM3VddoZvO1MnUBAz4JfBS4ckyb/4yIc6opx8zKVpseUETcDhxIXYeZVadOPaAiXi7pbmAv8J6I2JHXSNICsJBdfVlVxZl1VURomr9TnQ7FkHQi8MWI+NWc+w4Hno6IxySdDXw4ItYXWGZ9nqBZS00bQLUZgi0nIv43Ih7LLt8AHCJpdeKyzGwGjQkgScdJUnZ5I73af5C2KjObRW3mgCRdBZwOrJa0BLwPOAQgIi4F3gj8maSDwI+BzVGn8aOZTaxWc0Bl8ByQWflaPwdkZu3jADKzZBxAZpaMA8jMknEAmVkyDiAzS8YBZGbJOIDMLBkHkJkl4wAys2QcQGaWjAPIzJJxAJlZMg4gM0vGAWRmyTiAzCwZB5CZJeMAMrNkHEBmlowDyMyScQCZWTIOIDNLxgFkZsk4gMwsGQeQmSXjADKzZBxAZpZMbQJI0uWS9ku6b8T9kvQRSbsl3SPppVXXaO0QEalLsExtAgj4JHDmmPvPAtZnPwvAxyqoyVpIUuoSLFObAIqI24EDY5psAq6MnjuAIyStqaY6MytDbQKogBOABweuL2W3mVlDrUxdwATy+s25g3lJC/SGaWZWY00KoCVg3cD1tcDevIYRsRXYCiDJM45mNdWkIdgi8LZsb9hpwKMRsS91UWY2vdr0gCRdBZwOrJa0BLwPOAQgIi4FbgDOBnYDjwN/kqZSM5sXtf0zER6CmZUvIqb6bEOThmBm1jIOIDNLxgFkndT2qYemqM0ktFmVJP1cCPnwjDTcAzLDPaJUHEBmloyHYNZJHn7VgwPIOiVvqDV8mwOpOg4ga71J53ciwiFUEc8BmeXwpHQ1HEBmOdwDqoYDyCyHe0DVcABZq80SJA6h8jmArNVmHUo5hMrlvWBmI3geqHzuAZmN4N5P+dwDMsvh3k813AMys2TcA7LWG+zNLDescs+nWu4BWacUDZjhoPJ8UDncA7JG6gdCXqDM67M/40Jo+HH7x4/5OLLJuAdkjVaHnkr/MR0+k3MPyBovVeg4bGbnALKkBocuw0YNc+rC80Kz8xDMkhoXKKNOHlbXDX9UkNpoDiBLbrmNts6hM2hwLqhIO/MQzKw0g0PG5SbL6zS0rJJ7QJbcqI2vyUOa4bmtcQHT1fCBGgWQpDMl7ZK0W9KWnPsvkPR9Sduzn3emqNPKkzfp3GSDQ7KmP5ey1GIIJmkFcAnwu8AScKekxYi4f6jpNRHxrsoLtNJ1uRfQZXXpAW0EdkfEnoh4Erga2JS4JjMrWV0C6ATgwYHrS9ltw94g6R5J10paN2phkhYkbZO0bd6F2vz1hyhdGKb0e3qSfvbTZXUJoLz/wvCr8QvAiRHxEuBm4IpRC4uIrRGxISI2zLFGm0HTPtNTlq493+XUJYCWgMEezVpg72CDiPhBRDyRXf048LKKarM56Po7veWrSwDdCayXdJKkVcBmYHGwgaQ1A1fPBXZWWJ/NaNxR5l3UpWHnOLXYCxYRByW9C7gRWAFcHhE7JF0MbIuIReAvJJ0LHAQOABckK9gKGxU8Xd/w+rreM1TbXwiS2v0EG6Ltr7NpteUcQhEx1ROoyxDMWq7pG1hZuh7MDiCrjEMoX5fXSy3mgKzduv4uv5w2DMGm5R6QlarLG9ckuhrS7gFZabq6UU1r3Env28o9ICudg2hyXVlnDiCzGupKD8hDMLOa6EroDHIPyErTxQ1qFl0Zdg1yD8jmqosb0Tx1ba+he0AtNO4gx+HjsZY7IHLSQPE5bmbXpRB3ALXYcuFS5Du5HCbV6lqAO4Bapsi7Z17PZ/D6qJOHzbMGy9e1dec5oA6YNjym3Ria/HU6ddCleSD3gFom9QvXJ9maj66sQ/eArDDPC1WjS+vXPaCW6co7p7WDA6hFqgofh1y5urR+HUAtkOIF26WNpEreDW+NlCqEfKL5+eraOnQAmdVEl3o+fd4L1mBde7fsgq6F0MQ9IEnPlrSijGKs2bq28cxbF99Qlg0gSc+Q9GZJ10vaD3wT2Cdph6QPSlpffplWZ+MO4bDiuhjgy34xoaSvAjcD1wH3RcTT2e1HAa8C3gx8LiL+teRap9LWLyas28buwy9m1+QAmvaLCYsE0CER8dNZ26TSlgDKOz6oThu8A2h6TQ6evtK+GbUfLJI2SHr2uDY2X8NHrQ9er9vGXrd6rBkmmYS+Eniqf0XSaknnzKsQSWdK2iVpt6QtOfcfKuma7P6vSTpxXo9tlkobej+zmCSAfhIRP+lfiYiHgYvnUUS2V+0S4CzgFOBNkk4ZavYO4JGIeAHwT8D75/HYdeUeRft1PXxgsgDaI+msodtWzamOjcDuiNgTEU8CVwObhtpsAq7ILl8LnKGK/oNVDnnqOLyy+era4RbjTPJBxHcDX5L0VuAO4MXAA3Oq4wTgwYHrS8BvjWoTEQclPQocDTw8vDBJC8DCnGozs5IU7gFFxD7gZcBngGOAe+jtgp+HvLeD4W5AkTa9GyO2RsSGiNgwc2VU+47ld8f2cy/3/y3bA5KkyNZWRDxFL4A+M6rNlJaAdQPX1wJ7R7RZkrQSeC5wYIbHrDXv1m6/Lp16dZQiPaBbJb1b0vMHb5S0StKrJV0BnD9jHXcC6yWdJGkVsBlYHGqzOPA4bwRumTH0zJLr+ku4yAcRDwPeDrwFOBl4BDgMWAHcBFwSEdtnLkQ6G/hQttzLI+LvJV0MbIuIxayOfwFOpdfz2RwRewostxX/YX8Qsb3a0Asq85PQzwSeExH7JR0D/CnwI+Cywd3yddWWABpWt43dATS7JgdRmZ+E/jHwYUl30ZuH2QysBj4qac00D2qza/KL1ayvyCT0OuDXgF8H1kbEd7LbnwV8jNnnf6zh+mHoXtBsujgpXeRzQKcCe7I9YN8ZuP1w4PdLqcoayeEzm66FDxTbC3Y4cFfO7QvAVfMtxybhzwy1T9dCvEgA3QDszLn9YYY+D2Rm0+ta+ECBvWDQO1g0G4IN3ibguOwT0rXV1r1gg/pzB1W/gPN6X13ciOapqT3a0nbDN10XAqivyv/lqA2l7a+nqjQtiErbDW/NUeXxalaeLq1fB1DLdOnFa83n7wWzwhxu1egPY7uwvt0DapnUczD+aMB8dGUdugfUAYMv5gLH/uW+A08SbKlDsOm6Ej7gHlDrFHnx9nspg20Hr+ctY5KNoksb0Lx1bd05gFpsueHQuJ5K/+/cm6lW186W6M8B2Vy1/fVUhSb2gqb9HJDngGyu6nzStCZoYvjMwkMwK43DZzJdCx9wD8isNgYDuyth5B6QWQ11pffoALLSdeXdfJ66ss48BLPS+BStk+lK6AxyD8hK5RAqpovhA+4BWQWmPaSjK7oaPuAekFXI4ZOvy+vFAWSV6PJGNk6Xez/gIZiVzMGzvC6HkHtAVqrho+zHHXHfRV0P6OQBJOkoSV+W9K3s95Ej2j0laXv2s1h1nTab4VN/dF3eKVG6KPnR8JI+AByIiH+UtAU4MiLem9PusYj4pSmW3+23mJrI+9rh1K+9lNoWPI39Wh5Ju4DTI2KfpDXAbRHxwpx2DqCWSv0arFL/c1EOoJ7kQzDg2P6XG2a/nzei3WGStkm6Q9Jrxy1Q0kLWdtu8i7X569JwpB+2/ROPdSl881TSA5J0M3Bczl0XAVdExBEDbR+JiF+YB5J0fETslXQycAtwRkQ8UOCxu/0fbpgubpBtCN5an5AsIl4z6j5JD0laMzAE2z9iGXuz33sk3QacCiwbQNYMbRyW2PLqMARbBM7PLp8PXDfcQNKRkg7NLq8GXgncX1mFVonh3k/TAynv4wf28+owCX008Cng+cB3gfMi4oCkDcCFEfFOSa8ALgOepheaH4qITxRcfvf69A2z3MnxU79GpzE82Vzk65CarLF7wcrmAKq/tgbQKG0cbjZ5L5h13HIbY1OGMEU/5d2E51IVB5AlVeS7yYZvq+sG3MaeTdl8MKolNcmxYXXbuP3ljbNzAFnjpTjhWd3CsKk8BLNGGw6CFMEw2Itzb2gy3gtmnVNk3ml4PsfzO+N5L5hZAUXfcOvQs+oCzwFZ603Sy++3deBUwz0gM0vGAWSWw6fKqIYDyGwED8PK5wAyG8G9oPJ5EtpsBPeAyucekJkl4wAyy+HeTzUcQGY5PPdTDQeQWQ73gKrhSWhrvUmPlnf4VMcBZJ2Sdw4fB046HoJZJw33ijznk4YDyMyScQCZ4WFYKp4Dsk7yCcbqwT0g6ySHTz04gMwsGQeQmSXjADKzZJIHkKTzJO2Q9LSkDWPanSlpl6TdkrZUWaOZlSN5AAH3Aa8Hbh/VQNIK4BLgLOAU4E2STqmmPDMrS/Ld8BGxE5bdK7ER2B0Re7K2VwObgPtLL9DMSlOHHlARJwAPDlxfym4zswarpAck6WbguJy7LoqI64osIue2kQfvSFoAFgqWZ2aJVBJAEfGaGRexBKwbuL4W2Dvm8bYCW8FfzWxWZ00Zgt0JrJd0kqRVwGZgMXFN1lA+8r0+kgeQpNdJWgJeDlwv6cbs9uMl3QAQEQeBdwE3AjuBT0XEjlQ1W7P5MIz6UNvfDTwEMytfREyV6sl7QGbWXQ4gM0vGAWRmyTiAzCwZB5CZJeMAMrNkHEBmlowDyMyScQCZWTIOIDNLxgFkZsk4gMwsGQeQmSXjADKzZBxAZpaMA8jMknEAmVkyDiAzS8YBZGbJOIDMLBkHkJkl4wAys2QcQGaWjAPIzJJxAJlZMg4gM0vGAWRmyTiAzCyZ5AEk6TxJOyQ9LWnDmHbflnSvpO2StlVZo5mVY2XqAoD7gNcDlxVo+6qIeLjkesysIskDKCJ2AkhKXYqZVSx5AE0ggJskBXBZRGwd1VDSArCQXX2CXi+rLlYDderFuZ7l1a2mutXzwmn/sJIAknQzcFzOXRdFxHUFF/PKiNgr6XnAlyV9MyJuz2uYhdPW7LG3RcTIuaWquZ7x6lYP1K+mOtYz7d9WEkAR8Zo5LGNv9nu/pM8BG4HcADKzZki+F6wISc+W9Jz+ZeD3qNewysymkDyAJL1O0hLwcuB6STdmtx8v6Yas2bHAf0m6G/g6cH1E/EfBhxg5V5SI6xmvbvVA/WpqTT2KiHkWYmZWWPIekJl1lwPIzJJpXQBJ+qCkb0q6R9LnJB0xot2ZknZJ2i1pS4n11OpQkwnqqWr9HCXpy5K+lf0+ckS7p7J1s13SYgl1jH2+kg6VdE12/9cknTjvGias5wJJ3x9YJ+8suZ7LJe2XlLvzRz0fyeq9R9JLCy04Ilr1Q28P2crs8vuB9+e0WQE8AJwMrALuBk4pqZ5fofdBrduADWPafRtYXcH6WbaeitfPB4At2eUtef+v7L7HSlwnyz5f4M+BS7PLm4FrEtdzAfDRsl8vA4/3O8BLgftG3H828CVAwGnA14ost3U9oIi4KSIOZlfvANbmNNsI7I6IPRHxJHA1sKmkenZGxK4ylj2NgvVUtn6y5V6RXb4CeG1JjzNOkec7WOe1wBkq7/ihKtd/IdH70O+BMU02AVdGzx3AEZLWLLfc1gXQkLfTS+VhJwAPDlxfym5LqX+oyV3ZoSQpVbl+jo2IfQDZ7+eNaHeYpG2S7pA075Aq8nx/1iZ7g3sUOHrOdUxSD8AbsuHOtZLWlVRLUVO9Zpp0LNjPFDm0Q9JFwEHg3/IWkXPb1J9HqPpQkwrqqWz9TLCY52fr52TgFkn3RsQD09Y0pMjznes6WUaRx/oCcFVEPCHpQnq9s1eXVE8RU62fRgZQLHNoh6TzgXOAMyIboA5ZAgbfMdYCe8uqp+Ay5naoyRzqqWz9SHpI0pqI2Jd12fePWEZ//eyRdBtwKr15knko8nz7bZYkrQSey/ghSan1RMQPBq5+nN58Z0pTvWZaNwSTdCbwXuDciHh8RLM7gfWSTpK0it6k4tz3rBRVw0NNqlw/i8D52eXzgV/ooUk6UtKh2eXVwCuB++dYQ5HnO1jnG4FbRry5VVLP0PzKucDOkmopahF4W7Y37DTg0f7QeqyqZtErnK3fTW8suj376e+5OB64YWjW/r/pvYteVGI9r6P37vAE8BBw43A99PZ23J397EhdT8Xr52jgK8C3st9HZbdvAP45u/wK4N5s/dwLvKOEOn7h+QIX03sjAzgM+HT2+vo6cHLJr+Pl6vmH7LVyN3Ar8KKS67kK2Af8NHv9vAO4ELgwu1/AJVm99zJmj+/gjw/FMLNkWjcEM7PmcACZWTIOIDNLxgFkZsk4gMwsGQeQmSXjADKzZBxAVglJz5D0d9k5ml4k6SWStkq6Mju0AUnPlPRVSStGLGOVpNv77a35HEBWlT+k90nnT9P75O5pwF/TO/3GiVmbtwOfjYin8hYQvVNTfAX4o7KLtWo4gKwqx0bErfQC6PCI2BoRD9E7Edt3sjZvYeBYMPW+GeUzkr6h3lkuNwKfz9pZC/hQDKuEpHOAH0XErZIORsRKSX8M/DAirssOuvxuRByXtV8J3EXvOKgvSnoWvTMFPg58LyKOSfVcbH7cA7JKRMQXgZMlPZPeKYT/lt4BledlTVYD/zPwJ68FdmZ/R0Q8HhE/zIZnT/bPHmDN5sk8q0xEfAJAUtALn68Cl2d3/5jeEed9v0HvlLp5DgV+UlKZViH3gKwSkv5A0jpJLwAejogngF3AL0v6zYh4BFghqR9C3wNePPD3x2S/jwa+HxE/rfgpWAkcQFaVm4H3AP9O79sviN5ZDq8C+ifXugn47ezyJ4Fj1fsKoe30vrob4FVA/yu7reE8CW21IelU4C8j4q1j2nwW+Kuo0TeN2PTcA7LaiIhvALeO+yAi8HmHT3u4B2RmybgHZGbJOIDMLBkHkJkl4wAys2QcQGaWjAPIzJL5P2zb5mnJkRwPAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"plt.imshow(mandelbrot_set.T, extent=[-2, 1, -1.5, 1.5])\n",
"plt.gray()\n",
"plt.xlabel(\"$\\Re(c)$\")\n",
"plt.ylabel(\"$\\Im(c)$\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Mini-project: Markov Chain\n",
"\n",
"Key Skills:\n",
"* Calculating a transition matrix\n",
"* Defining a transition matrix in `numpy`\n",
"* Matrix-vector multiplication"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Wikipedia: A *Markov chain* is a *stochastic model* describing a sequence of possible events in which the *probability of each event* depends only on the state attained in the previous event."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Imagine you are standing in a five room house.\n",
"\n",
"*Drawing of floorplan of house*\n",
"\n",
"The doors between the rooms of the house are locked.\n",
"\n",
"Problem: Every hour the doors of the house are unlocked, and *you must* move to another room in the house. You have no preference over which room you want to be in next, so you move randomly between the rooms. **Which room will you spend the most time in?**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Now imagine you are in room 1. You could move to rooms 2, 4 and 5 (three rooms). You cannot move to room 3.\n",
"\n",
"Denote $m_{ij}$ the *probability* of moving between room $j$ and room $i$.\n",
"\n",
"What are $m_{21}$, $m_{41}$ and $m_{51}$? What are $m_{11}$ and $m_{31}$?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$m_{21} = m_{41} = m_{51} = 1/3$$\n",
"$$m_{11} = m_{31} = 0$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"${m}_{1} = \\begin{bmatrix} m_{11} \\\\ m_{21} \\\\ m_{31} \\\\ m_{41} \\\\m_{51} \\end{bmatrix} = \\begin{bmatrix} 0 \\\\ 1/3 \\\\ 0 \\\\ 1/3 \\\\ 1/3 \\end{bmatrix} $"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Exercise*: Write down $m_1$, $m_2$, $m_3$, $m_4$, $m_5$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The transition matrix $M$ has entries $m_{ij}$.\n",
"\n",
"*Exercise*: Write the transition matrix $M$ on paper.\n",
"\n",
"*Exercise*: Write the transition matrix $M$ in a `numpy` array."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"$$M = \\begin{bmatrix} \n",
"0 & 1/3 & 0 & 1/3 & 1/4 \\\\\n",
"1/3 & 0 & 1/3 & 0 & 1/4 \\\\\n",
"0 & 1/3 & 0 & 1/3 & 1/4 \\\\\n",
"1/3 & 0 & 1/3 & 0 & 1/4 \\\\\n",
"1/3 & 1/3 & 1/3 & 1/3 & 0 \\\\\n",
"\\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 232,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0. 0.33333333 0. 0.33333333 0.25 ]\n",
" [0.33333333 0. 0.33333333 0. 0.25 ]\n",
" [0. 0.33333333 0. 0.33333333 0.25 ]\n",
" [0.33333333 0. 0.33333333 0. 0.25 ]\n",
" [0.33333333 0.33333333 0.33333333 0.33333333 0. ]]\n"
]
}
],
"source": [
"M = np.array([[0, 1/3, 0, 1/3, 1/4], \n",
" [1/3, 0, 1/3, 0, 1/4], \n",
" [0, 1/3, 0, 1/3, 1/4],\n",
" [1/3, 0, 1/3, 0, 1/4],\n",
" [1/3, 1/3, 1/3, 1/3, 0]])\n",
"print(M)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Question: Which room will you spend the most time in?\n",
"\n",
"Let's assume we are certain we start in room 1.\n",
"\n",
"$$ p_0 = \\begin{bmatrix} 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\end{bmatrix} $$"
]
},
{
"cell_type": "code",
"execution_count": 310,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"p_0 = np.array([1.0, 0.0, 0.0, 0.0, 0.0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Where will we end up next?\n",
"\n",
"$$ p_1 = M p_0$$ "
]
},
{
"cell_type": "code",
"execution_count": 257,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0. 0.33333333 0. 0.33333333 0.33333333]\n"
]
}
],
"source": [
"p_1 = M@p_0\n",
"print(p_1)"
]
},
{
"cell_type": "code",
"execution_count": 259,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.30555556 0.08333333 0.30555556 0.08333333 0.22222222]\n"
]
}
],
"source": [
"p_2 = M@p_1\n",
"print(p_2)"
]
},
{
"cell_type": "code",
"execution_count": 263,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.11111111 0.25925926 0.11111111 0.25925926 0.25925926]\n"
]
}
],
"source": [
"p_3 = M@p_2\n",
"print(p_3)"
]
},
{
"cell_type": "code",
"execution_count": 262,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.1875 0.1875 0.1875 0.1875 0.25 ]\n"
]
}
],
"source": [
"p = np.array([1.0, 0.0, 0.0, 0.0, 0.0])\n",
"for i in range(0, 50):\n",
" p = M@p\n",
"print(p)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"*Question*: What happens if we are not sure where we start?"
]
},
{
"cell_type": "code",
"execution_count": 313,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.1875 0.1875 0.1875 0.1875 0.25 ]\n"
]
}
],
"source": [
"p = np.array([0.25, 0.25, 0.0, 0.25, 0.25])\n",
"for i in range(0, 50):\n",
" p = M@p\n",
"print(p)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Mini-project: Calculating $\\pi$ with Monte Carlo\n",
"\n",
"Key skills:\n",
"* Simulating random numbers\n",
"* Boolean operations\n",
"* Slicing with strides\n",
"* Counting\n",
"* Algorithms with no tight-loops"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Imagine you are blind-folded and you have an infinite number of darts to throw at the following square board:\n",
"\n",
"*Draw square with quarter circle*\n",
"\n",
"How can we calculate the value of $\\pi$?\n",
"\n",
"*Monte Carlo* (random simulation) provides a powerful answer to this problem."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Steps:\n",
" \n",
"1. Randomly throw (many) darts into the square (random simulation).\n",
"2. Calculate the ratio $r$.\n",
"3. $\\pi = 4r$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 1: Simulate from uniform distribution\n",
"\n",
"We want to 'simulate' the idea of throwing darts randomly at a square board.\n",
"\n",
"Computers can *generate* (pseudo-)random numbers.\n",
"\n",
"1. Tool: `np.random.uniform(size=())` to simulate N draws of uniformly distributed random numbers on the board.\n",
"\n",
"Uniformly distributed - equally likely to be at any value between $[0, 1]$."
]
},
{
"cell_type": "code",
"execution_count": 106,
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f92e178b240>]"
]
},
"execution_count": 106,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXtYFGeaN/yrbk6CnEGgAUG6RRREGjwFRXzVqPEwI37vZnaM7mgSo8647sQ5ZDWbd7PONWZ291vne80BlSSaZeJOZhJhgiRqogFRjAp2bEERuxEEukFOAoKcuuv7o3keq6qrDyDEKPW7Li/p7qqnnuepqvt+nvvwuxmWZSFBggQJEiQAgOxxd0CCBAkSJPxwICkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBAISkFCRIkSJBA4fK4OzBUBAUFsdHR0Y+7GxIkSJDwRKG0tLSZZdlgR8c9cUohOjoaJSUlj7sbEiRIkPBEgWGYGmeOk8xHEiRIkCCBQlIKEiRIkCCBQlIKEiRIkCCBQlIKEiRIkCCBQlIKEiRIkCCBYtSUAsMwHzIMc5dhmDIbvzMMw+xnGEbHMIyWYZjk0eqLBAlPOg4U6lGsb+Z9V6xvxoFC/YieI0HCaO4UjgBYbuf35wBMHvz3CoDMUeyLBCfwQxMiP7T+OIvR6HdihC+2H9XQdov1zdh+VIPECN8RPWck8aTev7GOUVMKLMueBdBq55AfA/hv1oJvAfgxDBM2Wv2R4BiPW4j80PvjLMT6/dKREsgFb9tQBGSqMgjvrFNj+1EN9p26ie1HNXhnnRqpyqARPWe44CoA8rdcBrx0pATF+mYU65ux65h2RO+fUOkcKNQjq0jPm1NJCQ0dj9OnEA6glvO5bvA7KzAM8wrDMCUMw5Q0NTV9L50bK+C+WESIbMkuxQtZ3z6yEHnUleL3KdRGAmS8h85W4bmEEGw/qsGrn3yHl46UIFUZgP84cfORFFyqMgjr50zE/jM6rJ8z0al5GM45tsbF/XvXMS12HdPSsXxZZsRLRy5j1zEtEiN8sSW7FPtOVeIZZQC2ZJfixSOXcVxrdPr+OfPsCJWvXAbsza+gyvdJWUT80PA4lQIj8h0rdiDLsodYlp3JsuzM4GCHWdpjAvZemqEIY+GLBQD9JjPO61uGLURstf19CUJ7GE3bPBlvuJ8Hjl6shcLXAzmaeiiDvXCmognPz4x4JAVXrG/Gny7ewY5FKvzp4h2rPo3UOUJw7yMR+LmaehzXGvHikUvYkl2KBIUPevrNyNXU468ltXjQZ0JPvxm+49zQbzKjp9+MTanRTo/XmWdHuGjILKjC7pVxyCyoeiIWET9UPE6lUAcgkvM5AoDhMfXliYO9l8bRC2Vvd7AluxSuctkjCRGCkVjpj4RQ42I0bfNkvF+WNSJe4YMyQwcCvdxQZujAujmRiAzwQnpsME/BObtzItd8Z50aO5dOofNqbz6Gc44YuPfxW30LAMBFLsOSqSE4U9GE3gEzPr9qxO6VcRgws8jRGMAwwBp1OHI09QAw5Pvn7LMjXDRsTlOO6CJiLOJxch99DmA7wzB/BjAHQDvLssbH2J8nCtyXZv2cifjTxTu8l+addWq8dKQEyxNCUFjZTH8r1jejpqULh85W8Y4nuwMPVxk+3DgLqcogzFUGPvJqi/vS7likGrJC4F5/pPpjb94Ai9JMjPCFtq4diRG+SFUGYdvCGDqfJ8oa8cHGmaJ90Na1Y1qYN87pWhDo5YqWrj5E+o/DsSsGrE1WIFdTjwx1OP508Q68x7kgs8ByHxxBW9eOd9apoa1r542DfNbWtWNrulL0HNJP7jlDnT/hfQSA/Wd0mB3tj0vVbZAxwHVDJ/pNDzf7J8sb4OEqg6tchrnKwCHfP2eeHeGiwXucC+/zXGWgpBiGiNEMSf0fABcATGEYpo5hmJcYhtnKMMzWwUO+AFAFQAcgC8DPR6svTyvsmVZSlUFYnhCCHI0B6bFBVCFsP6rB6hkK3iqM7A5SlYFw5XhDhYJnOCAvbaoyEIeLq3krRXur5AOFeuRdNVABcqBQj3JDO5bFh9D+DNeJaGveyA6K7AzkMgzawy8hs6AKKVF+yNEYsDwhxKagkcuAc7oWRAd6oqWrH5H+41Db9gBxoeNx7IoBu1fGobCyCemxQdibX4FtC2Po+OyZqLamK5GqDOLtWrifxXY65BzhznBrunLIc8cVvoeLq3G4uBoZ6nBcrm5DhloBs5lFjqYecgaYrwrCgIlF34AZv1oaix2LVXjpiIXEcijPk6NdonAntG1hDJ3TR9kZjXWMZvTRT1mWDWNZ1pVl2QiWZT9gWfYAy7IHBn9nWZb9BcuySpZlp7MsK1GfDhH2XppifTMKK5uRoVYgV2PAq598x1ulcQVjv8mMgxtScHTzXBzckMJ7kYgQGW7/yDW3D64ut2SX0mgUe2abxAhfnCxvpJ+JEzEm2IsKteE6EW3NGxGwgEV47T+tQ2+/Cd9UNGFamA/O61qQoVagsLJZVNAU65uRWVCFxXHBqG7pRoLCB3VtD7A4Lhia2nasTVZQ80aOxoA1agVMZvCu7ayJaigmOWfbtqWYSNTQO+vUmKsMBAAMmMz4+kYjdq+Mw9c37sI8uEGQMQwuVbdi3ZxIjHOT4wutEftO3cLOpZPpDsWZ58kZ05dwJ2QyA7tXxtE5HYlFzViE/M0333zcfRgSDh069OYrr7zyuLvx2MF9aZ6fNRHTB1/86RG+qG3rpr9tTlOipqUbOZp6LE8Ixea0GHr+H07chHqiH5rv92H1DAUiAzwRGeCJ6YOmk5nRAY/Ux/xrRmyaZ3EuRgZ4IinSD8e1Rtxp6cZ/X6ixK8xIP7Yf1aCrdwAfnq/GzqWTkVlQha7eAfzhxE275x8o1KPPZEZkgCdvzt4r0OHtMzrReUtVBtHPAZ5uuFzdit4BFkHj3VDR0IkMdTj++BM1qlu68B8nKlB/7wGWTA2hbR86W4U3Vk1DSU0bIvw9UVLThgx1OBIj/ZA+JQjlhk5M8HHHH07cxOa0ScjTWuaHO+9kvLbGd6BQjxAfDwR4umH/GR02p01CRIAn8q8ZefeLO37S9pbsUpy71Wxz7vtMZjoXkQGe9BnrM5nxwpyJWJMUgfxrRrycNglNnb3oeDAA9UQ/XLzdiimh3lgydQLKDR1gWaCquQszInxxqboNrz03BZvTlE49T6Tf2rp2+uwU65vpZ+5zOTM6gHd/Z0YHICUqgHedyADPR36Onxb827/9m/HNN9885Og4iebiCYU9ezH3N8uOoQkZ6nCcKGvArmNaZBXpqdI4unku3d4PZ3dgz+xBzBcEqcogbEqNdjqyKVUZhLhQb1EnYlyo97BWxwBszhv5TK4BABnqcDTf74OLjMHXNxpRrG9GTLAXevrNaOzo4bX9ygKLKeiVBTG4buzAjkUqFFY2ITHCF5vTlHhlQYzd1a8zkVYk+udwcTU15WzJLrVa9duKKiu2M/e2diGvLbdE9BTrm+lz8V1tO16cH419p25hx2IV/vaL+ViZqAAjYyBjAN9xrjg3uLOKV/g6barimsO4Js/ECF8cOltlleuRVaTHxsOXnGpbgnOQdgpPKISrJODhqoj8xt1NbE6LQXKUH/Z9VYmiW834pyUqrEmKQLG+GXvybmDHYhUa2nuHvKrqM5nx0pESuLsySIkKoNecPckfl6vbUFLTRlesZHeyPCEEfympQ3KUn9UYuCjWN+ODc7chYwBtfTtc5QwOnq2CjAGa7/chKdL2+bZW3pvmTbI5b+Sav8u/ARljiZmubunGb5bF4sqdezCZWeR+V4+S6jb8elksjmsbrFb19nZw3NUvt49k9UvmZ3PaJPzp4h26Yueitq0bx7WWeAx1pB++q70HAHSnZ2v8v8u/AbmMwdYFMTbbJud19Q7QXcjzsybanMs1SRFIjvLDnrwb9Pt/WqJCsb4Vbd391Al97Eo95k8OREpUAN0J1LZ1091Nsb6Z/m1vx9R0vwd78yvg5S5HSlQAsor02JtfgRfmTkRKlLQbcARndwqSUniCYctEQl4wrvkGADXhyGUM/ny5zuoFnxkd4LBNISIDPOHuymBvfgVqWrqRWaDHtoUxyCyowqZ50Qjx8cD2oxrIZMCevBvYtjAGn5bWY8diFfbk3bApnIhwzVyfjB8nhSNXU48zFU2QM8CHG2dh9QwFz9QhBjEBZw8k83hOTAD+7ccJkDEMKhvvY05MAIK93THOVY7bzd2YGeWPvRmJom3nXzMiaLwbpoR680w3eVcN8PN0Q4Y6wqqPRDDaUibc8RHzTfB4d+w/o8PWBTGI8B+HysZOasoiY9HWtWPq4E5LxgDv/2ymaNvce25LcduaS+H37q4yFNxsgpuLDHc7ewFYlOvl6jYkRfohxMcDW7JLcVxrGQcxdRIzmr37lhIVAC93OfbmV6BY14w/X67DC3Mi8ZtlU516Vsc6JPPRGIAjB6LQfANYTARRgeLx8hsPX4JcBl6bWUV6vHSkxK5Dd3Oaksakq4K9aJglcWi/s06NfaduIT02mP62OU1p1wnINYGlKoOQPNEfAJA80Z/Xrj0n4lBzHLR17di5dDK+q7W0+dbaROxYrMK+U7cQE+yFioZO7FikQkVDJ7KK9KJtb01XUoXFvd7J8kbROSTmN6HJTxiCSkDMN9xrxwR74WR5I73ermNabMkuhVxmOW6eMhAMgLyrljQg4dyR54iYFbctjEFhZTN2Lp1Mx2FrLoVRSf95ohKvr4zD1gUx6Ok3w8yy2Lk0FqsSw7D9qAbvntFhYNAT/K2+hV6PO057921zmhKzBncgU0LH48uyRpvPv4ThQdopPMFw1jkphKa2DVlnbyNDHY48rZGu4tcmhyOzoArbFsZgT94NVDR0IOvsbfz2uSlYkxRhsz1L1I0e08N9cKm6DcsTQrA57aFPIjLAE30DJnx0ocZqlWlrRcc1jxXrLc7RzWmTcOZmE13h2jvf2ZW38JopUQFWDu4di1VUmbV292OCtxve+6bKIjQXTYZMBuz6rIy3qnb2vhDnLtepyv0sHJ/YuMgOjJhx/lJShwGTGZer25C5PhmzJwXgi2sNqGy8T01ukQGe1LRHHOy7PivDjAg/HLtSj4MbUrAmKYLucsSc8+S5Id/XtnajqrkLc2IC8OH5amxOmwRtfTuCvd3x1lrLzurTK/VwkTNYFh+Cjy7UYHlCCD4trac7BUf3LatIj08u12F2tD+u1XVgbbIC+0/rhvT8j1VI5qMxguGYSPbk3cDOpZPxaWk90mODkHX2NnYunYzNaUoqZFTBXjhd0YQMtQK7V0yz2x5Z7Z26fhfLE0KQqzFQuy85htjKMwv1aOx4gMWDpo4DhXrkaOqQ/W0N1qjD6fFce/NQhTsAUdOZs1FVwjl1kctoW30mM/79y5t4flYEvD3c0Gcyi/pk8q8ZYbj3AMe1RnpfSPST7u59aqoL8fHAiulhPNv/iumh2DRv0pDG1dDeS01FWxfEYGKgJyob7yN4vCXaKXN9MlbPUPDGz402SlUGQd/UiVPXG7F6RhhV6pEBntDdvS96zVyNAW+smka/XzI1BK5yBvtO3cLBDSl4ftZEJEX64e0zOshkoIriyp02XKvvwJTQ8Tiva8HOpZOpf4tEcIndN01tG/bmV2D3yjj84f+ZAS93Od77pgopUf74a2mdU8//WIazSuFxZjRLGAEIt9qOMji5ZorOBxbBl8GJl09VBiE9Nhg5mnrMjvanMfm22tTWtVMfAml3msIH+07dQrzCsso8rjXi4IYUpCqD0NLVi48v1qLc0IHcX8xHbWsXPr5YixfmRNLxECUg7C/pnzNZuWLRU8TsZAskkxkAzySyKjGMtkdMTJkFVVg/ZyJViCYz/5pyGVBu6ICbiwyHi6vhPc4F+09bIpoObkgB8NBs8846NY068nCVYfUMxZDHBVjMfuQ5eGedmvoduNnA3ES9xAhfGm1kuecGxCu8re75UObSZAYv25tkg+87dQsfbJwJADhcXA0XE4ubDfcxXxWIzIIqAOA9Q2LX2nj4EnavjKMKa3OaEtXNXfjkcp2UwTyCYFhWlIPuB4uZM2eyJSVSnhtgTQMh/GwPu45pcVxrxKbUaCpEAOD9oip8U9GENepwFFY2WQl8MRABwxU4cplFQNS0dOG41ogdi1UwmS2CcNPhy+gdMGN2tCVCad2cSHxZ1miTduL7QrG+GVuySwE8FNzcz9w5To8NQo7GQJPZuH3m7p72n9ahu88Ek5mFp5sc7/9sJs9vQEJM+we1sqtchlWJYVg9Q8GbA3K8mIAWew5Iv7n3V9geOeevJXXI0dTDzUWGI5tmAcCIkslxnw/ybOw/rcPEAE9cN3RgnioIpTVtNqlDxNohY9iSXYpViWF4a23ikJ7/sQiGYUpZlp3p6DjJ0fwEg7uKJnHgXAcil96YgGSokrDGucpASoj38kclOHerBbtXxuGPP0nCO+vU1Mdgz6ErdGgnRvgis6AKiRG+1Fm7N78CNxs6sf2oBoc3zaLhirOi/fH7jMQfBIlZqjIIqxItJT2IE/TghhQc3JDCy2PYtjAGuRoDZkf7I1djoFQVBOS+bE5TYlNqNEyD6b6TgrysYu8BUBbR5In+2LFYheNaIy/zm2QU17R0ieaECHmsCFYlhtnMBiY7ri3Zpci7aoCbiwzuLjLebyOVCcx9PsizcXBDCvJ3pGGNWoFzuma71CEEwsAK4jgnOyspg3lkIO0UngKQ1Rd3RZ9VpMd/nqyEu4vMapW7LD6EvkiEGO5wcTVigr3w2vI4p1eowpUb93jyApPVPzFJ7Vikgvc4F+zNr8CswZ1CUqQvdE1dmB7ui2v17XSVTtqxdf3Rwr5TN6nZZefSKbzfHu4ULOPJGNxRiQllspLtN5lhNrPoN7FYI9hZcHdsh4urAQA7FqtQ1WTZYfWbzDCzwG+WxSJe4cszreVdNeBkeaNNQj9H97FY34wXj1xGT78ZOxapeIR1XDJAe20MFdy+cXdd9kgGubC1w/2+n5EnEc7uFCSl8BSAa67ILKhCemwQcjUW8jUiSGyZZogAtFA4JPHadPSicbfr2rp2K8X06icanoll/ZyJeP/cbXT3mfD6oG349RwtPr5YC3cXGX69LBb7T+vQ02+Ci4zBzqWxPBbR7+PFJ2Oyx6BKxkmOEfMpiJmifvbhJfSbWGSoFfjjT9SiSpqcs2TqBORdNWLAzPJMVOSYfpMZrvKHCn84sCdguf6OoZomncFwTZ/2FJlkMrIPyXw0hpCqDMKy+BDsP62DKtiLkq3FKyyrbFumGeKkzlBbksOyivT0e3vmCmKq4tIi3Gzo4LF+ZhXpkasxYL4qEDkaA55LCMHOpVMwOWQ8PN3kACwCNjLACy/MiURcmDcyC6qwZGoI+k0segfM2H9axxOENS1dTs/JcIrpOEPCVtPSRftFjtl/WmfVN21dO1YlhvGEtoerHPEKH5wYjK0nZiaS1wBYFMjEgHG0JoGFKqMZ2xZaKDK+1beIFq0ZatElcn8ObkjhjRV4aO4ZKvHeUKCtszDeEpDr5V01OKTEcJXL4OEqw6GiKmzJLpUUwghDij56SrB6hgKfldbhUnUbZkf74+sbd/H1jbvYsVglGp0kXJlNU3hjb34Frhs6eKtSsZ3AsvgQKmS4Side4YPMgipcN3QiV1NPGSujAj1x9GItooO88LdfzEdWkZ5Go3Bf5heyvkWOph5xIeNR0XgfPf0m/LWkDvlaI9wdROUIYW+lawsjWX9AuGsg/glbq2JuFNB1Qydc5Qw8XOW8OgQk4dDDVWZ1Px2NV/gbYS61N9ZUpe16Bs6aqOzNj9g8EHOYGLjz+K2+hfJTjUR/JDyElKfwhINQFJQb2lF4sxlucgYNnb0YGLBEs5AEptbufsye5E+pJbR17Zg9yZ/G1qdEBaCmpQs5GoMo5814dznNZ5g/OZhSFTyjDKB8RPd7B7AobgK1te9eMQ0zowOweGoIvNzl2HfqFvoGTPjwfLWV2YPwHAFAY2cv5qsCUd3SjYqGTsjlDD74mWN7Mxf2EshsUXk0dvTapKEgWDI1BLWt3dj3VSV6+kw0B2BKqDfe+FsZJvi4W7VrL/aetB0Z4ImKhk7kaOoRr/DBvp8k0R3EiulhiAr0xEfFNXCVM3CVy7BtoZLmNwgZXoXjFZsLkpwmNlYuPxHJLzlcXI3a1m5KpWGLUZVLV/Eo90gMJEcDAO2Xtt7iVF4QG/zI/XnaIdFcjBGQkMb/92Qldq+Mw8pEBfoGzGBkDJ5RBmJVYhhdSXIjibgRQsDD+gtCagGyWiQmqcyCKlqSccBkxn+dqgRg4SPasVhFK4sVVjbxTBYmM7A8IVS0FCV5gXcsVsHdRQZPNzku3m59OEgHfi9bphNbprNHrR29eoYC/SYzbRewrMTnqQJF2yXsqVykKoOsdhMnyhqQoQ6Hsb2HHkNMKpkFVVibHI4PN86iNS8AfrQZd2VPxkvmhvtbemyw3Qgd8kxtOnwZ2xbG0BoKx7VGZBVZzFEjZV4S67MtkPnimvgObkihdTdG09w1liDtFJ5wRAZ4ora1G/qmLnh7uCBXY0CGOhx3WruRPNEfb61NpMeRbOWpod7418+vY8fih0ypxFHtIpchaLwb9n1ViaRIP7paJGUoJwZ44rjWiK0LLHUZbjd3Y+X0UFQ1deHPl+uwc+lkBHi5Y9O8aN7KzRa1BuHI3zQvGg3tvXg5bRJ6B0woN3RCLmOwJikcurv38WVZg01WVFur1tmT/GkWLZcZdLj0IAS1bd3IuWKAXM7gyp02fHGtAZnrkzHB2wMmM4u3zwyNdoFrFtmcFmNV44FkFPcOmBHi40F3Bdq6doT4eKCxo9cmyyqXkPDD89VYnhCKXE090qcE2WQWJc/UzYZOFOtb8MW1BhzckIKYYC/sO3ULryyIsUuSJ4Q9ksU+k9khMywX9jLVM9QRQ8ruH2uQaC7GEJZMDcGtxk4a6fPHn6gpvQD3JeO+xKtnhOHT0nqeKYkwm04J9cZxrRGNHT348Hw1ZTb90YwwnK6wMGBq7rSh6X4fVk4PRa7GgGBvd/zLqmmUbZX7shIqiJ1LJyP7wh3MiPDD0Yt3KL0BKaqyNV2J2rZu/OeJSsjlDMa5yvHac3FYow6neRVcJlACMSHPTboTo8cYKj0IARHg/7REhW+rWtDTb4aMAVQTxmNP3g28+mwsrwCOM+1yBR2X+oKYl8j8JAqUBZc3iVtYiTteYnram19B63WTjGx7AnjJ1BAU61twu7kbMgYI8fagPFDE5CimhIaqtLncSdx7RBTGUEx8zvZnrEJSCmMIxHG4ekYYCiubeTZmrt2a+9LkaY2URG1qqDc+PF/Ns0EnRfrh3W/0WDLVQli2bWEM/ueSZScgYwDd3S6wLIs7rQ9obYEV08OsOP25FN5rkiKgb7qPU9cbMV9luQ6pxzBPFYjeAYvwWD0jDKsSwxDs7Y63z+iwYrolw7d3wGyTt8geXxH5fai1C8RAxjLB2wNfXGuAjAHMLFCsa8Ghf7CEng61XS75HxGgK6aHIUMdwbON2/MbvPG3MqxNDqd+gsgAT8hkQK7GgHA/TwR7u1F/0fZFkx3yQBEfDx2fvoVHXmdLCYmN19bOLFdjsNlnrr/EGR/BcDmyxhIkpfCUg+sMJKaHeIUvNV8QxcBVCGLsmumxwVbspYA1sykRsmuSImBs70GwtzvKDR20tsD0wcpYYo5WsrojbKrLE0Lx1fVGBHu7IbOgirKQmswsXn02FlNCvemqm6yYM9QRdonshEKeqxC4Y3oUkj3AIsDJnGeuT0aItweK9S1wdWHoboGwqXId+/nXjNDUtuFydRvvngi5/x05yLmlOJcnhCJe4YP8a0asmB6GN3LLqTOYEB++sWoaQnw8kFlQxVNU3GdDbC7J+FQTxqPgZhNcZAyuGzvp7s5ZwkHynKYqg6jSXp4QApMZWDE9jFdTg9tnewpQzBz1XoEOqgnjKZHgUAgQxwokpfCUg6woiSAFLA44riDlvgxiL7FMBmRfuIOfL1RarWqL9c3418+vY/WMMFpLmDiHL91uxemKuzwq61RlECb4uNtc3fHrRsfQSKflCSHYvWIapkf44u0zOgR4ullFzjh6qYv1luI4Oxar6CqY2NG5QtjeXAxFgJBCOgNmM48i+nZzF9Jjg7Fp3iRqMtu2MAYN7b2Qy4C9+RXUlm9v5WvLtNVnMmNLdim+q72HldPDkKupxxfXjHhl0L9zXGtEZeN9GhUlDEclCpBbKY87h8LiTIDF77M8IRTXjR2IV/jgm8H7naGOsFvBjoA8pw99GhYW3fQpQZSa25Zvx948CJ+zt8/o8OqzsaI7VQkWSErhKQcRZM4KUmH5TrIqIxTH3NUyEeA7Fquo6WhP3g3IZMAbueWobLyPzPXJdgvfC19yriDedUyLE+WNWDk9FCfKGpEc5QcAuN3UhU+v1FMB4GwVuPxrRswbZNsk/ZDJwHOK2psLMp/OFpZPjPDFgNmMvfkV2Ll0MpKj/NHY0QPNnXv4UZICKVEPfSpc8xyx5TtyQtsybXFLcVoiiO7RQjaZBVXIXJ9MmVHVE/0we1IArwQoiciapwrEvlO3aP0HoYIi85N/zYjZk/zxaanlnhTrW6iSc1bYksXH3vwKxAR7QVN7D79eFkvvFQBUCe67o3l41ECBsQpJKYwBDNdZCthfLTd29FJTEdfMlH3hDlQh47HhmShRO/AadTivT0QwcQVNVpEe7xdVw0XO4LXn4vDjJAW2ZJfib98ZUNf2AD9KUiBPa6SRM1uyS1Hb2g3d3fvQ1LbRiCUizPKvGbE1XWlVHOfdb/TU1EEgplCcBVFQJJon0MsNiRG+eL+oGp9fNaChoxe/XhYLkxm83APu/dm+aLLD+2XPtKWta+eV4vz5QiUUfh7UVxAR4EmF6MnrjbTkJTeZjdzX5Cg/h0KV7HaEJkfu7sYZxX25ug3B3u44e6sZZjMLuYxBRnI4Lla1Yt9XlTC29yA2xBunb9y1UlQjHSgwliEphTEAspJST/TDyeuNvJDbtyu+AAAgAElEQVRNRwLQ3mqZ+1v+NSPGu7sgR1OPny9U4u9nT7RpuyYvs1AwkX4V65ux67My/HpZLF6cPwnbj2oQ4OmG0po2DJhY/GZ5LG9nEhXoiYu3W1HZeB+NHT04dqWeV5BlS3YpQnzcaUSSWHRVdUsXPN3kvFrA3KLxYhATdJraNuz6rAw/TlJgxfQw7PuqEmWGDrAsi94BFlsXxGD7osm0mhl3Lki9Y3dXRjRElgt7yppEZ3ETyqqau7B1QQwOF1fjuNZId3BJkX44rjXiuNbIMyc5Ms042xexQj2RAZ7YdUyLfV9VIsTHHZ5uchqOnH3hDlbPCENl433om7pwsaoFlY2dYAGwLIvO3gH88tnJvOTK0QgUGMuQlMJTDu5KavakACoAkiL98F6BDvu+qsTLaZPsKglnVnnC/AKuoBYKG+HqTkwwkUxarlBS+I3DtoUxvMpv08K88cnlOny4cRaCx7vjuNYIVzmD0pp76Okz4Xf5NwAAsSHeVPhwhfCJskbsWKxC9rd3kKOpp7kEABxmuorZrEl1tT15NxDg6YbL1a3o6TfD1UXG88lw8wKIT4GE83JLeNpybpfUtCHEx4P3XW1bNxo7emm/3lmnRsRgvggAbFtoSeqqbLyP1TMUdCWdFOmHOy3dQzLNcOGMmU1oyvlLSR0AS6LinrwbuNPahfe+qcLzsyLw9Y0m/PLZybhY1YI+EwvTYFKim4scBzekoKG9l2bZb01X8naDpC+PGigwliFlND/lEBa2J0yc757RUWFBYCtj11FmryVaqAq7V8ahsLIJ6bFB2JtfgR2LVdiUGm2VhSrGHXRwQwoSw32tjiVkfDsWqdDa1Yf9p3U08zY9NgjndC1YPcNS24Ac5+EqR2+/CfvP6NBvMuPghhRKBSFWdJ4Q7PX0m9FvMtMaCY7sz7aydTenKWn2bb/JDA9XGVzlMsxVBmLbwhi8dMTC3vvOOjX2nbqF9NhgmisRGeBFuaC41xBmFtu7J9z51da182o9vLU2kVf3geBGQ6dVlrozxH9DATcreVNqNA5uSKFsvR9frMU8VSC+LLNwGsUrfOHuKkeE/ziwLNBvYimxn1iWvdhza4+jSsKjQ6LOforArQNASNQcVTMTUkUTGmcS+kdeyHfP6HBe34IMtQIernKcLHeuUpoYFTVgTdD2kDI6BLmaeqxRK/D1jbsAQCuRlRva8fv8CgCAq5zB/06JoBW3XjpSguUJobzaBoR4b3PaJBwqqqJ0y8IaCVyQMWvr2nGTw0WUGOGLmGAv/OfJSrCspTZCUqQvfrs8jvadVJfbmq6k9yJVGYjtAjI5LlGbGJEbt9/DrURnj5raUa0ER+RytiqgTQ/3RUVDJ95Zp6aEdYTckDyTW7JLaR2Nnn4T+k0sxrnK8MHGWbx+Pu4qfE8jnKXOlsxHTwEOFOqhqW3j2aujAj3R1TMgajrgQmhbXjw1hG7FM9QRqG3rxpbsUjR29mLrghgc09TbjD4Sbt1tbfNJGC3XXpwU6QdtXTvO3moerO2gxsnyBjR19mJFoiUG/7yuGTIGCBrvjs6eAeju3kdylL8lu5eTU/H8rIm86KqIAE+aaKatb7dJlwE8NB1N8HbDX0vqMF8ViCs191Bu6MDZymakxwbh72dHorSmDbVtD1DR0In2B/1YlRiG3y6famXvFvOrEPNV/jUj5INmJjJ/RCEsTwjBRxdqsDwhhNYkJnPqjLPckW/CnlmIaz4j+RVcBzPxryRH+dHnAwB+tyYBK6aH0bDZ2dH+uFLbjvmqIJyuuAvDvQeobX2A2tZuyGQMPtg4C3Fh3jina6E0Jtx8BsmBPLKQfApjCJraNhoeuX3RZBoC2HS/VzQHgQuhbXnF9DCaTdrVO0Bt9yR0tba128p2bSvG35Zg6h0wW1EV1LZ146PiGkwJ9Yam9h6SIv0wXxWEL641wMyyGDCxuNvZB7mMAQvg18ticbm6DQDg6Sa3so8TRyUAmoj146RwALCi/+CCRFO9900V5qkCcV7XgnmqINS0dgMAkiItIZr/uFiFi7db0dDeg60LYhDs7cFLJiR2f8Bi6xdz+HJzGfbk3UBFQweyzt6m9ncS0+/lLneY2yCEs2G3Yn6l2rZumgTJZcclDn6uf+V2UxcaO3sp6y0Jm50V7Y9r9R1Ym6zAF9cakBLlh7OVzUibbEma+9XSKUhVWnI2UqL8AQC9A2a7XEjOhihLEIekFMYQLle3IX1KEI2BJ1TWa9TheGNVvM3VvK2V/IrpYTRrdmaUP363JoEK9iVTQ+iqnht6KSZsnOW35xLCLYgNps7p1TMUiAn2QvaFOwjzG4f27j70mR5G+iRF+qGysRNvn9FZjYEoI6KYCHncpnmTqNL47EodTt9o5PEpFeubsf9rS22Ic7oWzI72x4WqVsxXBaK334ySmjZMD/dB/rUGuMgYbF0Qgz9dvEOzl8WSCVfPUIg6fLm5DKpgL5yuaMJ8VRAuVLVSH4aXuxx78ytQ09KNzAL9iJtTbPESEQ6njy7UIEOtwKel9bzw1TVJFvK5T6/UY+uCGDqm/GuWMFiSAb8xdRJqWrpx6nojMtThmBHph39+bqpVktmSqSE8R7rYLnQk6LrHMn4QjmaGYZYzDHOTYRgdwzD/LPL7RIZhvmEYRsMwjJZhmBWj2Z+nFVvTlTwn6PRwX3ywcRZlSLXliLPlsMu7aqDO3YqGTqvrCWmfxTAUemp7TvPMgiqsUYej3NABmYzhOU1TlUGICvSy63QkVcS4/SGfifIR9vEZZQDO61owXxWIy9VtmK8KwnldCzp6+jE72h+XBu3h3KplhJb8ZHmjqENbzOFL+pseG4xL1W2I9B+HS7dbafU6AIhX+GKawgc5mnqH1NLOQEgznqoMok5yrlMdeOjgL6xspgV+uLTn3OJNpE0y3+R/CyV702A7TXbpyR05kEeKrluCfYyao5lhGDmASgDPAqgDcBnAT1mWvc455hAADcuymQzDTAPwBcuy0fbalRzN4hgpB509B+VwnZ3D6RNx1MYrvFHV1AUWoDWJyw3topXbHO1O7Dm9ufWWMwuq8FxCCI4ORs6c17XAzUUGGQO4yGWYGOCJ64YO7B6sM829TvegPZw4tG3NJ3HolxvasTe/AmvU4ThRZqnJ7OEqp4qR2Ou5dZQfRQja6k96bBByNAbRusdZRXrax8LKJh4Drb1nZCSfJS64ARX2ggYk8PFDqNE8G4COZdkqlmX7APwZwI8Fx7AAfAb/9gVgGMX+PLUYyRDDkQz3G0oBFS7IKjRDrUC5oRMDZha/WhqLgxtSsCW7FPtP67Bz6eQhhXPa6o/wO5MZNIR0jToc53QtmKbwwTxVIFzkMqxKDEP+jjTsXhmHfadu0WuR3Ydw9WxrPgHwiiP93cwIuMhlcJExMJlZvHNGRxWCsI7ycENHudfnrrZJGC/p97+feFhrm4Qlr5sTibbuPhpuy93NiD0jBwr1yLtq4I0976oBSZG+vOOK9fbrZgtha4ciYeQwmkohHEAt53Pd4HdcvAlgPcMwdQC+APCPo9ifpxYjKci3piuhrWu3MjEkRvjyXl6xamc/fvccXjxyiX4u1jfj/XO3ERXo6fQLzFVwU0J98PrKOHi4yrH/tI5WfFuVGIbNaUorE5Yj84KYQBF+R3YaiRG+1OxhbO9BiI8HDm5IoSa5zWlKfLBxJp1jW4pZuHMh/XxrbSJWJYbBzUWGzgcD1KfywcZZSIr0Q7G+BdPDfXllS0cqHp+rCLm5FKTfVU1dNG9EW9eObQtj8GVZI60g98HGmTTfgtsm934kRvjSimhkfo5rjbhc3eYwD8EWdh3TYkt2Ka+vW7JLseuY9pHmQwIfLqPYNiPyndBW9VMAR1iW/S+GYZ4BkM0wTALLsrxHjmGYVwC8AgATJ0ohagDfVEJeRq6phKyChwNnit6LHVPZ0ImefjOyivSIV/ji5Y9K0N1nwvq5ExGv8HXKdCBM0IpXWATjO2d02H9Ghwx1OKICvWyezxV43GLzQtMFiZkHQAUvMZsIzSPk+9UzFFbXspe4RwS4rfG+tTaR8hhlqBV0HisaOpGhDseJsgbR8T2qDZ2rCLOKLJFFwoTDvKsGm6Y/Z/rAVdCkDWISc8akKGYKbOzoQd+A2epYCSMMlmVH5R+AZwCc5HzeBWCX4JhyAJGcz1UAJthrNyUlhZXAsud1Tax6zyn2vK5J9PNQkVmg4517XtfEJvzrCfanhy7YbJdc879OVtBjDp3VsdGvHWfn/+E0G/XacfbQWR1t+7yuic0s0NFzyd+OxnjorI5V7znF/vLPGjZ6sE1H53D7JDY+lmXZf/7sKvvPn121Ov9nH160OtaZ/g4V53VNbNy/fMm+kHWBjX7tOLv72FVWvecUu/vYVTbuX76k4x7uPbV1TWefm/86WcFGvXac/a+TFcO+nlgbzrRrq59kToT3V4JjAChhnZHdzhw0nH+w7EKqAEwC4AbgKoB4wTFfAtg4+PdUWHwKjL12JaXwEFwBSISI8HdnBZnYSzjlX75w+PKKveB/l3mejXrtOPt3medttu3sC02UzC//fIUnFOwpqZFSlCwrrkxszetQjhUKud3HrrJRrx1nV/7fszzFN9LKyFYfhYqQLArW2VkUOIKYgraltJ09n2VHRlmNRTirFEbNfMSy7ADDMNsBnAQgB/Ahy7LlDMPsGezc5wB+BSCLYZhXYTEtbRzsvAQnwDWVZKgVyCyoQrzC16bJx1Fb3O3+4eJquMpleCXNEoffdL+X0l8QZBXpkVV0m9rj5yoDUW5ox+XqNswepDLIKtJjc5rSypTgbASKyQysUStoZMzmNCXiFb6iZhliwhGGMJLPYjkSjsA1k2nr2iGXgZqWAL7Jzhmzm7CvqcogalqbEuqNMkMHMtQKGtU0EuYiLr6taoFcBl6b5YZ2tHX38fpKTGu/WKQCAKdMf1wMxVwn1i4xH3FNgYDFr3CyvJH3zEkhqSOL0fQpgGXZL2BxIHO/+z+cv68DmDeafXiaQWzDqcpAfH3jrqXyGEeor0oMG9ILw1UyHq4yfDjIR0Ne6ONaI32hSZgiCcucqwzES0cuo6ffTL8jx1zQt+DltBirF/xAod6pfIdDZ6ushIDYuLhtCYXNUBSkcE6IQkuPDUauph67V8aJCn1yrBgPkzBxj9tXkquQo6nH7Gh/FFY203yKkcY8VSD2DvJHce/R7pVxVDnFhXoDAM/J7cg/IoTQx6Kta8esaH+E+Hjw2lgWHyLabmKE70NOqUUqHC6uxvvnbkMuY5xSKhKGj1FVChJGD0KBREI1l0wNoUJd6Bh1hF3HtDiuNWKeMhDa+odRNYSRk+t8zCq6zYvTT1UGITbUG4FebvQ78v9xrdHqBT9cXE0dj86McShCQMzJ+SiCg78jC0dmQRU6HwyItpuqDMLyhFDkaOqRoQ53ateWVaRH7uDxJA9gtIQduSd78yvwVXkjLle38e6jmJOejGsofREqe6Lc31lnKR06nJyFAROLV5/lO8WHqqwkOIakFJ5QCFdiBzek4KUjl5F31UApnYWwl9xFMnyBhyYD4XY/VRmE4PEPE4e4RG0A8LdfzLe6JjH3kLYeZYxDEQJcQT5PGWj1uy3KDTEIw1bTY4NEBSc5trCyCRnqcORq6gGwKKxststSu+/ULSqYibDctjBm1ITd5jQlvipvxKVBMx83AY87zpE0zQxVUZOFCGFbJUl1wnDckTavSZCUwhMLMWHWb2IxYGbx84VKuqpOirQUfnlrbSK1eW9bGAOT2dpeLgwZBMAzQQ1XaNh7we2dLzZGZ4UA6es8ZSCu3GnDluxSnumLZEQ70w53jkimc4ZagT9dvAPvcS6ULtt69csiR2Og4aa25oabmc1VfEP1fziLrCK9ld9HGDI8FNOMszxXtsKFxUDmU/i8jdacSHgIqcjOU4K8qwa4yhl4uMpwuLgaALBtYQzO6Voov0+q0sJzsze/AjcbOnkvPTe3gVswhSRrPUrWNHmRuS8493shxBLjnMl8Jedx+7owLhhmFhgwmbEluxSvfqKhjLLOKjQyR3IZcPRiLdbNicSUUB86l2RTxj3WsmNoRoZagRNljTbnicw7F8JEsEfFxsOXkFVkmTviQ1gUFwxPdxfsXhmHvfkVeL+oyuauTAjh/SH2f5JEZispbSjZyCNdCEiC85BYUp9gECrh2rZuvH1Gh4MbUqCaMB5nK5vx+VUDSqrb8MHPZmL1DAVeOlICfVMnPi2tpzbv5QkhvELzgO0yjc7U67UFW2ystuirbbFhBo13o6U3uW0T6mRyHpepdE/eDfzTEhUuV7fBd5wrLt1uQ4Zagd0rpjk1x1wKasJG+z+X6jA11Bsfnq/GzqWT6RySY7nj3ZymRHKUn93xjjaa7vdgb34FvNzlOK9rgTLYC2cqmvDC3ImUifVO6wP846LJvPPE2G8B6/tD6LLFSrQSDPUZeJTnTYI4nGVJlSqvPcEQkqsBoFEyOZp6zFMG4uPNcwEAr36iQY7GgASFD/RNXUiJ8sd5XTPWzYlEZIAXEiN8kXfVgJPljSNOYDYUGm3h2BxVbBOjshCzW7/6yXc0ukfX1GU3DNJRPx0Rsg1nvKMNskMgVc+4zuXhQGyeueZB4bz8EOdkrEGqvDYGQFZPb5/R4XZTFz44dxvbFynxaWk9rfpV29oNTzc5MguqoAz2QpmhAxH+43Ctvh2L4oLx+VUjogLHYf9pHVQTxltVRLO1OhtKwRNnC74IfxdW4BIWiRdbkYqdl1WkR9bZ28hQh0NTe48WtRGuUp3h67e1k3rU8Y42UqICUKxrps7l6KDx+OxKHW/nVaxvxnsFOvzp2xpM8HG3e2+F8xwR4Gl3Xn6IczLW8IOopyBhdEFs7OvnTMR5fQt6+k3441e3sCw+BHMHI27+9p2BkoitmqFAQrgPqlu6EeE/DmcqmjBPFYhPLtdh28IYvLU2kbfiPlCot2nfHkq9hOHAlv2Z6/MQY14VnkecyrtXxuGPP0ni1T4Qi2RxRKj3pNq5hc7l2tYuGipM/DAkF2WeKtDhveXO8+Hiaiuiuu9zXg4U6rHrmJZ3vWJ9M3Yd0w6JgVWCBVL00RMMYYLPewV6dPeZ0NNvpqybeVctbOREsB06W4Upod642dCJuJDxOKdrEc2GfulICXYu5duYhYR7I5kLILyOrUgYADYjoMTOI+Pg5lPYi+6xFyHzKCGyjxPCREPyed2cSHx+1YgXj1wG8LBeBTfLWuzecp+PzWlKNN3vxXGtEeWGdjqv3+e8JEb44t1vdDS5EuCHU0sYGiTz0RMM4uADAHWkH8oMHQDLotzQQU0nS6aGYMnUEF78e2FlM1Ki/HGl9h7mqwJRUvPQpELMMjsWq5BZUGXXlCJmqhkJ2HIy5l01iJbeFDrDtXXt6DOZkaoMQnKUHxrae9FnMlPzBzFbHCjUW5lQsor02H9ah7kxAThzs4lnBnmUuse2TGvfB94+o6NOZcBiSvJyl6Pc0Inl8aE4r2/BgJnlldXk3lv1RD80dvTQecq/ZkS4nweyzlbD0P4Ab61NhKucwb5Tt/DKghhau5s71tGck8gATyRF+uG41ogcTT3+9p2Bl/kswQLJfDQGQOL/N6VGY/8ZHZZMnQB3VzlSlYFWIX+EF5+YTq4bO/DCnEiUDiqEzIIqXslFLl+RrdKHo1XwhOxEuKGPqUpL6U2u2UcYNmmv9KaYaYsk7BETSlaRHr/Pr4CLnMEvFqmGbQYZbdPaUHFk02wrp/LmNCVeWRCDw8XV8HCV0VBmbp/Jvb1W345cTT2dp8QIX3x+1QgXOYPjWiP2nbqJzIIqq0p4XIz2nKQqg7ApNRo9/Wb09JuxKTVaUgjDhBR99ISDvFzpsUHI1RisMmOJ0EwcdBiT/+Uy0AQ28nnfqVvUUUgUgK1IG6GpRviZRJuQa3L5f8j3jqJOHF3DHghlB7eMJWBNikds6f0mMwZMLNxcZHj/ZzN5ppLhRMjYioL6oYCMG4CVyYXsErnzviW7FAMmMy2IQkxN9iKOxK45WnPCvY/c/v2Q5vxxQ4o+GgPgCkmTGUifEkRNPqnKIBo5RFZpxCTTZzJjT94Nq88HN6TwzDIyGfDh+WrRiBJHceQkkmf2JH/sybsBmcySM0A+c81QtuBMtBGB0Dzh6SbHp6V1uFDVSqNjhOYvco2ePhPO61tgZoFt6Uq0dvfTtogZxJapw5ZZRFvXjqmh3iNuWntUEJNZZWMnYkO88atBQZ531YDYEG/EhnjjurETb6yaRuc5/5oRzygD0Ntvxu3mbgyYWaycHoaqpvvI0RjsRmJx54drkkqe6IdfLokdkTFxFdz7P5uJHyeF47jWiONaI5Ii/R5LbsgPEc6ajySl8ATjjb+VIczXA8lR/shQRyAlKgAyGXD4fDXaH/SjsaOXrtKJcK1o6MS+U5W8VZSYgJcN7hyEioK8+I7s6/nXjJg9yR/7Tt3CjAhfHL1Yi5Qofxy7Uo8di1VoaO91ypbMFSTLE0J4ZhCxxDXSvxPlRhTebIarnIGm9h6+uNaAzPXJotFKv8u/ARkDuMgZaOvb8YwyAHvybqC6pQuebnLUtnVThVLb1s1TDrbCWGdP8repUB8n+kxm7PuqEpWN96lC2JJdSj9vmjcJa9ThVqG6b+SWo/7eA7jIGcgYoNzQAf3d+zi4IQWt3f1U2XPnIf+akS5Ipkf4oratm8518/2+ERPY+deMCPFxx6+WTkGqMoj6GACgd8Ashb0OwlmlIEUfPcF4ZUEMtmSX4nJ1KTUB7D+tAwBqq+WaW0hS23wVnz+ImHO435nMsMnJwzUP2UpGIsIgJcoP53QtiPT3wDldM+arArHv1C27kU0EBwr1kMss0UYZagVyNQZ4uslRf68H81SBvLoGgIWGWUhx3fnAolAGTCzKDfwxZhXp8Z8nKuE+SBMOPGSb3bFYhf2ndcjR1FNTBHc+CcSisGyV8/whmJBSlZZym1uyS0WjjmxhwGSx1a9RK/D1jbtgWRMYxlJxNzHCFy8dKcHaZAX173BJ/Ugt5X6TGa7yh3M9UnPyKBxZEqwh7RSeYDiKuuCaXyoaOpCrMWC+Kgjndc3wcpdbEppEoooAx5E2jhK9yG7j6MVaRPqPQ21bDyL9x6Hc0IHnZ0Xgfy7VOaSx0NS24ff5FXhGGYDM9TPR1NmDjy/WYryHC46V1mPn0slYkxRBz3/12VgEeLrhows1yFCHY+GUCTSh6mpdO87dakFKtD+95q7PyjBPFYh/XR1vtcL0GeeGBIUPzutbIGOAEG8Par4i0U3cKKyKhg58dKEGm9MmwUUuG1GKho2HL6Hpfg9Soh6em1Wkx9tndFijDnd4vtDEFRngiVuN91FW32EVdSSG/GtGhPl6wMwCZ281Y+uCGPzj4slgGKCysRO9A2bMUwXi3W+qALD46EINVYyzJ/mjob0XPf0m3G7upteSaCu+f0jmozECrk1c7AUnAouwdR76h5nwcpdjb34Falq6kVmgH9ZqzZG9v1jfjD15N5ASZQmVtSiGB5ivCsSFqlZsWxiDXZ+VQd/USVfVU0K9sSW7FLWt3VgyNQSfltahqqkLxvYe9PSZcExjQEqUH0pq7iFDHY5PS+t51wZAlcAxTT2Oa43IXJ+M52dNhHqiH74sa8BxrZHy8xzckILtiyZb+RiWTA1Bn8lM2yq904bz+hbqGxAqxIcZ0wrkafmmOG67wxV+XO6ilKgAmmcQE+yFYG/7mceAtQLPKtLj0NnbcHORwXXQZGbPlDMzOgCebnJ8fPEOZAygrW/H6hkKTAn1xttndNg0LxprkiJQUt2Kc7oWRAWMQ2FlM08xnChrtDKlSRnN3y+kkNQxgmJ9s82wQvL7ibJGZKgVtKLX5jQl1qjDkaOpF80Kdhb2sotJCGxpjSUXoq7tAeargmgIrMkMLE8IQY7GgPRY/lafhDmeLG/E+z+bSUNu02ODcN3YiR2LVCisbKJ1DQjNN+GBmqsMxKrEMN4ckPDdxHBfm9nQ3OOJaaPpfi9kDMOb27yrBiRFWhTiq598RxPBpoT6jEo27+Y0JXavjMPv8yuw7P87SxPRXk6LcYqdlGvievUTDX6fXwFXOYPfLIvlmc2yivSiGcCk3YMbUujxLx65TLOYSYTSdaMlIbKi8T68BqlVuGHQ2vp2WkCIG5oqZR3/sCAphScY3KiLDzfO4r3gXArpDzbOxB9/oqaCIatIj8LKpmHnF3ApqrlUB0Q4ARY7c2ZBFXYunYzrxk7sXhmH0po2rE22ZE/LZUBhZTPmq4KQozFg/fvfUsFDlAApC8n1KWxbGIOdS6dg28IY5A7ufv508Q7yrhrwzjo1Vs9QYPtRDVbPUPCqxRFBeaOh0+G4iR0csCgouYzBr5bGYlViGLZklyJXU29hTI0NshAPqoLwZVkj9bHYopx+FJjMQKT/ONxs6ERsqDc2pylRbmhH74CZKlF7NnqiwHM0BsQrvPHb5VOQWVAFwBKSOivaEhQgljfAzeTm5gMkhvMz4J9LCMHd+33wHeeC2rYH8PFwwXldC7YtjMH+0zrIGPAoRuzlKgyXPl3Co0PKU3iCcaBQj5qWLqyeoeCZbfKuGhAV6AUAVs5gboGZ4TKh2otxJ/4MsTwFcu21yQocu2LAzqWTkVlQBYWfB8rqOzBfFYif/y8VdUrKGIb6SEguhXD1KSwWxB3TUBlWxeaX67Dn1r6OCfYaUdZRR3g9R4uPqX/mARIUPig3dPCc6fZyBezNyVDyBsTa0da1o7a1i9aa+PyqEV29AzCzQILCBzWt3QBgVcTJ3jVJngl5nsgztyoxjNb4kDA0OJunICmFpxS2ooMOna3CKwtiHFIYO4oucjY5TAiyqlyeEMKzO08L88Y5XctgoSA5XeVzBcNQxiBMuntU6mZueySaiERzkdrKoxVdxOUq+rKsEV5uMtS29b20JCUAACAASURBVCBB4YPdK6c6FLL2kgBtJZ+JzZe9BQVXMZAIJaIYXOUMPnpxNm3LHvU4VxGThcaSqROQf60B7i4yq0WHRMXtPJxVCpL56CkD2XZzaQUIY+T2oxorYQqIV/pyREvw1tpEauYh9nlbjKrCa21Om0R9CcTJ/KeX5yJe4YN+E4slUydQ88uOxSr6N9kFzY2xLgPKvbYY/cajVDgTtpd31TDIIdVE/RtirKsjhfO6FotPISMR6bHBg5FcHrjR0OkUO6ktIr+8qwabNCVi95+EEouFKW9NVyIywAvTFN7I0RiwZGoI3Fws4sVkZikxoyNqFHJdwLKz6Ok3IUdjAMuyvMXBD41K5GmClKfwlIFrShHGhw8l7d8RCyrZKXDZSgHndgpEKGQV3aZCpljfDGN7Dy1fOU3hQ8n+Dm5IETVZ2WrfFsPqcFbxYu0JTWVN93ux/7SO16+RXLUe2TQbgGWlnsvZmaTHBuFydRs9zhZrq60+cAsqCedJ7P6LcRtx8wEIW6mbiww5mnp4usmpLyhXU4+YYC+H+Rvc66bHBqPfZLFkyGWM1XVHi6V3rEPaKTxl4L4s3+pb0G8aPkFYqjIIcYNUDWQ3cKBQjzXvnsNfS+oAAHOVgdi2MAb/8MElbDp82e5KjStgdy6dgg82zkRmQRWyivT0+z/+RE2/37FYBcAS6UISrRwpNnv01sOBWHurEsOwKjGMfkeq3nFXwyO9aiUrdW5diO9q23k7KdI/oRlQzGF76Kzjmsy2osvE2swq0uOlIyU4uCEF81WWBcKAyYyvb9zF7pVxcJHLcFxrdHjNA4V65F01YFqYN3I09XCTM8hQKzBgYvHC+xfxes7DYIZUZRCmhfk4jCaTMDRIeQpPIbjUEDIG+MX/UtH4cCFNgz0U65vxwbnbkDFASU0bPFxlkDHAX0ssWb4ZagXe+qIC53XN6DezSI8NwnZBnV8ubPEl5WoMPK4d8n1Dey9NICM5GFxeIm4/yZhGusIXtz2SBLZp3iQsmRpCr62ta8fLaZPw9hmdXY4mR/TR9n5v7OjFKwtisCYpwmqOANg8j0szwU0UTIr0w5RQb945tW3daOx4SD9SrBevMifMe9h1TIv3z1Xj18tisSYpAsb2HgAsbrd0Y2aUP/ZmJCIp0g8+Hq7IUEfYvTd9JjP+40QFqpq74SIDXF3k0N29Dxe5DCzL4mpdB5o6e7B4aghez9Ei/1oDEhQ++PrGXSRHSTxH9iDlKYxhkNwFNzmDfhML73Eu1JS0JbvUqRWsMDZdxgC/z6/AvlOVeH1lHGQM8PHFWvSbzOgzschQh+PDjbPttmnLrn9k02zR7xMjfK1yMOQyjLgt2dnwR3t2bHs5G86cL/z9QKGe7qASI3zp6p/bJ7IrcNQvMQp0Erprqy/kM8n7IG3sOqZFuaGdUorsO3UTuZp6AEC8wnKuXGbxg2Sow1HR0IlifbPT/hsAcJHL4CJjIJPJ6E731Wcn479fmgN1pC8+vliLGW+exMcXa7E4LhiG9h7sXDqZl7MhYfiQdgpPEQ4U6qGpbcOevBtYMT0Urz0XBw9XGd79pgpmlkVVcxdWJYZh07xJDtviruojAzzRN2DGxdutAIAZEX4ovdOGATMLFsDsaH9oau/xSN8etaiKLebLi7dbsX2RklcQ6FFtyc7UZgYsq9rqli7s+6qSZkWTqKv3CnQOGUMdZYFzfx/v7oKss7etqDzE2GWdaXeo9a7J/Z8S6o3tRzVYMT0MK6aHoehWE7Iv3MGPkhSIC7GYFl3kDH69LBZ78m6goqEDWWdvY/fKOOxeMc2KSNER8q8Z8XLaJIT5eqB4kLk2Q61AgJc7MtQR+PvZE3HimhH17T0ALMR6meuTMcHbA8e1RsSGeNNdnAQ+JJqLMYg+kxm7PivDjsUq/Hb5VNS2dWP/aR1Sovxw6vpdbF0QgzdWxTvVFtdswjUjaGrv4by+BSxYsIPhhvce9FNBbcvEYE+oicEe86XPOLcRpaUeCkW3p5scOZp6SnsREeBJWUYJpYY9QchnfQ3F5rQY+huXcpvwNwmpPGwpP3tV8GyZgeydQ+4/d24CPN1w6vpd7Fw6GftP63C5uhUucgauchlenD+JwzulwO4V03hzK8ZxJLZw6DOZkXfVgL+U1FHm2pqWbrycNolSdORqDFD4eqCzdwBmsxlhvuPwhxM3kbk+2akFz1iFpBTGICIDPJEc5cdbRW9bGIMvRXhn7IH7shJhvm1hDOrvPcDNhk7LDoEFXpgTid8uj+Ot4Akl9lAErdi1M9QRWDI1hEfDTL7j8hKNFC21s6VFa9u68cW1BsgYoPROG7641oBViWFUeQEPacO59OBkDKTvyxNCkauptyIm5FJu52mNSI8NokR79pSfLcHPde4LFVZtW7dT8yikL184ZQJyNPXo6Tdj9YwwTFP4YN9Xlfiu9h62LohBntbIa8uWT4e7cMi/ZoSmtg1v5JbjuqEDMhmDV5+djFAfD1Q23sdxrRH1bd1475sqrJsTiYqG+0iJ8kN1ywNcvN36g6pZ8UOFpBTGKIQr0U9L623WNCYQrtj6TGZKTNc7YMbsSf7Yf1qHysb7+FGSAhO83TEpyAtfXb8L1QQvvDjfsjrzGeeGxMEXnCiGodZwJoLCVi2D9wp0dus0OwOxFWpWkR7vfqPHzxcqbQpIImAz1ycjxNuDMqj+dnkcT9lxixhxBbPJbMbBs1XIXJ+MzWkx8HKX4/f5FfjqegP+UlLHo9x+ftZEXLzdQnmriKAVCxSwJ/i1de2izv1/P1GB7As1vN0Nl4xQOG6LIgtBrsaAqqb7aH8wgJXTQ5GrMSAxwhc3G+9jVWIY3lgV7/Q94ZvL5Mg6exvzVIGYGR2A5QmhyCyowqvPxmL1DAWuGztQVt+BHyWF4csySyitp7sc5241g2GA68YOSuo3FDPlWILkaB6j4OYBnChrwLaFMXZDAAFrByjBca0R3b0DyCyowqrEMBzckIK31ibiw42z8eHG2di5dDL2nboFwJLMJqyF7ChRSQykj8e1Rh7pGgCa1PSoIafC8ZKM4Z1LJzuVBAaAjstVLqOhqMIxCJ27IT4evOPiFb5wc5Gh3NCJ9XMmwmR+OLZifTMuV7cNOtnldgMF7IXh2nLuTwvzoZ8PFOpRbrDM33VjB014JA74LdmlWBYfgj/+RI11cyKhqW3HjAhfFFY2Y/fKOBy7YsCOxSpKP2Er1FTMma+ta6ecTGvUCnxX2447LRazJ7eNQC83tHT1WarjLYzBG7ll+H1+BZIifTHe3QUTvN3x8kclePHIJfoMSlxJw8Oo0lwwDLMcwP8FIAfwPsuyfxA55nkAbwJgAVxlWXadvTYlmgvbsEdnMBxOG2fq75Lz4kK9ca2+3SZXzVA5lggVgoerDK+kxYxKTV8yXpJEJ6zqJpZ8NpQ5FqNz4F73cHE1APCoQri5ALZ4l0aK+4f0hVuUKF7hy0vOE1KNHCjU48+X7qC6pZuOi8u3ZSvCSDhPNPlxsB40Gd/EAE+UGzrgISh8ZDKzYFkWDMPAzLLoN7EwmVm8Ptjnl45cxoN+M1xkDP77pdl03qSEtod47DQXDMPIAbwL4DkA0wD8lGGYaYJjJgPYBWAey7LxAH45Wv0ZC3iUxC1hSCUAu6t8svIj5xXrW9DVO4D3iyzMm2T1HBPshQOF+iH1hbvDADDs5CR7oabc8W5Om2RFZmcrhNLZOba1S+Jet99kxsENKaK7E7LC5x4/PdwXq2coeOPiruiHuiombedo6rFGHY7Mgip8q2+hv3+rb8HJ8kYreomWrj4rmvaT5Y12w4KFu6fjWiMGTGa6I5irDMSAyYzrhg5kqMPhKpdZdolndDCZWchlDJYnhKK7z4SefjMYAK+vjKN9JktbFzmDb/UtkkJ4BIzaToFhmGcAvMmy7LLBz7sAgGXZtzjH/AeASpZl33e2XWmnMDoQW8FyV/3Cl4zrgM4sqBqkkbYoAkLDICxLOZR+EDONI5oOe8RoYuypLx0pwc6lk1HV1EUJ/cgKfPUMxYhQU9jbTQAQ3Vlx+829vth9IavrpEhfXK5uo5+dISXkztfDnUIQTpQ10voWXGUststx9t6Igbt7arrfy7sHZGyE+fbFI5fR02+Gh6sMzyWEIUdTDznDwMSycHOR4cimWXQ3S3aUh4qq0NNvtru7Hat47DsFAOEAajmf6wa/4yIWQCzDMOcZhvl20NxkBYZhXmEYpoRhmJKmpqZR6u6Tj+Fy0AvpJ7gFagDbFAjbFsZgb34F0mODUVjZjNdXxlHeG5WA58ZZcO323OS5VYlhonZ+R4lby+JDsCW7lNr2LX6QSnw6SNPhPc5C/5WrqXcqsc/eHJPfuLsJ7mdS2+GddWoc3TzXUvmN03fh7kR4X7h1uNNjg/FNRRN6+k10tU3mzN4YyHyRxDgLsV8z1iYraH2Kw8XVOFxcbbXLcVRXwRGEu6fVMxS06FFiuC8ObkjB5jRLMl7eVQNc5TKkKgMhY5hBhQCYWBbzVYFwd5Hh5Y9KkFVUBQ9XGVzlMniPc4GrXLzYlATnMZpKgRH5TrgtcQEwGcBCAD8F8D7DMH5WJ7HsIZZlZ7IsOzM4OHjEO/q0YLjMkUKTyFtrE2kNAwIxc4rJDF4Ft3iFL9xdZIj0H4dLg0Vohrp9J2YToQB6a22iqJnGllOXy0vUbzJTE1S8whcm1mJmWDJ1AvbmV2DJ1BC4yGU8PiNbsDfHXEc7dyVOPkcFeg3JvCdmqiLV43I09ZgV7Y9+E4t+k9lpkwm55r5Tt5AeG0zrUnxZ1ojdK+Pg4SoHgP+fvXePiupK08afU8VNFJCbQAGCFIIKIiWiNkqwNe3dtGQmyYyJv9EoMfa4/DrOdGdiOt98K71i36bNt+xkUGmVtB2/TjoRu5X2kmhAlEQBK5agiFUlCBRgURJEkFvV+f1RtTf7nDp1U0lPtJ61siJVp87Ze59z3r33877v82JgyIKgMT50Uq3UddL6GJU6q+ouMfDXbFnLYohpLjIJBfqPZNhfumV1pmtaR5RwNx2soX6MLQuTIZdxkHGAmQfmJ4fjalsPnpkRg75BM8b4yfFvi1OwdVEydpTWY+uiZPzb4hQkRY6VXER44Rp/b/poN4CveJ4vtv19GsB/8Dxf5ei8XvrIOaQcxqPBq+4u19GiN4TaGDZb8D1lOL5u7kZeSgSOqA2jWnyGpUIILZGvikVqdJBARpvQHAAo1UFoh9mJobjY2OUR3eBsjB053sl37lI77LXIb9gxJ87heckR+Epv1YbypA8sjRPo70OvS9pQZ+imUtm7TmuRnRiKr5u7aRU1wDm9yI4FmfhI21n6idB2LBVFPiPtkcuA35xowOSosWjrHsDmBUk4rzVhXnI49MZenKzrwJK0KCRFjqV/u1vf40nC/wT6qArAZI7jJnEc5wfgnwD8VXTMEQDfBwCO4yJgpZP0o9imxx7uaPA8CshlwI7SeloeMzsxFPeHLPhSZ6Jqp9tXTMHOUzdGbbXGUiHWkp2xOKJuhdz2VLP6Ta/kJqF/yIIhswV1hm5a4rOqsQv5qliPypI6G2PW8U4mIrYt7lA77Op6Q3E17Q8Zc1V8CE7UtmPNnHic13ZaM39lHN77QitQEa3UdeLl4otYd+Ci4DpiGofdAZCdWkGukoYcPz01Cl/UG2n9CxKeTIw2AMGOh3Xmk10cCW0mEwepmf2LZzMEVNT02BD6GRmTwjI9il/ORunWp/DeGhUKy6xFlgpylXQHebKuAz33hwVy4J7oLXkxglGbFHieHwawBcBJANcAfMzzfB3HcW9zHPeM7bCTAEwcx10F8AWAn/A8b5I+oxfuwFHUizv+BneOWXfgIooqdDBbgO226I+Xfn8BZxs6sXBKJFarYukLXZCrxL51sx6q+IyrCKLNC5LwTmk9psUEo7zBSNtUVKGj8tDASCSVjOPwmxMNlEvfvmIKdYpL0Q2OZKKLKm5KRma9cVhD+XgSQfPaR19jQ3G1010ba1wJFbbpYA2+pwxDYZkelbpOmC3AmjnxOF1vRFZCKP56uQ0BvjL8Q1YcXl+WCo7j8OGFZrxZokGlrhMbP6jGmXoj5tmkrMnYSQndyWUQCMq9cViDXae1NuezlaoqURswNTrIzmiT/JFX85R2kx87gbJ0Iivwxz6z9SIqyp1or29rIfSkYFST13ie/xvP8yk8zyt5nn/H9tn/5nn+r7Z/8zzPb+N5fhrP89N5nv/TaLbncYfYMcmGObpahQJAk6kXGz+opoaBUC9Npl56zLzkcOworYdcZjX602KCcU7biblJ4di/brZdDP3DrtZctdtsAeYnR+CcthMxIQEoyFVi84Ik/Or4dWoMyZgY7w0gJiQA/r4y6I29eG+NCmmKECxJi6KJY64S+5wlulXqOmlhoLnKcOxZm4WBYQtK1K1Ymh7t1FixxpUNWd2YmyRYbR+v7UC+KhbntJ3IiA3BvnXZ+MWzGSjIVeLghtnw95HhwwvN+P/2XUTfoBlviug7YmSJSipg7bfeaL3HxzRtNGS0f9CME7XtyFfF4mJjF3xkHC7d6hI8H3vP6umESnw6pAodUXolu7IjagPeLNHQiWDLITUu3jTZVY9jJyd3KuY9SJKkF47hlbl4jMAqm+4u1yEqOADLp8dA09KNfFUcZDLgjU9rMThsVfgkoYxEjuBm5z2cqTdCb+zFwJAZPy+9BgD4t8WpVK4gKyEMY/3l2FFaj+O1bVDf+gbzkyNwte2uxxpE7iipEimEDcXV0Bl7aPGdwjI9ZDKgsbMXFTesRqC9ux9Npl58UtOKH/9gMgrL9DBbeLz2gxQAwM7PGtA7aMbWRckIHuOHqOAAbDmkxms/SEG+Kg7xYYGobuoStCk+LFAwbu9/ocNPl6VSQ8sKvnXcHcDG3EnU4PYODKPOcBdpimDUGpyPDznPpoM1AqE5wq2z2kPlDZ0oyJ2EM9eNWD49RtDWIZuarYW3qtf+8h9mCK7jSOiuRG1A4UszETnOH7vOaLFiegz0nfdwf8iCGx098JFzGLbwkHGAvrMXrV19+NXx63ghO46GJH/wZROWpkfhk5pWrJ+XiM57A3QC3b58Gow9/fjwQjMAHh982YT31qjQ+s19NHTcw6oZCsSHBaK5q88jtVN2IfSgsidPCrzaR08AxEZ1VmIY1cYhK8/l02OQr7JKL7999BqenjqBCqxtWThZIFq3/3wjvqcMQ2NnLyr1dyDjrLLVAASGOishDMdr23C9/R5So8eh5EfzHuhl9ESyWmfsQYnagKXpUdi+fBpkNn69vr0HPnIZ9v3LLFh4XnDM9LgQ/O6MFmGBflRFc9UMBd4+eg1To4MkRfqk2sSO248WKO0KCRHBN9bg1rffRYktxPMPG+bQ8SGaTlITYVRwABWa+9ECJTYvUGLLITVkMmD/+UaqPbRt8WRsWTiZyng33+mjOlFv/aXOmuzFAS3f9NOCNOJrSWlTxYUFUoG8j6tbsG1xCu70DqKjZwAcB/jKOIDjkDVxPEqvtON12+QokwFFZ29idmIoztQbqdx3VWMX8lKtfojegWEcVhuQlTAe57QmqoX19NQoZMaPFwgneqJ26qhwk5Qq65MO76TwBMCZUdW0dGP2pFCBYuqy9CgcvdwuEH1jV6EFuZOgnDAOZdetK28fOYfkCeME4m6AlUIpURuQGj0ODe33MNZfjtWZcR6/jO4qqVbqOlFYpqfJVdWNXThe246kyLFovzuAFdOjoTf24nS9EUvTo3CitgMzE8ZD09KNcf5ygcpoc1cfbhp78cmlVkmRPqk2bV6QhE9qWt1WZa3UWUtnrpoRg/KGTjrO0+NC0NDRg9+d0Ures6OXDdAZe/HqU1ZZj+XTY5AQHoidp25gz9osmC2gRnZ6XAhSo4NwTNOGkDG+2HNWj09qWjBo5rFoSiRu3xsEb+HxdUs3nRjEk26lbkRZ9UBlI45p2qhA3syE8XjrSB267w8hMy4ELV39kMk4qOLH40v9HSqPTSbNpenROF1/m8p9T48LQb4qDlkJYYJdTnXTNx5JeLvCo6609zjD3UnB59tojBejA01LN+VzSXgk4XPJTiEvJRK7zmgxPzkchy400xBRUjR984IkysceqGyE2cJjjK8MPACzhceO0npsXzGiAko4dcJVk78Bq4/BUycf6yTcujBZckJgwx2NPYM4p+1EuiIYhu5+ylUH+Mqwb102cpQReG6W9TfLbCvrfJUCf7xwC0FjfGg4JeGf5yrDJTlr0qZ8m/wDub5UsXmp9pIi9+L25ygjsGqGfUgrAIGkBLnOkrQoei5yvTRFCA21JAlwE8PGoM7Qg/TYYKibu2mi229O1OPj6haEj/UXhM/ajastu1gMEoqar1KgRG3AxcYuzE4MRXlDJw2XJZnrZEzJM0iuw/oUxM8fGwiQowzHgcpGwT1xFcbrxaOHd6fwHQaRaM5LibTxudGUz81RRgi29V/q7+DFOfH4yZKpAEa4crIKfT57Iqob76Ch4x5eX5aK6YoQfKm/A185h5iQMZTf/d0ZLV6cO5Fy6sTHcF5rwmqVOGFdGuJ6DaS+wMfVLXZ1dll6oFLXiQ++bIIycixqDXeRlTAe1U3fQDE+APcGzHgmU0FXnrfu9OLj6hb8x/Ip+KSmFTEhAfj0kjUrtuhfZuH57InUVyC+JruC/ri6BVsXJdvVRt57Vo8Jwf52v9t7Vi9Zb5rdQUmtjB3RIAPDFqd1jcm5jmnabRXwugXV1f5p9kSrxpBoFS6+HqFxSDvf+kstshNDcU5rooa74kYnOAB3+4dpUSWzxUL9EYTTf+tIHaKC/RHoJ6eGX7zLITuno5cNVAp99qQwHNO04ZimDZnx4wWy6Y52Zg9b4e9Jgpc+egJADPsIn3tbUL7Ruq2Pwul6I/JVCnx+zSigPqoau/DKUyPS2m3d/chLjYDe2EtLS9Ya7gqcfqtVschKEL5sWQlhbk8IwAjtJZMBbx+9RumZrYuSBdXbgBF6gF3ZhgT6oaWrD5rWu1iaHo0tC5Ntxe370dbdj0GzBb86fh3/viQFBblKm9Fsg4wDcidHYMvCyXR8ti5KtiuGwzouSdEicdGYCcH+2FBcDX9fTlAk59mZsYLzkePF9Q/ExW2Io1t8f10ZtpFJNQpn6o3IV8UKitxIXYuc19n1JgT741fHrcEIE4IC8PonV+An5/DTZamIDPLHgfNNeGZGDJru9KGrbwirZijQ3NWHo5cNaOi4h5SoIAwMW+jEMysxDFkJYWg09aKho4cWUNLevkcLEuWr4pAZPx7HNG24ZerDH2zOaGe7z4et8PckwTspPAGQMvyEz917Vo+ZE8fj82tGWsVrWXoUyq4bqQEXG4ZZiWEYGLYIithkxo8XcOCPAmQV/MantXh66gRaCMiZX4Jd2Q6aLfi4ugXgeWhv38NqVSySIsfiD1/eAols2bM2C6sz4wTFc+ra7sLwTb+gvvLqzDjBtdx1XMaHBcLfl8OO0no0mXqpXERhmd6pQXpU0TJsPW4yqT47U4Gjl9vo5Eom3Qe5FlvF76axF529gyj6l1lYnWmtfufvy2H/uSb8+h8zaEZyiboVOmMv9qzNwvp5kyQnnkA/ueB5Ehckig8LRP+g2c7n42hHQIoIufJLeeH+pDCq9RRGA16ZixGIpSaKKm7i2ZkKxIeNpdmva+bEIz5sLC7eNOGLeqNAdkKKr3UltfAo2iyWpkhTBCMjLkSQ4+Doms6UOmfEheCc1oR8lQLvvqCi/g7SZ2uGb9UjVdF87SM1StQGzE4MhdaW++DMILkzvu7KXRAZCrMFaL7Ti8OXDPTvDNvCgN0JOhtXR5CqCcGei/hGnKmTivtD8l+mx4agvr1HUiZkanQQNDaZECKRseu0ltaTKKrQYeepG9Tf4qydXljxP0HmwotRBpEAIIk/2xZPxqELzZDLQLNfD11oxvX2HlrFiyQpOZJdcCdZ6GHbLJSmUOCq4S6OqFvdEvJzpNQ5MSwQV9t6qEPztY++xq9PXMeaOfGC5C0OQEJYoMdJTo4ym0/UdlD9pLyUSJcrVHfG15WwIbnv2xZbczGut9/FoQvNtFAQuUbx+tke3UtxHyt1nThQ2YgcZbjkeLEOeQAOk8ekKvsNmS2o1JkEGcjshP+vNvnuTQdr6IQwbLbgmKYNr330Nc1/0LR002eJXL+oQuetuPYQ8NJH32GIqQ7i9P3V8esIH+eHz68ZsTQ9GiXqVqyYHo302BCcrOtwuc0eTecd8YPsKK3H0vRolDd0Ytviyahq7KJZtGzbnOViDJot1EldcaOTJkk1mfpQom5FjjIcX+rv0NrGmw7WQCbj8H//KRPLp8d4RNuIuWuyC3k+Ow4XG7uwND0aR9StGOsvt/O5PMgYOQvVJd+/ffQakiPHUupw+/JpLs7sfh/JeAHAz1enS45Xpa4TPy+9BhlnFRrcvEApeZy4Pz8vvQa5jKOht+RY9nmODwuk/oW+QTNu9wxAJuMwNTqIhr5uXz4N6uYuOkFsWTiZPlt5qREPfR8eN3h9Ck8ApDjbrIQw3Oy8hxJbolB1UxdNenomU4EpUUEu48EfpfNOaoL5pKYFFh44e6OTJtFlxo/HLVMfPrnUiqXpUXR1P2i2YNPBGpysa0dCeCCNSJk9KRRvHanD8unRUE4YRyNbZDLgk5pWOlEQfl1v7MXtngEa8ulpkpPYsL3/hQ7PZ8fRIvIFuUkY6y/HzlM37KKZHI2Ds4mWjVBamh6Ngtwkwe+sORg+OF1/G3GhY3CtrYcWrnd1bnf6eNPYiw4n40WeieXTo/HTpVNoFvfy6TE0i17sgyH9IUmRYj+H2NnO+hdefSoJEeP8cLreaIuy+gbT40LQ3j0gSJDbf76RUmje6CMhvHkKTygqdZ0ob+jE/ORwnNOaMD85nAq/OYvRZ3lfIjq26WANMmJDcE3E+3oCqepnzrI3mwAAIABJREFUJB6ebQcAXGsfoX+mKYIFtE9Dew9eLq6i0tek3CeRgSBgeWaSr0CqwonzINjYf3cgzqkgdQFYAUCSQyA+r9Q4sL4RMUbi+63Kr9MUQdQvMpKD0Yp8VSw+v9aBYdvkSfITnJ3bkz46Gi+iecTeo80LkrD3rF6StiL9SQgLxO2efsE5ye8c/WbrwmQU2XwW+SoFyhs6aX4OGc+e+8O0zaMl1f6kwOWkwHHcFgAf8jzf9S20x4uHgFi/HoDN8RqLNIWVj16ZEYNti1PtkrDERguw8r7ndSbJpDJ3wconS5X6nKsMFxSKz1FGYJoiGDtK63HV0IPyBqOg/gEgXTsYsPpRyIRArr15QRJ2nrrhNFnNXYiF16QmSkcTjXgcnNW6IEXtR8YjCDtK63FEbUBbdz/tE+tA33SwBsNmC947o7Vz3j5MHx2N1ytPWY1ymmKkoBBbK0F8TnFwADuBSf1O/JvdZ60V1p6bFY/nZsULhPeOXjbgmKbNbpHhTXp7MLizU4gGUMVx3CUA+wGc5L9rIUtPCMQyw3vP6pGvisWJ2nYE+Arr6LISxOzugDXevnJr3VtPjKlU5AwATI220lbzlOH4V2aSyVFG0PKf7Ir7qqEHJepWWi+YGKq9FXqH2c9SkUqFZXo6UbjKRnYGsaEmkxmJhnEHrrK3HSFNEQIfOYc6w11sXZgsOfntWZuFNz69gkrRJE7URqUiuwBIRgU5WjiI++LuJMc+l7vLddi6KNlaQtQ2gbEZ0OLfaFq60WTqRfH6bPp5RtyIsm1GXAje/8K6WJirDLdbZHjhOVxGH/E8/zNYS2buA7AOwA2O43ZwHOedgv+HgUSdsKusd1/IxL51s3CyrsPueHEkCmu0hswW7FmbZScP7armglTkzKaDNdC0dtPSjWL84tkMO6NV3mCk0htEWnmuMtyjGrzuaPF/m3BX4pmUQ2XrKgT4ygVFcaSM752+QcHYvHFYg798bV1Fs9LfG4qrBeVDyXcsJQe4Hi/2eXFWx4CNuCKRU09PjUKlzoS8lEjsPHVDIOFOjiMTAHl22R3tqhkKvJqnpMV6yHh9pbOWY3GntKoXDsDzvFv/AZgB4P8CqAdQCEAN4Nfu/v5R/ZeVlcV74RyFZVr+vNYo+Oy81sgXlmmd/u681sir3j7F//PeL/n0/zwhOAf5PTmGfHdea+Sn/Ow4v/esVnBs+n+e4NfYzsOeS/x7R20g3z/zXgU/5Wd/4/ee1dLP1x+4wC/8ry/s2uGqf55Aagz3ntXyU352nP/tyXqnfXCnX67Gged5/rcn6/mE14/xqT/7m9PfsZ+RsU/92d/4qW8d59P/8wQdux//6RKf+Poxu3ulevvUQ/XJ09/uPavlE18/xj9XeJ5PfP0Yv/3wZadj4+51yHj99mS92314kgCgmnfDxrrcKXAct5XjuBoAvwZwHsB0nuc3A8gC8A+jNFd58RAgKzN2VU92BeJKagTs7uJQwVy68hL/nqUNSFEVEi8vjkM/rzNhemyIQ9pKjN3lOhy9bBCs7ldmxMDCW4u/EH65qrELt3sGKO3gLK/hQSG14yFKrWRlTK4tHkep8fV010J2FTm23RGB1O8c5W5kxo/HnrVZKCzTIzlyLErUBqxWKQSOWHdX+1Ltc1TQydXvCsv0WK1S4GJjF1arFDhe22FXqEc8Vq7a6C208+jgTvJaBIBneZ5fwvP8n3meHwIAnuctAFaOauu8eCi4SoJi4YnREr+kBblKLEmLwqaDNVhT9BU2flANX7kM+apYqG99gzpDt8BQimkrMnmxVAGhPwrL9LRA+1c6E7YcstZc3rooGTtP3cDZBiOll1hunFzPnRKjUpCa/EgJT2J45DK4Pb6eJAW6mqAJrUL6x/6bTTirt1F1eSmRduqm7LUexJg+KDVHVFXJOJJIIrNlpOLc1OggQX9In+bZVFSJn0RqvDyZoLyQhlfm4jEHK0XgzBn4sOcEQCUkAODFOfE4XtuBZelRAsluFqxBIy91naEbvznRAJkMNPyUPffWhcnUAZqXEokSdSsCfGXYb5PNFktCV+qs5TtJti/5nhgiV9EpRD4hXxWL8gajXUgp0TtyNr7rDlzEvORwQf+LKnQ4rzWheP1su2u6krmQiuYBQB24wIiTdeMH1egbNNP2k/aS3zoKkx0tPl7q/rDjOCU6COpbXfCRy2h/BobM4DgO/74kRdA/qXBqqfHywgp3ZS68k8ITgEepC+Pspd51Wov7g2YAwLCFFxgiKQPsSMfIwgODwxbkKMOxhZE7IN+BB36yNIUaYxLmuj4nUdIwk+zj1UyMO1sjwVVfia4UmVjY7zUt3egbGHY6vmINJvHfD3ofNhRXY2l6ND6/Zt1dZcSG4NKtLqxWxeIXz2agUteJdfurMH+ytX62eEIE8K0bUykDzuoYAaDhtT5yGcLH+qHR1If5yeG42tZDnxVWB8o7IbgH76TgBQDXOwVnLxVgbzSkQhzZl5rkE/jIrDV9nU1EYkE/IqrmI+PwowVKHKhspMZhz9os/Lm6GSVqA3xkHMb4yelqUSx0F+jvI4heyVFG4KXff4VzWhNSo4Nwy9QnCOl0Nm6uVtHu7sTIRJCdGIqqxq6HmhAIXvvoa5TYEtjiQ8dg1xmtYNf09zSYnlzbkWBeyBhftHTdR2r0OFxvv+dUlO/b3O18V+EVxPPCLa7Vmd9B6ruTdR00ZJGAxM0DoCJ3wxYeaYpgSZ6a9SGQou+7zmgxOGxdvvrKOcxVhltLUPLAsNmCP1e30CpqMo5DdmKo4KX3lcswzybcRrh+8v83SzQ4rzUhZIwPrrf3ICthvKT/gYU7nLknXHZBrhLZNuG87MRQtyYEZ/4QErZLMp4Ly3UI8JXBwgN1Bmsb2RDl3eUjInGsUXbHv/IgbWwy9Xrkb2Gd9jnKCDw9dQJauu4jcpwfGtrv0ep54mtJ+X68E8LDwat99BjDndoAzsTX3K2hTETqCDXxSU0rti2ejC+uW6kjcZEaoq1EahDvPXsTAMDD6ov49yVT6PerVbGobb2Lszc6ka+KxbsvqJCVGIr95xupcNuWQ2oUvjQTP346hQrFEYpoWkwQSq+0I10RjJau+/CRcWgy9WGsvxwDwxaHmk7u1P71pGh8UYUOH1W1YLZtp+COcJ4jDSpSe3tJWhRSosfhS70JQ2Yeq2bEIDM+BO9/oafnZ3WrooID3NK08kSnyVEbX/tBChXHc6fOAXueE3VtKDp7Ez4yoHfAjDVz4vH5NelniYz7g9Z4fpLgrafghUdw5ndwxydBKACWsiF0AflcvEIlfoKhYQvMPATOXGAkyoZw52JHryOKS8z1x4cGoLnLWs/5uVnx2FBchSEzL6Cg3IG7VBvrQCe1AHaU1uP7UyKxf91sj3wKUvQUG320/kAVZBywND0apZp2+PvK8MyMGBy+ZKDV1qTqFTiju5xRMuz9Jf2sM3Rj56kbKMidhAOVjYIsb0/8WeSZuD9ohq+co8WTSAEjksEs9Sw96mCKxxHu0kdeQbwnHITXZ8MSg8b4UMewu1o45CVlvyNx8+LPyd/TY0NQqTMhwFeGf7XJaZDcA3K+LYfUlP8nL/+StCia0UrA8tXkWlsOqZGvikWJupUKAz43Kx4FuUnYdUaLjFjpzGBHcCVqR/6dERcikFr41Yl6BPjKsNGmdFqQq8SXOhOOadokHdeOssyJfAVr4P18rAxwfGgg/H2t/16RoUD4WH+HonaupDZcSVhI9ZPkbwT4yii96O6zw143IzYE53Um/GiBkk4iRGRQ6hkTT2API2XihRVen8ITDrKK3bwgCdsWp2LzgiTsKK2HXDa68d+Vuk5cae1GgK8MvnIZ5irD8d4aFQrL9IJVthSvDzjPD2CjbMobjHhzxRRcbbNq7Gw6WIMDlY1UcsOTvjjjr9nviNQCYBXva75zH/tszl+CjblJaL5zn3L6RRU6uz4Qvt5RHgGReJgeG4JdZ7RYn5NIFWRJXL/4N67yEsg12ckjLyWC7oak+jlstuC4TV+LJNo9yLNTqevEtfYeu7Y5yucgY/A/ScrkcYCXPnrCIY4AIqv10QxZFK/4AaHUszvnd0YZOKKyfnWiHg3tPYKQTU9F7QD3qTYATqkT0oe8lEgcUbcKVE/Z2H1nkTUsDefLxPYDEMhoL0mLolQM+b24pCXbJnLtvJQIHFEb7KguNn/jeG2bXQ4JubfuPjtiwcEHvTdeOIaXPnoC8ChCDslxRI8+RxlO5ZDZc7L1flkD8iAThHh1B4Cu7qQyf6XgjAZxRGVNiwmG3thrFz3lCdiVdlHFTQSN8RHUvGYzbwHp2hVSfchXxaKwTI+e+8MC3wGrLpoRFyJQtpWK7//V8evwlXOCncl7a1Q4etlAazoTo0tKerIKpTnKCLpbXG3z8WxfMQWFZXqBTDYZgz1n9ZBxwn6ybRT390EoHW8uwrcLL330HYYnMhbOwL7kV1q7selgjd055yWHP5JrAQ9eB5oNfyRtzlcpUFRx0y0aSKw+SiQz3F2JiimRbYsnY0dpPV4uvoiiCh3l1xdMiYTZwmNgyEKNpBR1UqnrRFHFTepgJ6G5eSmRdhMkuddk/Cp1nQIjT3w0wxYeMyeG2k3cv3g2A/vWzUJhmV5AfbETGoHZAsxLDkeJupXKmBBDz47BXGU4/Hxk8GHoP7aNziAOZf3FsxnYuigZG4qr7e7No3rOvXAP3pDU7zDcDRl1FmJIQgHfW6PC89kTaV1ccb3k1Zlxbl3rUUGqzermLrzxaS38fTkadvpJTSstuelOvWV3whcdjdfes3q8tXIa7TOpif2nqhacu9EJDsAzmQrUt/egrbsf2xZPRnv3AKKCA6yTxLCFhnQSw7Z1UTI+qWnFsvQo/Lm6FfOTI/DZ1Q74+XCIDPIX1DmWyYA3Pq3F4LD1vuxZm4XVmXH0fPvO3YSMAzrvDSIzfjwN1yUhp1J9J/QhO34l6haUXmlHvioWRzVtmG5bpc9KDBOE4ZZeacPG3ElYNUMBTUs38lVxbpc4lQplffvoNTw9NQoffNkkuDfuPudeOIe7IamjulPgOG4px3HXOY7Tchz3H06O+0eO43iO41zyXV4IwVIQjhQkna20SFlFlj7YuigZYYF+dmqgjpyPowFHKqXbFk/GzlM3kJcSSflxdiUrBekdRiyKKm7ijcMawYp1d7kOJ2rbsG5/Ff28qEKHdfurEBMSYDe+BblKvPpUEoYtPHgA/UMWfFFvxNZF1rKQYv1/AkINFeQqsXlBEg5daMa85HDUNFmznasau+x2bFaV1mi7e03u5561Wdi/zlqM5uXiKmw6WIMlaVH0mlK7K5JASFRKX/r9V/jwQjNenBOPd1/ItNvlsLsX8m92l+fOjo8cJy04aJR0grvznHvxaDBqkwLHcXIA7wNYBmAagH/mOG6axHFBALYCuDBabXmc4SqaBHAeNfPKU0kC2etKXSd2ndbiTt+gnRpoUYWOGpUjaoNdYZRHCUdtLshVoiB3EqU22MnMkTEihplE+BDjs23xZBzTtAmMr1wGfN3cDY6zcvSvfaTGO6X18JFzkr4IdvwBoETditU2/4CzDFvWuJotwGqVAue0JhTkTkJBrhJ71mZhZUaMW0bTkXR2RmwIVs1QSPS9E9sWT6ZUD4n6So4cS+t6v5OfQcd1SVoULcDD9lsqE1oqw/mNwxpJZVNNS7dgkUEmealoJXeecy8eDUZzpzAbgJbneT3P84MA/gTghxLH/RzWWg39Et954QRSYX8biqtRVCGUjRa/gFOjg+xC+IjxYePryTkLy/RYlh6FHaX1yEuJRHlDJ3U+jubLKbU6fBDjQPpo3WFEYOepG8iMD0GaIoQa300Ha7BiVwV2nrqBhVMiIeOA+4NmlKgN8POR4bUfTLbbiYj5dVIV7vNrHdQ/4M6qNiMuRCDJTXZlv3g2g5YxZY3mXGU4lqRFUaPJTjDs+JAqdyN9l95d5SgjBPLaV9uEobqrZihwsq7DbqfZZOq1G3+5DNhQXC04ltCR4t+z+TEnajuQGT/iIyD37OhlA944rBE850SmXRxq+zByHV6MYDQnhVgAzczfLbbPKDiOUwGI53n+2Ci247GFVIw2oVcINUCKpMtlwIHKRgT4yqBpFRaHYY2vo6I4rd/0Y7UtEeylORNhtoAmmhE86hdTPAGQ1e6D5E3kKCNsOwwDlqZHU3oGsBq9/iEz6gx3sTQ9Chtzk2C28Bi22MK1eR67TmvtHJtk/AEI6JvsxFCq0+Rq4nIWz1+p64TGlstReqUdmxck0WutmqFwW4sJgNPdVVGFDkdswnpXWu9iWXqU3biSSYjd/ZBdiBTFJ3bmk4mX3fXsOq3FkrQobFucin3rZtlRZgBofQ32OSc7NrJ78TqeHy1GMySVk/iMJkVwHCcD8C6sdZ+dn4jjXgHwCgBMnOjVNSGQoksKcpXQG3ux6WCNNRrFbAHHcfj1ievgOA4/WZICvbHXLjOXNb5isJms5Bg2hp7NdSBG8mFDBsUx+XOV4bQuglSikqvVuLiPJJ7/5eIqWCw8hs088m3S2oF+cgyarY+qXAYMmnnIZRa7c5K+7S7XCcby6+ZubF8xBWYL8NyseEmJCIKjlw1YkhZl16ejlw04WddBcw02HazBb0810BoTUpnijhK5yLkchdH+18kGfH9KJN59IRPTFFYJjjVz4umEQ9ofOc4+S9pR5jMJcWaPLVG30s/0xl4AwnrQJPHOlWQFOXbLITUix133Sls8YozmpNACIJ75Ow4AS0wGAUgHUMZxHABEA/grx3HP8DwvyE7jeX4vgL2ANXltFNv8WGDVDAVK1K0gEhJRQf5Uk54YbhIxAsDO+EolSEkdQ5yTbIKTpqUbdYbuh54gpAzcvnWzsPesXpBHwSasOTq/o/Y/PXUCStTWR9IqtpeJogod3imth7+PDMunR6NEbYCvnAPHcTh62SBpeKQcyFI5GMS3wSbtnazrwHtrVIIxIveFPQ9pKwk3lUrQI3+TdpB2kWvkKCMQNMYHO0rrAVgXEL86UQ8ZB4EER2NnLz6qasGPFigFSYVSkhVS+SLiCXiuMhwAaLbz3gq93eRG7mWOMgKR4647leEgx4qv681neDQYTfqoCsBkjuMmcRznB+CfAPyVfMnzfDfP8xE8zyfyPJ8I4CsAdhOCFw8Gwm8Pmy1oNPVhStQ4nNd20kgjQh2IE6SAESMGjIRiSq1ASQlFUvu3sEyP6+09VDaDGIgH2do7ymV45akkt2LWWYcnS/PsLtfRBK1STTt8ZBz85Bw+v2blzM0WYOGUSMyzaSXlKMPhI+MwNykMCeFj6TUd0WRS7WYN93trVDimacO6A1VYt79KQD+RPrBlSQHY6B0rHVVvk+YgEwwbBMD+7UgqpCBXie9PicRvTjZg56nr0Bt74SOX0ZKplbpOHK/twKoZMdQncvSygZY8JbTUpoM1NHLLFcVHKMw9a7PwSm4S+ocsGDLb77zI2LrjM5I6rsnUaxextelgDZpMvZLn8EIaozYp8Dw/DGALgJMArgH4mOf5Oo7j3uY47pnRuu6TDjY8cVl6DIYtgIwD6jvuUcPNvmiuEqS2HFLTTGYWZGVKXszyhk5aHnOebUfC8sePKnyVdYyvKfrKYX1mNqSVFdfLiAuhEVb+vjL8YcNsFL9sLYm56WANMuJCsDE3CV83W43ploXJMFuAr/R36G+J8XXXf8K2JUdprRUwOGwBONC602xyGNvHl37/Fd6xKaq++4JK4CcgQQB5KZG2IIAIgYwFuY/ie7cxNwkyDlQvaeuiZOworcf19h5BhBIxuB13pWNAOu72200AO0/dsAtxXpkRg5UZMQBGdhu+cplkRJM7PiNHx3Xc7cew2UJ9F6SCmxeewat99JiBlW3eUFyNrITxOK81YZoiGG3d/S5rE5MXTkqqmV3xEomFbYsnw2wZEdablxyO81oT5iWH45zWRHn6R835Eu0dcX1mqXrMeSmROFHbTjV+dpfr0GTqFWjzVOo6cfSyge4GWBqCUEpRQX7oG7Jg66JkAT3miVYT0TlarVLgeG071QzatjjVjuYildXmJ4fjjxvnCs5FZEdIpbvZtgI+OcpwHCqY65RKIQEIZLXuK5fh6alRtIqbWKJ8Q3E1np2pwPHaDoE+1nmtyW7B4IiuEfdN/DfgvpyFo+OOXjbgmKYNA0NmDJp5+Mk5+Pt6Jo/+OMNbjvMJBnnhiDNY/H9XBlos+CZ1vp2nblBDkRkfgqrGLmxdZF1VN9/pxYcXmhExzg+me4Nu1Q3wpCwooQWmx4bgSqv1e1KfWdxHtmTluy9kPvCYvvaR2q4UKAA7w+b8HCNteW5WnKDu9E+WpAgmsmkxQbZJVWik2bEh0VNPT43CEXUrfOQcAnzldkJ4rOggifpZmRFDHcdyGRDo54P1OYkoqriJZ2cqEB82lhpisgBYmh4tOXG4g2+L7ycTOMGbj6Ds6eMCbznOJxiERzZbIIhJJ387o3KkuFpCZ7BUxdL0KByvtTowo4IDAFh17zPiQvDXy23wkVnlFkgil6uwUXfLgrLGMCE8EFsXJWPIbKH8d5oihIZPvvaRmoZaljcYXbaBcOriDOg3Dmvw+bXbiA8dg2GbXAVL+0gZRnESV6WuE3+70oY0RTA+v9ZBOfb967IxPzkcO0rrUVRh9XeQCWF+crhkVjELq2x1GwJ8ZQjwlWPromQBvXRM00azm99bo6JRP0mRY/HHC7cwPzkcZguQHDmW6jkdutAsSEwsyFVi2+LJOKJuxezEUBxRtwooInfwoHpXnoDQgr5ya+Cjr5zDrtNab6Kbh/BqHz2GIKUkq5u6EBUcQHVvSNnMjrsDkto07Jb++eyJVG+GaN/0Dgzjgy+bMDsxFKfrjVSf5umpUciMH48th9TQG3vRZOqDr1yGHy1Q4qimDZsXJOGC/g60t+851MRxVRa00dSLnZ814JapDx09A9izNgvDFgt+efw6OABzk8Jx8moHjmnasDQ9GpdvfYPqpm+Qr1Lg3RdUmG6jTJrv9OHpqVGSbSB6PKTcpUwGvHWkDlcNd8Hbvl8xPQZ1hru4cPOO09KPrLZPc1efLVeEw6//0Zop3NBxDwCQGh2ELQsnY6y/HDtP3cBX+k46IVxt66FjL9YUIrpDLXfu42ZnH360QInNC5Ro7x7A+nmJVIuof9CM8zoTZBwQFRSAErUBWxYq6W7K31eOhPAx+OzqbTSZ+nDgfCMWpEYiZvxI6dE3SzTYe/YmnslU4HS9EfkqBT6paXVLa+rbxH+XaXHVcBd+PnK8+lQSag13MWy2QMZxDu/5k4T/EdpHXvx94am6pLOCJWQHka+KRVVjl11iFgkRrNSZwHEQKGfuOq3FMU2bywgkqQxmglUzFBgyW3BeZ8L6nEQAwK7TWsg5YMjMY4JttzJstuDdz27A1DsIXzmHz6/ddrpSZFf0pL+7TmsRExKAHaX1mBgWCI7jIJdx2LM2C8/NikOgnxw+Mg4HKhsdnpt1Fr9/ZqS+AclUXpkRg4s379BomYJcJZQTrDITieGB+OPGuYIdgnhVTf7NFqUhn5Nj2V0fAJrNrjf20vv8ap4SKzIUmKYIRom6FXMmheFMvZHuFIoqdPjwQjOemhxBnc/lDZ0ugwek5C4qdZ1Yd+Ci5OePKunRxxbqum1xKvaszYLPaGqxPKbwjth3HI5ePhJ66UjzSAqOtvhkciHaO9tXTKGGgaV1/njhFuYpwwU0BqnOtTIjxuMEM3G/SJjt3go9LcBS/PJsmmk9NToIHMdh2MzjDxtm4wNbVBGhT6RkstmJc3e5DnWGbgyZLagz3EV2YijqDHfhZyteA1h5+d//yyy8viwVSZFjnWZU5ygjMDU6iE5krE8EAG73DMBs4bHpYA1e+v0F1LbehZwDTL2DAtpOyvi6itRhv2clOC7d6sIxTZvgPJsO1uDWnT5sXZiMq209iAkJwDul9Xh+dyV2lNbjxTnxqNSZEB82RiB94mySd7QgeZQS7GIkhI+1S+zbszaLBg944R689NF3HFISxFsOqRExzg+BfnJK++w6o8XS9GiYLbxLWWMxiFwyoSaIjDb5++hlA353Rov31qjw46dTkBk/Hm8fvYa8lEh88GUTXn0qCW+tTKPnk5KmLqrQ4Y1Pa7FnbZYddUUkoAtfmomooABKh/x06RQAQGGZDtNjg3GxsQspUUHY+cIMSjsR+mRWQih+/HSKXd9Y2mqcvxx7z96En48MKzMUOFN/m/LTVY1dAIDXfpCCo5cNOPjVLfzyHzKwfHoMNC3dGDRbUHqlTTC2lboROWtNa7dAzvq1H6Rg1QwFSq+0oXdgGE13+iDjgD9unEPlI2QyoL17QJJ3J/eEXJulmAbNFvyspBbr5ydiQlAAHbvkCeNg6h1EV98QDl9qhd7Yi12nbwCw7mLiwgJhtvC43NKNITOP1m/6MTsxFDdu38PtnkFszJ2ErIQwKuN9RG1A+91+SZlxTUs31s9LtKMDH1aC3ZkM/Kt5Sjs6i9CmXnjpoycGjnYDrDomoX2OqFvtlE2d7TQIWEqCXYWRvxPCx9rRTpsXJOFEbbvkql9qFckWjGH7pWnpFiSfsXHuv6/Q0x2M1tiLfJUCVw13UWcYSbwjx9c7qcdMaCsifsfzvMB5+9oPJgMAjmna8JXOZA17HLagztAt2EmxuQtsvohYzpqMFclZIBJLFh70nCTCKyMuRPIeEdlzkqxGdhZEkO6f58SjsEyPo5cNdOwKy/R4fekU7FmbheQJ41CibsXEsEAqL7HpYA1WzVAgXzWiBnuxsQs6Yy8C/eRIU4zUjy4ss4bEOqMoHdGBzmhCV8+mlOCeV/fo0cK7U3gMIFU8hazmSKQQkUsuLNMLHISfXmrBzs8akBk/nu40Nh2sgal3AAnhgbQQD/mO/M2uioljm6BSZy2Ywq76NxRXw9+XoyuOerP+AAAgAElEQVRN4vg9d6MTf/iySVAwhu3XrMQw6iBnneCZ8ePx/hc6PDMjBv/vYguNsiIOW1KIR8ppLl5NvnFYg4+rW5CdEApT7yASI8airbsfsyeF4eer09HePYCNuZNwy9SHTy614tWnkvD9KdYorM+uduDDC7eo83b9vEQ0d/UJCvI42rEUVejozsRXbvVblF03osnUh09qWikV4mg3+OzMWBoi/PbRa6hv70HR2Zv46bJUFOQqMT0uBL87o0VYoJ9gRd7c1YePq1uwND0aFTc6YeF5nLBlUMs44JCtnsK1th46YT0/Kw67TmudBgGwRZkAq+O3RG1AQe4k/PHCLUH7f3niut3nYjgqxEOKKnkL7ngGd3cK3knhMYCjl6yqsQuRQX70xdyycLJdFEugn1xQae3npdcAAGu/l4C3j14TROKwf5NqXlJgq3MBVuPu78th56kbmJlgnXyau/pQom7Fzc4+p1E8zs45M2E8yq4b7aqhzUwYjyNqg+BzMhGJq4JV6jqx87MGAMDbq9ORFDkWhy+1wk/OwdQ7aFs5x6G5qw9/sFUE++OFW1iXkwgLz+PsjU6YLRZcuvUNCl+aCcC6cs+MH4/U6CCBEVyaHo1zNzoxM8FKI/30z1fgb0u++2FmLI7XtsNs4XG17a5k5bENxdXQGe+hsExHqZgLN034qKoF02KCcLr+NvJVsViQOgGlV9qQr4qzWyywvoaC3CQ0mXpRojZgxfRoTFME46PqZuQow/HFdSMGhi2YEjUOnb2D6Bs04x9nxklWrAv0k6NE3YrzOms9iLiwQGw6WIOGjnsofGmmYFImz5E7k7WjiLTVmfb98sI13J0UvMlr33E4yxQF4FJxkpzj5eIq9A9ZaIawpqWbKp+SLFxrZbBvaGawu2CzrEl7iir0sPDApqeS/q4ql2zbSDIZB2C1KpZScOKEODZXYmp0EC42dsFPzuHVPCXtC+mrowRCUhGNzaouqtDhv042ICshFOpb39DMbIKXfn8B57SdNBGPtOP+4DCGLcDsxFBaQ4FNYJPKTmeflbyUCJyo7cC2xZPx21MN6B+yZjq/OCcex2s7sCw9Ch9eaEagnxwb50+yu1+kHWyG9MqMGEHfyHEkE9uTJDZHyZSunmsvhHA3eW00VVK9+BbgSi7ZmfqpMxCumOgZpUaPo7IVnr6A5FzvrVFRLhkYyTb1tG2O8CBZs+znGbEhOK8zUeMDWBPAxIKABNmJofi6uRv5tugnKVnpDcXVWJoeJZhU0hQhdm0iPP2B9dlURoRVMy2q0OG8thPzkyNwRN0KgEd5QyeemRGDDy80w0dmdWZzsIZllmoMOHzJQCdw8RiLFxPPzbIadg5WCknGcfjr5TY6uZSoDZgcNQ7bFqcKZMzTFCHUd/Ln6maqOiueEIARFVQxxJ+z97FS14miipuYnxxOZb/J5EqSMR/Fs+PFCLz00XccYj4fsG67tbfv2dEtjugTklj16lNJ0LR245imDatmKJAQHoiiszeRGh2E6+33MD85HNVN33ictMT6EKoa7wAA/H1k2JibRBPr3C34TiAVhaJu7sIbn9ZSiooYPimqS/x7Eik0KyEUpVfaBf6P1apYNHf14a2/1GJCsD80Ld34njIMxzTWwjc9/cPQ3bYmo9Ua7lL/THxYIAaHzZKF6MX9FNNjWQlhGOsvx69PXMfFmyYcutCM7SumID4sEN/0DeLsjU5kJYzHmXoj8lIi0GS6j0EzDx85h3yVAh9VteAnS1Oon0Y8xuLrNXf14ZimDc9kKvD91An4Un8HMg74YWYsNC3dKHgqCT9dMoWei9CBMs4akfX7Cj1O1HYgX6VAk6kPgJVWEkdkuQPiSyBU0zMzYvC3K+14PjsO+881UT8Sab+nz86TCi995IVbeOOwBsc0bdSpSSYJsgqeFhOMc7YV6tW2u27rJ4nBUlRbFyY73B24u9qXos0cCbdJCQCKaTa2DGmdoRs7Suvx/SmRtM6AWDcoIXws5DJrAh35HQDBDo387kFpjt3lOpy7YaQZzn/cOBcbii/idL0R6Ypg1Lf34IXsOPz1chuGzRYMmXnw4CGXyVC8fkQk0B19ISmK70BlI1ZmxOAXz2ZI3hdWE+mIuhUBvjLss0VakfF0RCO5ahO5nyRIgjx3hOrylML0wqt99ETDnTBTAkcJPxbeWm6zpqkL+apYOiG4o58khaOXDfCVywTZt1LncTcLWyoUd9viyThe24G8lEhbXeNIh0lWzjKOC3KV2L5iCs5rTYIwUhKSeUzThr6BYRSW6bEyI4b+jmQrE/rOUXKZu9m+chlw3pbhfF5rwspdFThdb8SiKZEwdPdj1QwFDl1oxrDZgn3rsrFqRgzMFoAs9DwJ12TlxUmb96zNorWZpe5LYZmeiuStVsVi37psmrA4MGxBdmKooGQn0ZFyp03W+2Atn/rSnIkoyFXSsOGC3EneCWEU4aWPvuN4EBqF/Q2hn0i4Kfl7QrA/DSs1W3jMnhRKQy5JmKS71EClrhO/O6O1i0RZPj0G+ao4QXtYqumjqmb8uboFhS/NFDg12XayUShbFk6GTAYUnb2J2YmhOFN/G9sWT7YLdSUgvydhpmwUS1ZCGIaGLQLdoF+euI7Cl2ZSddGC3El4a2Wax/Rdhijihg0xffvoNZqw99aROsg44O79YSjGB0Br7EW6IhgfbcpBfXsPStStSFMEIy81EqnRQTYjHYWrhh7Ut/fg4+oWj3YnUhFe0235EFkJoVg+PYZGAv289BqyEsZTDaijmjYsnx6DsEA/7DqjxfeSwnGl9S6WT4/B8ukx2FBchU9rWnGzs1dwP509M2xEnUwG7D/f6DKM1QvH8IakPiF4kFhu8ptGkzUpieQAkBj70itt6Lg7IJgA3j56DZsXJKG9e4D+3llYKgtHxobwwOI+kHBV071BymtLTW5ShoMYRneE25zFy7Pf1dzqEoRb/vLEdcycOB4nr3ZQ/wH5DZmwHPl6yOeOQi3J5zeNvbh15z7+fUkKhs08ag13ER86Bre67uPMtQ6UXTciXxWLWsNdLE6LomGeBblKNJn6cPZGJ5amR6PARn+5A0dtnhDsTydxYvR5nkdL1338r6eTaajzpoM1qG68g9ToIFxq6sKzMxW23AYzrrTehYWH3eQLWBc2n15qQaCfXJArk5UwHokR4zB7Uih2lNZj2+LJ9FqOwli9cAyvT+EJAvELkJoChM9+/4zWLpqGgA0jJPVyAcf1AUY7DJA9/4HKRgDWGgnsv9nrPkqfgqNQXtbfQEJVfZixYv0QzsbOEcShluLP81WxOFHbhv4hC+YlR+C8thMcZ818fnFOPN7Jz6D9dlRYaLVK8UCcvtT9YcNOLTzwXFYslU8HgA3FVTDzQPH6bOqXSYkeh+vt9+xCdsXtYcfy6GWDrT6EdaxJeDR7H0kfAHgccfakwl2fgnen8BjAUfJQR88AXrXlAYhXVWyWrYwDmu/cx75zN92mah51whB7fhkH/P5fZtHMZbZv5LqOEuT2n2sSZFI7SrRztnshuyQAVDfoh5mxqG/vQVffEDWymfHjcUzThiaTNbHN3Qlhd7kO6uYuAR0ik1n1lQbNFrpDOaxuxcCwBf4+1rj/S7e6YLYA6YpgZE8Kp/dlZsJ4wQ6OJKbNTBiPnZ814JimTZCxTnZbbLY6AXvPWZDdW/+QBT9aoMSWhcnYdVqLzPgQ7Dmrxy1TH4z3BrE6U4FhiwXWuYOH+la3NcRVJsOWhcmUgmKfx/iwQDqWJepW1Lf30IVKjjICsxLDkJUQJmgT2XU5yvZ2dxf7JMGrffSYQ+ys9JXL4Ocjw3+X6bDxA+tOikgISxVpEcsqV+pMgmLqYielKwXThwWrsuorEmjylcswTxkuuK6UoqvZAkFUijOVUWdFX8h3bA5IjjICf/nX+XTlSo5fn5OISp3JqYaPGKR06eYFSdi2OBWbFyRhR2k9mu/0Chy9KzNiMMZPjqXp0TbKBpifHA6dsVdwX8iqWCpnZc/aLKzMiJFUyvVEWt1RoEBUcIBA0jwpcix2lNaj6qaJOsl5265iy6GR4j/ie0LGsn/Igv4hi0BV1hkcaX95HdEPDi999B0FeYHZTFmSjeoj4/C6Tf+GPZ4YD6lwzCFbMRK5jHOLqnmUL5+j80tlEj/KJDdxdi9bi9od6uFBKbXd5TqaLe6q5nFRhQ6/OdkAGQdKX5Ea0Wx2tCuKyBFV5U4fXFFt4mxmklBHCgWxz6ij+uBSWdGe1FZ21D8vRuANSX3MQVZIO0/dQF5KhK0MoXUlN8ZPTssusseTl5FVHWWVPH+YqRCUtiQvpLPiO2J4Eg5LsPesXlDeMUdpVQn9fxea3bquu9ckx7GqppsO1uDl4ov0b3dDOFnD6Gg35qhtGXEhMFsgUAotyFWieP1sO8NeWKbHP8yMxf512di3LhsDwxbsPNWAvJRI7Dx1A++tUSEpcqxT5VBnu7wcpWvFUvH917R0Y/OCJBp2u2dtFpalR6N/yIKBITMsPJCvisU5rQl5KZE0jNnVhAAA+9dlU1VZUoDInXsxmrvYJw1en8J3GGzGLOHh7/QN4XvKMPypqkXAs0qpmrK8enxYIAL95PjblXZkJ4TizHUj/b2zSBoxpDjeDcXVmJccjqwEYSY1yRDOSggVhGL+d5kWf6pqwX89P8POSEld1xmvzPLm5LiE8ECEj/XDwa9uYWDIDJ2xF0vTo/BJTavL1T4JnyX1AthdBqlvwLbPUdvMFgtKr7QLfD4k8ostubl+XiLWz5tEw3UDfGU429CJq213MWeSNeu5sExPo83q23uw81SDIBmRVZcVR+64o1gqvv8kGi15wji89oMU1Bm6UXT2JuYnR6DR1AcfGYfLLd1IihyLihudyEuNwOrMOMxKDJP0WZReaUNUsD/+bXEqfRYz48cDAAaGLU7Dnl31z4sReKOPngCQF2JqdBA0rd2CKBiyVWd1h5wZu0dJEYkpCUc0EPs5ALtoKE8zpqVoEHE/iKbQalUsjtsie2YnhuJiY5db1IOn4yRFFS1Lj8Kfa1rh7yMdyeTqPpFxMlt4DJt5bLdpSL320dcoUbdSwTxyfUfROeJnw5N7zkY5HVG30ja8WaKh1NHllm4Mmy2CiK1Hzfk/iN7Vkwp36SPvpPAdhTshlZ5w3Y/65XKlbEl4cLatheU6DJl55KsUePcFlcftcJc3JyJ/Ab4yLEuPwRF1K1arFChv6PTIILozviPGMwIlagPyVQqqSErE5MSSEq6u+94aFb7SmbDrjBa+cg4BvnI8PXUCjqgNWK2KRXmD8Vu552zo7LsvZGJ3uQ5Npl403+nDOa01FDpojA9+c6IB2TYFV68T+O8Hr0/hMYcznt8dnphg3YGLKKrQCaJxiip0WHfgosAP4QlYjreo4iaKKnSCNuWlWK+zobgadYZu+vmQmUe6IhglagOKKoQVzFzx/O7y5nkpEfj8WgcCfGXgAHx+rUOy5rQzeDK+xD9yRG3A7MRQHFEbaE4Be571OYlOJwRA6AsifQ3wlSN8rLVmxmqVAu++kOnQvyGGswgsV2DHu7zBSH01xzRtuNzSja0Lk3GgshG7TmuxIiMa5z2M0PIED+LH8sIxvJPCdxTOXmhPHG/zksOxo7SeGmFCr8xLDn+gdokdsNsWT8aO0nq8WaKxlQVV4IjagKTIsfS7/y6zXtvfR4amO314cU48dpTW47WPvvaY+nJUxJ7ILx9RG5CdGIr967LxPWU4BoYtSFOEUEeoO7pOnowvcRavVsXiYmMXVqus1dIqdZ0eO0il9Im2LkpGo6kP6YpglDdYz0kmor1n9U7P96BwNN6kDCqLYbMFn1+7PapOYE9Ca71wDa+j+TGDp443ItG8o7QeldpOfFTVQvlhKTgrnC4lyZyVEAZjTz8+vNCMpenRgrKg2Ymh+EpvwpCFR2r0OPQNmgEA2xanwsIDJerWB67Kxhax33JIja2LknH4kgHPZ8fhnNaEhPBAHNO04389nYz27gHkq+KoQ/VROjZLr7Rh9qRQfFLTSjWCNi9IwgX9HfzujFZwHrZkqdTYivtaqbNKmjw7U4FK3R3qbCZyH2+tnDYqzla2DbvLdYgKDsDy6TE4ojbg/zyThv4hCz74sgkrpseg0dSHlRkxeGtl2qg5gR3JhnhpKiG8yWtPCMRbZxIuyCZYOQsffeOwBmmKEGTbnK0p0eOgN/Y63Hq7WpVJ7WDiw8bSQjQk/PK9NSoc07TB31eO2YmhuN5+D09PnUBlDsobjB6toB3tmgjlUpCrxL51s6iKKgnnJM54Fs6oB0e03d6zekkKo8nUS53pZFVdWGZdwYvPs23xZOw8dcPpipftK2nLO/kZ2LduFpWWJn0Tj8mjolnYNpDnAQCK189GnaEbR2zO7hO17di6KJnSYs6eRSl40l5PKD0vnMO7U/iOQxzySMIF2TR/Z+GjOz9rwKc1LWjuuo9Um05NQ3sPNuUpJVdzzlZljnYRF2/ewen624KwRwD408VmbFmoxKmrt7E0PQpH1AYE+MpwwlaPwJ2VuKudCxtOSaQ0PviyCT9aoMTz2RMdhoxGjPOjAm3seTvuWncV4jEhonHi85CwTfEuZmDYYnceUl/a3RWvq75J3W92N0LaOHtSKKoaux6oSA37PNS330XR2ZvYvmIKti+fhpkJ42mosatnUQqeSFi4E1r7pMOrkvqE4GG2zvFhgWjt6sPXzd3wkQE9/cOQccCQmceUmCABjSH+nZQOkvglfuOwBr8+UQ+tsZfKZl+4acLvzmghl3FYnBZFM13DxvojLzUCe8pv4sc/mOywYpgYD2s4cpQRkuOXGh1k15ednzVgY+4kSVVUR/eB5BiIx0+qL4SKIUqkRMfKHYlyZ0aRTJw5ygj4+3LYUVqP6sY72F2up9nRD6MVxE5I+SoFti+fRj9/mKpohHpj1X6JUq+4eqA3V8E1vPTRE4SH2Tq3ftOP+ckRGLYA/UMWvJo3UmDGEYiDNEcZjgOVjXSLTxycG4qrsfPUdRzTtMHMA1sXJVMOvKqxCzLOeh6zBTRXISMuBAW5Sjw3K84uG9tV393RvnHmjJYaP/a8a4q+wl++NgjOV1Shw4biagG187AURoZNfvpAZSON3tl0sMatyCtnjnaW8ivIVWJecjjOaU2YGB74QFX0pK4/EonUaRf5RRzkntJXGXEhlBIjkWNSRZM8ybj3wjVGdVLgOG4px3HXOY7Tchz3HxLfb+M47irHcRqO405zHJcwmu15XOEoikX8Eu4u16GoQid4Ceclh+Ni4x0E+MoQ4CvDgcpGpClCULx+tsNrEQO0xSamR+QISKTN0vQoGmZZvD4bhWV6arD3rM3CvnXZOFnXQauXsS/0qhkKWu2LvZ4zw/ggUg2s4XA0fuS8lToTLDyPrYuSseWQGq999DXV9xfH+P895BZcGUV2gnvtIzXOa022uts9yEuJcDghuGPEXU1ILDyNEpIK52XlUAgeJrTWC3uM2qTAcZwcwPsAlgGYBuCfOY6bJjpMDWAWz/MZAD4B8OvRas/jCmcvpfglJOqcRISUhJ/KOatWzb8tTgEgNPKsAdhdrsPRywaBcujWRckYGLLg/TNamqVc3tApUNOUWoU7MuLurvzFY+DKGDsyHGxWr3j82PP6ymXYdVqL5MixtPykWHCQPc+StCg77R5Hq2KSK6JpsWalr89JxK4zWoQF+glUWR3BHaOYo4ygCXTzksNh7Bmg4cEkHFkMd4y4J6t0T++ts3Beb27C6GE0dwqzAWh5ntfzPD8I4E8AfsgewPP8FzzP99n+/AqAdN1ELxzCVRIb+xIWlumxfcUUunLfeeoGvj8l0lbfV0E1dFZmxFCxM9YAZMSF4GRdB/2bvLQrMmJwXmeiNZFZA8vSIcRguzLiziYNsTEgsg9L0qJcrlQ9GT9xjeU9a7MwMGzBxcYuzE4MpQlbjs5DsrWPXjbQdjpaFZNckZHJ2hqdZPjmPuoMwuziBzV8lbpOHL3chnRFMM5rTdi8IAnvvqDCmjnx+NXx65Lj9SAT9O5ynV2+AttmcSKheKfF9o1E0pFItPIGI42s8+YmjB5Gc1KIBdDM/N1i+8wRNgA4LvUFx3GvcBxXzXFctdFofIRN/O5DapVIXhpA+BJOiQ6iBdCJI3P/utkCDr2wTI/Icf44aYsAcrWKZ1/aE7Xtgu390csGmC08VmbEUIO9obgK6w9UOaUb3jiskZxIiCwDezwxuqtmKKhBYVeqrowoGT92sslRRiAhfCwy40Po+esM3RgctiAxPBBXWu/aZT+L70OOMoIWvndlUAtyrX6cHaX1WLvvAvqHLHhzxRT8dFmqILHwQQ0f+d3ry1KhM/ZizZx4FJbpUVShw/HaDry+LNWh8qzYiIuPc7UbFbeZLAjEuxSpvhGfgjicl0hzeDpheeEefEbx3JzEZ5JCSxzHvQRgFoA8qe95nt8LYC9g1T56VA18XMFSIgBwoLIRAb4yXGntRlGFTrBKn6sMF6ySiQHYujBZ8gVjj8lXKQQ+gbnKcGw5pEaaglQw64dcxtFVc52hG4NmXrBCFO9sKnWdOKZpAwDMVYZjrjLcTiyOGIOX5kzEyboOgYgc6bdU3Qh3x4zQSu9/oQUAJEWOxc5TN/DinHgct+kWkagpTUu3oC4DQaXOqiHkajwJCnKV+KCyCS1d9zE7MVRATe08dQM994cfuAwqu4shekskV4MUJWInATIemw7WIDsxFF83d1Mjvn1FsODc4vvxxwu36G5U3Gb2fuQoIzBNEYwdpfW4auiR1GtytQt253n1wnOM5qTQAiCe+TsOgEF8EMdxTwN4E0Aez/MDo9ieJwbk5WFVR/evs9bNfae0Hi/Oice2xanUiLOKqo4mDAKW+imquClwtrIvbUZcCPae1VPnbF5KhM2o2GdLkxccAOXVgRGhPABYmRHjcvIi199QXI2l6VECgTtipB05H6WMG2mH9XzRtB4xMa57z+rxylNJtK1k8jl62YCTdR3YvCDJ5XgSFFXo0GqbEKoau1BUoUNBrhIFuUr03B9+KMMn9i1IjZ14Ujx62YCBITO+qDditS0RbeEU60SSpgixG1NyznnKcKQpQgTXACDYxZFrFuQqcdXQgxJ1K20LK9InrsdMdmPk92Ia0tn4euE+RpM+qgIwmeO4SRzH+QH4JwB/ZQ/gOE4FYA+AZ3ievz2KbXnikKOMQEZsCC1tqGnpht7YixfnxOPwJQPd7mfGh2DnqRu0wMyStCjMVYbbOVzJSpLl2UkWrVQIIktHWZ2zVsE2R/IZBOyL70gszplPIkcZgaXpUShRG+iOxF3aRezLIFRJQe4kmo0NWA1cjjICrzxlpZCOXrZGxVgL9lThmKYNy9KjaGaxK18HcfhvXzEFH7+aQ6mkogrdI49ochZpxdIxxzRtAMdhXnI4StStyEoYjzP1Rjw7U0EjtsiYsufUtHZjQ3GVZFitmGKr1HXaZa676ytgn8VAfx87Os/rdH5wjNpOgef5YY7jtgA4CUAOYD/P83Ucx70NoJrn+b8C+A2AcQD+zHEcANzief6Z0WrTk4RKXScu3foG+SoFrWlAVq+rVQrBLmLfull0q370sgGbDtZgZUYMra51TNOGlRkxaDL1YklalMPtvBjWiBerTLXVOdspoChctV9qFSimIMhuh90RlDd0UroD4OyoCUeS0WSFT65JDD0AauB2n9XjJ0tSBP3fUFyFITMPmYzD4LAFy9JjcPiSQbCLYuVH2JWupqUbX+lNgh0U+f8xTZskPbckLQoAaAlOUrNBb+xFQvhYSp2Jd0auxk68iwga44MdpfWYnRiK81oT1tjos5fmTBTsjNhzkN8EuHF/HbVFvGOToszE1BLZ8ZKJ3F3a0At7jGqeAs/zf+N5PoXneSXP8+/YPvvftgkBPM8/zfN8FM/zmbb/vBPCIwB54bYtnkwloQvL9FiWHoUdpfXoH7JgyCwskE5WcYT/P6JuxW9PNeCIuhWAlVc/WddBvycQhz6yKKrQUR0crbEXmfEhboVpkvZL7VqOXjY4nJhYQ/PuCyqstukt5aVECoyK1Gp008EaHNO0CVb2u05rMWyrGUwg54Bdp7WCPvAAhi08Bocttnj6ViqPzV6TTbxiV8DF62fb7aAKcpVYlh4jyakD1gmDjCVx7h5Rt9KV+4biaursJXA2dqRNZCImsterVQpbOKjCphsVIYgKExtnswXYvmIKVBND6S5PKqz2YaXf2V1HjjICS9KisOu0Fn0Dw4IJwbtb8BzeIjuPGUihE7KKJAZvYlgg9MZeLE2PpgVmXslNklyJVeo68XJxFfqHrAYx34MCNOw5NhRXU+NI2nF/0IyE8ED8fHU6gJEVHrvKJSt58r2Yq3e06md/R47NS4nEidp26lBl28euRsWrb8AaBdVxtx/9QxZU6qxFY+Yqw+m5SXGcYbMFPAALDwwOWxyOF7kmWylPvGtwN+GKjCcpdM8B8JHLsD4n0Wm1O0djR74n/qUmUy+OqFsF51yWHoXDlwxUSsPR8yAeW9bPJOWMF/dZ6veunjv2mSX3yRuRJIS7RXa82kePGQbNFvzujBbLp8cgPiwQzV19KFG3oq27H6tmxODza7ch4wBfuQybFyixfHqMnU7M/9/emUdHVeX7/rOrEoLImAQzEQMkBIEwhMEBRbxAI6K00qtvD7S00g5oX65te187dOvq9+jb9ur7VtNv2SoNqOjjPtu2VVCICjYaQYIKIRDABEiFQCZiEkIAw5Sq/f44tQ+nTmo4CZkk+7NWViqpU+f8ap+qPfz27/f9lTc08e7uKpp9ErcLvqo+5UjC2krO3moevHk4d00YYmr6zBufTP0359h19ATv7q7i3d1V/GJWBs9tLuFgzWlTV0hpCVn1hGL7GIVk7F9yq/bR/OwhlDc0sXhNPgdrTrP87ok8MG24KTJnFyZut4wAACAASURBVGaz6jc9c8cY+vRyBxzXp5ebFVtK+frUObOW8tyxSSy6cZj5Wpc/8e/2cUnkHzEkPI7UN7FkRnqARo9V1+itXZW4BGRc1ZcdZQ3me2iN/lBqbB/OnveyzVNPs0/y81vSyU4daL6fJTNGRNTEsrZdYUUj1w4bZOogjUzsx/t7j5my1y4XvPhJKY/fNjLg3HZ9IevgY9Uh8vp8rNhSyoTUgaZG1eI1+ST0j2HWqASzjV7MLeG5zSX8YlYGS2aMwOWCx/+xl4LyBqobz4bUUCpvaOL9vcdwCcg/2sD7e4+x/O6JekCwoAXxehBWpVDVkS5ek89nh+p4+bPDuF2C28cmss6fzfq/7sxi3vhklrxewNyxScwdm2SKlqkvq/GaJIqqTxHtFuyrOml+ocNdX6GK208eGmt2PnPHJrFkxgh2ljVQVm/kLH5x+DhulwhZmziU+J71eXvnN3dsolkE3nqMVZgtzxNZHO93OUWAEQpr7eBcLnhlWxnZVw+k7vR55mQlsjy3lOV3T+TOCUYqzhs7KgI6+fNeH4vX5LO7/AST0wZR2XCGTw/UMbhfDMtzPUGF3sKR56njdzlFuAREuQW7jjZQWNloDl7q/Thtu1GJ/XhlW5k5cOTsreb+acNYdOMwAHaUNTB9ZDxeHwECgHaxu1C1LSpPnOFgzWk2FFZz9rzXbNv/mD2S1Ng+5r7IG19WIATsPNJAZUMTq7aUIQQcrT/D4unDzWOtn7c8Tx2LVu/g+uGx3JZlJFJaB922iPFdjuhBoYegvkxWieIP91eTe6CWI8ebcAl46Z7JeH0wfWQ8GwqPMXdsktkBFlY0mgVmAF7MNWbtj8zM4K18wzeef+SEUa1sW1lQ5clISqWpsX0oq/+GZR8dJKewivwjJxiZ2I/a0+do9kluH5vEtpI6ruofE/ClL29o4sXcEtYWVDEnK5E3d1YwMa3lwJQa24dPD9by1q5Kc9av7LCqmFoHhFCqmtaOdEraIJbelWUW65maHo/LZeQNPDIzg6tjryShfwwvf1bGL2ZlcNeEIeYKw+uTnGv2mdcsb2gy8y/GDRnI/qqTuAQUVZ9iTlYib+VXOpawVgM3GPc246q+fFJci0vAv83IMFd/XxyuJ2fvMR66eTirth4mJlpwrtlntoly34xK7Ndi4LDKcqu/J6XFBtgWTO3V/jp13KxRCUxIHcjagkqz037pnosuPSX5/otZGew80sCZ8172VDQiBPSOdge4/+yftxdzSzhw7BQVDWfMgbGg/ASfHao3BxKNVkntMagNTBWS98u/F/D7nGKEgDHJ/Yn27zY+ND2dMckDuHVMQoBQmt2fmxZ3JSsWTjJLUz4wLZ0VCydx7bC4S9K0mTc+mbMXvOyvOkWUS3C49rT53NqCSlIG9g7QbFq8Jp/7X9vJBn+lsk8P1vLY7BFBwzrzPHXsrWw0Bf1UGG2oMFSn4nhFx06Zbayu6/VhJq+NGzKAeeOT6RXlMpVd1XXnjU8OaFurrtHagkqykvvjkzC4by/WFlSSndpyI1ptktqlPVZuKSV98JVm7obXBwuuSyUzsZ+5UfvwLcP5zK90e316HI/NHsHvc4q5/7Wd5mb0ktcLcLtol5DXS9UisoYwj0rsR7PP2Ots9kkzGMLalipoYdmmA2zcX8P3JqbQ7JOc8++DuV2CmGjdvbUF3WrfcoLlA/SKcvH4nJFUN541k8dWbfUE7azsqKgOe3SH9f+h7IgUMSKEwO3yR+p4JTdlxNOnl5sol+C9PdXmwPa5x+jMmr2SWaOuMjdMVcU268CkOrcVCyfxyr1TAPjZqztYvCY/5CajU3E8JWoHmIPegWMnAyqbTU13Jmeh2k1JPOyvOklWcn9qT58ndVBvNhfXcltWQtC8Cnu01I0Zcewpb2T44CvN5z/YV8MTc64xr+P1wauLprBi4SSWvF7AqTPN9Onlptkr+dxTb24q22UkrNdpTUcfKb9ArW6i3S4emZGBlHDfqztanD9pQG++LGvA7ddD6BXlCpBnV9faUdbAuWafKb/xj/xKolyC28clho160kRGu48uA1Jj+1B87BSbi79mwBVRSCm5f9pw5o5NYun6Ivr3juL9vccClutW14pTwlU5O+/1hSzyojqIlT+dRPWJs5Q3nDE2ZI838fD0dP595ggA+l/Ry3RlTE4bRPbVA1lbUGW6NYJFq1h92NbN1ylpg3h0Vmar2tHuD+/Ty226fBbdOIziYydZW1DFvPFJASGkkfY9rG3w/AKjBGha3BV8XFzLyMS+lNY1kZXcnzzPcc43e1tsCtv3TV74xMMPpgzhb19WhCw+Y92sV661h6eng4C3dlUyJyuBgzWneeaO0eZ7n589JGCfwGkBI2udZut+zNyxieaehHJLqv2Z3tEuPimupaHpPHdOSOGpdwr5z5wiak6eI8ol8EqYec1gyuqa8PokH+47Zu5ppcb2IdotyD1QS5RL8FX1KWKiXDw2O9Osha2CAuzV7Xoy2n3Ug7DmA/ikMSNXM9zpmfGU1TeZhW2g7cJq9tngU+8U+jelL4aOXp8ex61jEgKOU+4agKJjp5ifnYyUhntLyWv/4XvjAmQ2dpefYG1BlZl8p1Y6dputs3676yecKyTYLHicv0NUWFcBv/z7btb57bEXknGSdWx1WamZ/YLrUjlaf4b52SnsrzrJpLSBIVda1pWYktyIVHxG2aZcay99dpiCow1Eu4UZdAAEtKt1NejELajabcnrRp1mZeMFry8gp0WtWq0SFwuuS2XLwTqWbTrAuoJKvD6J2yW4OTOen1yXyubiWtwuwbQR8dwxLikgn2J5binzs1Muupm8Pp7bXOIog1wTHr1S+JaT56njqbf38fhtI/n13NFMSB1Izt5qvD7J3spGth4ysnuPHj9jRn60pmSnFfuM9c2dFQAM7hfDL79jzMqXvF7AL7+TGRDRNHloLOUNTabLQm1gf3LAkEJeur4Il3+z/PkF2WYJyiiXoKy+iRnXXMWqLYbOkirTGawdQm0eB9todDoLVquwtf5B988/zA44t3pfka5r3YBVZSb/9mUFKxZO4oFpw6k9dZacvceYn53C+sLqoKGeaiW23i+j8Y+dRqb4x8W1QdtGvafld08k46q+fFxci8sliHIJotyCPE99xNDNYKsg+4oxNbYPLhc88dZedh1tIMotiHa7mDc+2Tzmqv4xAcEQeZ46nttcwuwxiby2/QjXDYvlp1PT2HX0BEfqmyiqPkWvKBfzs1P4y48nMmtUgrkKuth+5Wb0lU/CDelxLJkxIuCz2tZSoJcjOvqoh2DNB1B/35AeS83Jc+yvOsn87BT+dbKhS7i/6iTbPPWtzjmwYu0kHrp5OA/fks5f/AVhrIONPTJFuWaONZ5j0Y1DuWuC4apQf68rqOKZO0YHhEN+d0IyR+ub2PRVDfOzU4i9MibkF9zq+rG6M6yhtlZ3mdPa1nmeOpZtOsi88cl8erA2IHS1sKKRd3dX8b2JKQE1pV0uWFdQxV3ZwZXiJw+NZUdZg2mv6iAfmz2C2CtjiO/bi2UfHWwR0z93bGJAzsCNGfFsLz3O/OwU3sqvbDGQWNtkR1kDg/vFsL/qJNcNiyVpwBUcrmsKcLMFcykGC91N6G8EBZTVf0OfXm7KG5p4Zt1+mr0+zjVLrh0Wy+/uygoYHIO1t5ogPDBtGB8fqOXeqUMZ3DfGzL14eHo6z9wxpkX7nff6eGbdfsCIYLpzgiHYV3XibEDYdKha2D0Vp4NCR6qkajoB+8avitwxMmtT+GdRDf8squGRmUb1sIlXD7okRclgmkRO5IutLgmF2qi1/9+qjll07JR5rX+dPCSkbpH1HNYN43AS2natn2ADwpLXC8xwSHtWsHVz2qocqjZvw2GPTLp1TELAOTYUVvPS1lIKKxo5Um9ENil3jIo2emdXVYBOk12DynoNU7XWL2HR7PXRO9pFYWWj2X72NoqkT7R4TT5rCyqJdruYMnQQn5cex+UyVqhwsa6F9R6Hk11fvCbftAsMyfdgn9PCikbuGJcUkH2+YuEk1u+pCqnDpXGOXilcZqh4eLdLMD1zMLvLT+D1SXaUNbBi4SQenZUZ0bUSimAuGpWQZU2autS48L9+6qGgvMF0J/1gytW4XPDU2/u4MSOuhRvC7vZpzSog1OY4hE7EsroknF4rHJOHxgZkU09NjyfaLViz/SiD+8WwufjrABfPeb//3JpUt3R9UciMaOt9GxLbhw2F1bhcguuHx3K03sheV24k9b4nD40N+/7nZw+5mFHt9VFx4izRbhcv3TOZeeOTzT2tJ28bFbS939xZwSMzM8wVVnlDE+sKKvFJWL1oCndOSOHt/Eo2FFZT0dBEn15u837XnDzHzZmDW9wHq4tJ0xK90dxDsdf5XTR1KHdOSG5RjyBUzkGkcwfr8KyV1dpjc2/cEEPOW1VxU7Pvx2aPMPMnIm1+RgqRtXaUoWx3UvtYhWdar2X9v1Psm7qqNrGS7LbH6YfKswi2gW4VwlOfjxULJ5HQvzcSOHvBx9iUi1pTatM53Pu3rhjB0HyaNeqqkIOhvb3tsuuFFY3clZ3C6kVTzFXYr+Zk0uyV1Jw8GxBWrQIbdOnNjkEL4l2GtEVQrC2EcuW0RtgtFJHew7JNB0y3z2OzR7b69U5sd3KMNbt40dShrM4rA2gh2+G0rdT7UlFOrb2HdpdPMCE867FWUb1otyuk3EioawAsXpPP2Qtemr2Su7JTglZRa81nJZhQ3+ikfnxWUs/8EOfXRMapIJ5eKVxmhJoBP/VO4SVlnAbDyUy6rYSb6UcKAW2vVUB7Fod3ci71vlQtiIdvGd7qFZjTMFJr0t+D04abcupOsIYYq3O89rNrGZ3cP6hUObTus2Jtq6np8UzPjOezknpGJvYNunLStC96T+EyI5Qf+GDNKf7ycUnEEEwnhEti23mkIeRzrfH3hvL3Owk9dbIX4AQn+wUqUmpw35iAiCz7tSKdy57cNn1kPMs2HSImWnDXhCGU1X/DwZpTbC6qYXNRDbNGJQQkDtojqyIl06k2Asx2LvRvECvV0lCo8FprO5c3NPHmzgrmZCXy4b5jQTWq2tLuxcdOsa6gipsy4th9tDFkyK4mMjoktYehOur52UMCMolVpuqsUQmXvCGqCBfjr8IVL2XwCdfxF1Y0RuzwQ4mytWUTMlIHq3IwrANYqEzacDH/1vd13uvjWOM5bsyIY9mmQ0xMG8jIxH4s++ggX1WdpKT2G6LdgqXri0xJaiU7DkYy4wufePj5LekhN/+tuSOqnSekDuQvH5eYoaaRBnbVzoEDmgwIBsjZW01BeUOA0F+ep45n3t1nCiAGu4Y1P+SmjHi+qj7FY7NH8FZ+pZnbYn9f4SYregNabzT3OJy4J5zoEzkhnIvCqfsiHOE2Uh+anm4K11ntGTdkQIdU2XLiqrrv1Z0tXD2rtnoC7Pnrpx5WbfUEnGvVVqMgkqqSNjU9PiBL3OuDl++dbOg4fVwCGIV0Zo0yKuglDYhhu18nSvGbtYX8PqeYx2aPiFi7OFg73zomwdzYVa9ZtdXDfa/uDOk2s2drK4HGwopGsyqcqgKnPpc3ZsRF1Er6cJ+RzJd/pIGHbxlual+pYAN7oER7uvt6MnqlcJmgsjyXri8KqYcTKQTTTriZ1/zsISFn0E7cF+HOX3PyXIuZtnWmf97r475XdxITLZiUFmt++Z1KTzvFqavqxow4lueWmuGkSl77wZsvyjYXlDfwrL+zVsVjns0p5rsTkrkhPY5nc4o5Uv8NH+6vodnrY0dZA/dPM3SD1hZUcriuiYduHk526kBe236Ea4cOYnd5I1FuwS+/M4Kl64soPnaSN3dW8pPrUvnVraM47/Xx1Nv7+O74JM5e8JkrPNVOD01Pb3H/+/Ry88q2MnM2/t6eKt7eVckTt400w0fts2/rykyt3J56ex+D+kSb2evLc0sp9osJrlg4yUxeDLZyte53qCJJamUwNT2+hRS69TPSXqvhyxHtPuphKD366ZnxvLb9CHOyEngrv9J027RWBkKdM5grKL5vLw7XneaVbWXmAONyYXbITgcfp1ITdlJj+xATLfwdaZNZpEZVDWsvX7OTvQlVZ8DaGb2yraxFFI8qUrM8t9Q8RoXYLr45nSP137C2oIoJQwZwvOkCAIP7xvC7nCLcLsFDNw9ndV4Zu8tPcPvYJD4u/pr52ckcqW/ii8PHGZnQl83FtczPTua/vj8hoJ1e/KQ0oJjPsk2HuDEjjklpgZnLVlG8peuLyBh8JbvLG+kV5TIHOCf3KDW2D57aU/46GAn8eu7ooGKCoSYPz7y7L2iW+HP/LGFEQt+w7iGnE5KeiB4Uehjqi7Nqy+Ggejj2Dk6tLOwrCSdSEM0+X9BZ7/SR8Zxr9jkefC5lZjcpLZYj9Uap0bEp/dn01deXNCuMVD3OarN9hmottak6I6XfZB08dpQ10DfGzWvbj5glM897fbyYW8Lm4lrGpvTny7IGbh+bxPTMwWa5z5fumWwmnTV7fZTVN/E/bjUUQZfMSGfLwVqONpzh2qGDKChvDGjrYO30yMwMc2UTrKM3/Pkn2Vxcy7VDB9HQdKFVulkqr2ROVgLrCqrYWXacj7762gyztV7XOnlQexmT0gaZKwNVaOmNHRUsumloxMTF1q6GexJ6UOhh5HnqWLq+iDlZiWwu/rqFHo7KUFUdn1pZKFdCqDrBwWZe4Wa9NSfPtSryp60zO6Pj8Zgd6ZyshAA569bS1lWLeq01s3t1XhkbCqsDNn/BcCGt2nKY+dnJrC+sxuWCZ9bt52DNaZbMSGfTV18zJyuRdQWVFFac4LphsdSdPs+88ckUVjRy/7RhuIQgM6Efj88ZxdghA3hvdxWe2m8Yk9yfo8fPtNiEDdZOv547OuxgvGqrx7SzoLyRJTPS+aykjm2eeiZePZBrh8WGnK1bV6QPTEtnZ9lxPiup56aMOFb+dEpASVNrxvrYIQNY9tFBNhRWM298MnPHJpkyGp7ab1h+98QWLqffvvdVQFa0XSOqrZn7lyta+6iHUVjRaLpQQunhHKn/hhc+KTFdGw/fMpxnc4oZndyflVtKg84Ag2kdqdjyU2eazQSycB2yVePITrDzO02eUu9XxfWPTu7f5oHBukHeEUl/avb869uvYXluKdMzB/NsTjHDB1/J9cNjA/SScgqrOO+V/Js/W9i+kW9lu+c4v779Gh6Ylh7QLmoTNlw7BdN9MgQADwWcc/GafFxCMDU9jsLKRhavyTc/Q/ZENuumc56njvwjJ7gpI478IyfMvIPnF2S3+LxNTY839YvUPbjg9XH2go8Hpw0PoZ+UwvLcUlMzav2eKuCiRpQ1QEHvKzhHDwqXCVYROLt4mWLe+GQ2FFazeE2+mX0b7RbsrzoZVhAu1Dlb25nbCXf+cOeyDoDq2NHJ/Vm26ZDZQbQFa4cTTtwvmD0rFk7ic0+9+drr0+MCOiNrZ6kG0/nZKfSOdvHOrioemz3CjD6KiXbzyMwMM9oqVMdWWNEYULvYHqX11089IdsJWt4/Jb738r2TzQiv/VWNnG/28b2JKcwbn8z6PVVsKKxm0eodzB2bFJBdbBUmdCImGKz9p6bHM7ivkdXdO9rV4vNln0So6Kq7r7uajftrWuzlhJuQaIKjZS4uE5zKCOR56vjZqzs4e8FHtFvQO9rNoqlDWyUFsX5PFRv31ziSUmgPm9v7taGwSmOs2mrUb7CuPMKd36m0SLDjgA5boQRrp1VbPSzbdKhFh20dQMCQrzjX7ONXt2YyJnlAwITg9zlFpjT7n384ocVnINL9CffZ2lBYzQWvz5TdUO1jH+DUNadnxrO2oCqk5InGwKnMhV4p9GAueCWPzzFcP8Fm6cE6v6np8easV7korDNUwHHHHOr8TjrES3ltMOydWr8rong2pxggwDUTTBI70opHdYAQKE9de/qc+XdbVih2gnW0qpqc9X8q/yHYCsPqQgOIiXJx6kxzCz2l6sazzM9OYV1BJSD59GAdt465mAlt1YZSnwfr/bGvbK0aUkoW29peoVxOKprqUlasGhtSym/Vz6RJk6SmJdtKamX20k1yW0lt0L/V/7J++6HM+u2HcsHK7TL9qRw56pkPAl7z5Nt75PLckna7ZndleW5JgJ3Lc0vkyi0lAe995ZYSec3TH8g/bSwO+77s55LSaAt1LtUuT769R24rqQ1oJ9Xm2Us3RbxOJNrjfizPLZGPvrFLpj2xQf5pY7H808ZimfbEBrlg5fag51THPvrGrlZfXz2v3rdqH/sxoT6P3+bPX1cA7JQO+tgu7+Rb+6MHhdDYv2T2L8eTb++RWb/90Pz/yi0lMu2JDXLR6i8CXm/vLJ10eJfaoTkhki2twWmHojrFP20sviTbQrVTe3dsl3o/Fq3+Qg71d/JZv/1QXvP0+3LEb96XI379vvl+rLZe8/QH8iertstrnv7AHOTUpMPJ9Z22bzDa8/PQE3A6KGiZi8uISDIWaXFXBmzEPTAtnd/cfg3bPcdDSlKoZb5SWbVKByjXQLBrBtP1v1RV1vaUMXAixxFJ4qI1toW6N6EkPVZuKW1T+12KlEmep44dZQ30jnaxcX8N5y54OXPBh1vAr+Zkmu/H6s9/bPYIU5fovld3klNYxQWvj22e+oi1JZ56p5DVeWUB7duaz0hHqvT2ZPSgcBkRqRML9iV6YFo6D0wbFrITUZ3UhsJqfvbqDkOHx7I56nYR9JodoUPTHrpK9vOFk+eOJL/dGttC3ZtQHduDN7fULHLSfq0ZyOyoKKqX753C+WYf572SXm7BXdkppu6Qdd9I5SI8vyCb5bmlTEobxP/7ohwBZtnPxWvyg9qc5zFKjgJcnx7H8wuM8p6hjtd0Hh06KAgh5gghDgghSoQQTwZ5PkYI8Xf/818IIYZ2pD2XM63txKyvC9aJWGf6U9PjWTR1KGcv+Dhz3svnnnozGkRVSLOLwVk3LdujA1dcykzYvnrJ89SxOq+MqelxLTrQcKJ8oc5rtW16Zrx5bFvuTVsGwLZ+BhTBBiiXS5ibvqGOVe/7s5I6ot2CKHfkbsVaAW7J6wV87hf2s1YI1HQNHTYoCCHcwAvAbcBo4MdCiNG2w+4DGqSUGcCfgT92lD2XG/YOTsXuh4oICka4TsQ601edZ+9oFy6XMDtkrw9T7Ex1imqgUG6G9lBltdvc2pmwaivre1q11cOi1TsAWDIjo0UH2paiMKu2eli19TA3ZcSxrqDKVAa1lsNU57Hfm2DuNoBRif0ct19rBrJQqCggI1zZhcAITVWlMEPN+v/7i6NMTY+jd7SbWaMSzFKwKxZOCnp9azSS+owsmjqUP3xvnGNbNR1DR64UrgVKpJSlUsrzwBvAnbZj7gRe8z9+C5gphBAdaNNlg909oySL7VLZ4fyr4ToR9Xjxmnx+9qrRef7H7Exiolz0jnaxOq+McUMGmO4DNaNdnlsaEP/eVldGMNo6E1ZtBZjv6Y8fHMAlLpbNbEsHqlCD4bM5xUxKG8S2knoWXJfK8txSVm31sHF/TYvZtv3eBHO3LV6TT2Flo+P2aw8fu8oKfvneKbxy7xSi3C7OXvDyvzceDFvB7fkF2bz+wPU8MjODdQWVzM9O4b+/OGraFYr2/oxo2gEnu9Ft+QG+D7xk+Xsh8LztmH3AEMvfHiA+3Hl19NFFOiPyZ8HK7f6Qw4KgYZTqmvYoko4IF7yUaBNrW418+v02R7yEs+3RNwrM8MzspZvko28UmFE5TrDaqEKHOzvc0t7G6r7+2B+SGu54ZaMK7XUakqpDSjsHHEYfdWTyWrAZvz192skxCCEeBB4EuPpqLYWraKssg1PyPHUUHTvFIzMyzAxf66pi3vg6c2YdTDIh3CqkLVxKwpq1rYLJJ1wq44YMYOWWi7pT0zMHs7agslX3xWrj1PQ4llhe2x7t5wR79rt9Fm+/tvV4+z0HwtrcEZ8RTTvgZORoyw9wA7DR8vdTwFO2YzYCN/gfRwF1+KU3Qv3olcJFOnKl4HQW922Z7an4+ZFPv2/OwNvLVvt5Vm4p8cf6F7Tq/J2Z8+HUlu5+XzXOoauT1/ydfCkwDOgF7AHG2I75N+Cv/sc/At6MdF49KBh09JfWqavm25BAFCmj+FJtvRQXit3G7tIJfxvuq6Z1OB0UOlQQTwgxF/g/gBt4RUr5eyHEUr9x7wkhegNrgGzgOPAjKWVpuHNqQTyDjhCEu1zpzLZq67X0/dR0NE4F8bRKqkaj0fQAnA4KOqNZo9FoNCZ6UNBoNBqNiR4UNBqNRmOiBwWNRqPRmOhBQaPRaDQm37roIyFELXCkky8bj5FY113pzvZ1Z9uge9vXnW2D7m1fd7YNusa+NCnl4EgHfesGha5ACLHTSShXV9Gd7evOtkH3tq872wbd277ubBt0b/u0+0ij0Wg0JnpQ0Gg0Go2JHhScsbKrDYhAd7avO9sG3du+7mwbdG/7urNt0I3t03sKGo1GozHRKwWNRqPRmOhBIQhCiFghxEdCiEP+34OCHDNBCLFdCLFfCFEohPhhB9s0RwhxQAhRIoR4MsjzMUKIv/uf/0IIMbQj7WmDfY8JIb7yt9VmIURad7HNctz3hRBSCNGpUSFO7BNC/MDffvuFEK93J/uEEFcLIT4RQhT47+/cTrTtFSHE10KIfSGeF0KI5/y2FwohJnYj237it6lQCJEnhBjfWbaFxYm+dk/7Af4LeNL/+Engj0GOyQRG+B8nA9XAwA6yx41RqnQ4F2tTjLYd83MCa1P8vRPby4l9/wL08T9+uLPsc2Kb/7h+wBbgc2ByN2u7EUABMMj/91XdzL6VwMP+x6OBsk6072ZgIrAvxPNzgQ8wqjxeD3zRjWybarmnt3WmbeF+9EohOHcCr/kfvwbcZT9ASnlQSnnI/7gK+BqImBjSRq4FSqSUpVLK88AbfhtD2fwWMFMI/aVgpAAABBVJREFUEazcaZfYJ6X8RErZ5P/zc2BId7HNz+8wJgNnO8kuhRP7HgBekFI2AEgpv+5m9kmgv//xAKCqs4yTUm7BqMUSijuB/ysNPgcGCiGSuoNtUso8dU/p3O9EWPSgEJwEKWU1gP/3VeEOFkJcizGL8nSQPSlAueXvCv//gh4jpWwGGoG4DrLHjhP7rNyHMXvrDCLaJoTIBlKllBs6ySYrTtouE8gUQmwTQnwuhJjTadY5s+9/AncLISqA94F/7xzTHNHaz2ZX0ZnfibBEdbUBXYUQ4p9AYpCnftPK8yRhVI+7R0rpaw/bgl0myP/sYWNOjukoHF9bCHE3MBmY3qEWWS4Z5H+mbUIIF/Bn4N5OsseOk7aLwnAh3YIxm9wqhMiSUp7oYNvAmX0/Bl6VUv5JCHEDsMZvX0d9H1pDV34vHCGE+BeMQeGmrrYFevCgIKWcFeo5IUSNECJJSlnt7/SDLteFEP2BHOBp/9K0o6gAUi1/D6HlEl0dUyGEiMJYxodbVrcnTuxDCDELY9CdLqU8101s6wdkAbl+b1si8J4Q4rtSys4o8ef03n4upbwAHBZCHMAYJHZ0E/vuA+YASCm3+8vsxhPie9PJOPpsdhVCiHHAS8BtUsr6rrYHtPsoFO8B9/gf3wO8az9ACNELWIvhr/xHB9uzAxghhBjmv+6P/DZasdr8feBj6d/B6gQi2ud30awAvtvJPvGwtkkpG6WU8VLKoVLKoRi+3c4aECLa52cdxkY9Qoh4DHdS2FrmnWzfUWCm375RQG+gtpPsi8R7wE/9UUjXA43KNdzVCCGuBt4BFkopD3a1PSZdvdPdHX8wfPGbgUP+37H+/08GXvI/vhu4AOy2/EzoQJvmAgcx9i1+4//fUowODIwv4j+AEuBLYHgnt1kk+/4J1Fja6r3uYpvt2Fw6MfrIYdsJYBnwFbAX+FE3s280sA0jMmk3MLsTbfsbRuTfBYxVwX3AQ8BDlrZ7wW/73s68tw5sewlosHwndnbmfQ31ozOaNRqNRmOi3UcajUajMdGDgkaj0WhM9KCg0Wg0GhM9KGg0Go3GRA8KGo1GozHRg4JGo9FoTPSgoNFoNBoTPShoNJeIEGKKXxO/txDiSn/Ng6yutkujaQs6eU2jaQeEEP+JkVV+BVAhpfxDF5uk0bQJPShoNO2AXxdoB0Y9hqlSSm8Xm6TRtAntPtJo2odYoC+G6mrvLrZFo2kzeqWg0bQDQoj3MKqSDQOSpJRLutgkjaZN9Nh6ChpNeyGE+CnQLKV8XQjhBvKEEDOklB93tW0aTWvRKwWNRqPRmOg9BY1Go9GY6EFBo9FoNCZ6UNBoNBqNiR4UNBqNRmOiBwWNRqPRmOhBQaPRaDQmelDQaDQajYkeFDQajUZj8v8BBxSqIc09DSgAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N = 100000\n",
"darts = np.random.uniform(size=(N, 2))\n",
"plt.axis('equal')\n",
"plt.xlabel('x')\n",
"plt.ylabel('y')\n",
"plt.plot(darts[::100,0], darts[::100,1], 'x')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 2: Calculate the squared distance from the point to the origin.\n",
"\n",
"1. Tool: Slicing to get all x and all y coordinates. \n",
"2. Tool: `**2` to square entries of an array."
]
},
{
"cell_type": "code",
"execution_count": 107,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.50799735 0.62220291 0.86391143 ... 0.23784649 1.12973904 0.22278458]\n"
]
}
],
"source": [
"dist = np.sqrt(darts[:, 0]**2 + darts[:, 1]**2)\n",
"print(dist)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 3: How many darts landed inside the circle?\n",
"\n",
"1. Tool: Boolean selection."
]
},
{
"cell_type": "code",
"execution_count": 111,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ True True True ... True False True]\n"
]
}
],
"source": [
"inside_outside = dist < 1.0\n",
"print(inside_outside)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 4: Count the number of `True` values.\n",
"\n",
"1. Tool: Use `np.lookfor` for `count`"
]
},
{
"cell_type": "code",
"execution_count": 109,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Search results for 'count'\n",
"--------------------------\n",
"numpy.bincount\n",
" Count number of occurrences of each value in array of non-negative ints.\n",
"numpy.busday_count\n",
" Counts the number of valid days between `begindates` and\n",
"numpy.nanmax\n",
" Return the maximum of an array or maximum along an axis, ignoring any\n",
"numpy.nanmin\n",
" Return minimum of an array or minimum along an axis, ignoring any NaNs.\n",
"numpy.count_nonzero\n",
" Counts the number of non-zero values in the array ``a``.\n",
"numpy.ma.count\n",
" Count the non-masked elements of the array along the given axis.\n",
"numpy.fromfile\n",
" Construct an array from data in a text or binary file.\n",
"numpy.fromiter\n",
" Create a new 1-dimensional array from an iterable object.\n",
"numpy.nancumsum\n",
" Return the cumulative sum of array elements over a given axis treating Not a\n",
"numpy.bytes0.count\n",
" Return the number of non-overlapping occurrences of subsection sub in\n",
"numpy.frombuffer\n",
" Interpret a buffer as a 1-dimensional array.\n",
"numpy.fromstring\n",
" A new 1-D array initialized from text data in a string.\n",
"numpy.nancumprod\n",
" Return the cumulative product of array elements over a given axis treating Not a\n",
"numpy.ma.count_masked\n",
" Count the number of masked elements along the given axis.\n",
"numpy.bytes0.replace\n",
" Return a copy with all occurrences of substring old replaced by new.\n",
"numpy.ma.frombuffer\n",
" Interpret a buffer as a 1-dimensional array.\n",
"numpy.all\n",
" Test whether all array elements along a given axis evaluate to True.\n",
"numpy.any\n",
" Test whether any array element along a given axis evaluates to True.\n",
"numpy.npv\n",
" Returns the NPV (Net Present Value) of a cash flow series.\n",
"numpy.ptp\n",
" Range of values (maximum - minimum) along an axis.\n",
"numpy.sin\n",
" Trigonometric sine, element-wise.\n",
"numpy.sum\n",
" Sum of array elements over a given axis.\n",
"numpy.flip\n",
" Reverse the order of elements in an array along the given axis.\n",
"numpy.ipmt\n",
" Compute the interest portion of a payment.\n",
"numpy.prod\n",
" Return the product of array elements over a given axis.\n",
"numpy.size\n",
" Return the number of elements along a given axis.\n",
"numpy.angle\n",
" Return the angle of the complex argument.\n",
"numpy.ravel\n",
" Return a contiguous flattened array.\n",
"numpy.fliplr\n",
" Flip array in the left/right direction.\n",
"numpy.flipud\n",
" Flip array in the up/down direction.\n",
"numpy.geterr\n",
" Get the current way of handling floating-point errors.\n",
"numpy.select\n",
" Return an array drawn from elements in choicelist, depending on conditions.\n",
"numpy.seterr\n",
" Set how floating-point errors are handled.\n",
"numpy.unique\n",
" Find the unique elements of an array.\n",
"numpy.average\n",
" Compute the weighted average along the specified axis.\n",
"numpy.nonzero\n",
" Return the indices of the elements that are non-zero.\n",
"numpy.reshape\n",
" Gives a new shape to an array without changing its data.\n",
"numpy.convolve\n",
" Returns the discrete, linear convolution of two one-dimensional sequences.\n",
"numpy.digitize\n",
" Return the indices of the bins to which each value in input array belongs.\n",
"numpy.errstate\n",
" errstate(**kwargs)\n",
"numpy.gradient\n",
" Return the gradient of an N-dimensional array.\n",
"numpy.chararray\n",
" chararray(shape, itemsize=1, unicode=False, buffer=None, offset=0,\n",
"numpy.histogram\n",
" Compute the histogram of a set of data.\n",
"numpy.is_busday\n",
" Calculates which of the given dates are valid days, and which are not.\n",
"numpy.seterrcall\n",
" Set the floating-point error callback function or log object.\n",
"numpy.einsum_path\n",
" Evaluates the lowest cost contraction order for an einsum expression by\n",
"numpy.histogram2d\n",
" Compute the bi-dimensional histogram of two data samples.\n",
"numpy.histogramdd\n",
" Compute the multidimensional histogram of some data.\n",
"numpy.result_type\n",
" result_type(*arrays_and_dtypes)\n",
"numpy.ma.dot\n",
" Return the dot product of two arrays.\n",
"numpy.ma.sin\n",
" Trigonometric sine, element-wise.\n",
"numpy.ma.diag\n",
" Extract a diagonal or construct a diagonal array.\n",
"numpy.busday_offset\n",
" First adjusts the date to fall on a valid day according to\n",
"numpy.chararray.count\n",
" Returns an array with the number of non-overlapping occurrences of\n",
"numpy.ma.size\n",
" Return the number of elements along a given axis.\n",
"numpy.datetime_data\n",
" Get information about the step size of a date or time type.\n",
"numpy.ma.angle\n",
" Return the angle of the complex argument.\n",
"numpy.busdaycalendar\n",
" A business day calendar object that efficiently stores information\n",
"numpy.ma.ravel\n",
" Returns a 1D version of self, as a view.\n",
"numpy.ma.average\n",
" Return the weighted average of array over the given axis.\n",
"numpy.ma.ediff1d\n",
" Compute the differences between consecutive elements of an array.\n",
"numpy.ma.nonzero\n",
" nonzero(self)\n",
"numpy.matrix.ravel\n",
" Return a flattened matrix.\n",
"numpy.histogram_bin_edges\n",
" Function to calculate only the edges of the bins used by the `histogram` function.\n",
"numpy.chararray.resize\n",
" Change shape and size of array in-place.\n",
"numpy.linalg.matrix_rank\n",
" Return matrix rank of array using SVD method"
]
}
],
"source": [
"np.lookfor('count')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"2. Found tool: np.count_nonzero"
]
},
{
"cell_type": "code",
"execution_count": 112,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"78599\n"
]
}
],
"source": [
"num_inside = np.count_nonzero(inside_outside)\n",
"print(num_inside)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Step 5. Calculate $\\pi$!"
]
},
{
"cell_type": "code",
"execution_count": 113,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.14396\n",
"3.141592653589793\n"
]
}
],
"source": [
"estimate_pi = 4.0*num_inside/N\n",
"print(estimate_pi)\n",
"print(np.pi)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Summary"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment