{ "cells": [ { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "# Numerical operations with Numpy\n", "\n", "1. **Elementwise operations**\n", "2. **Basic reductions**\n", "3. **Broadcasting**\n", "4. **Array shape manipulation**\n", "5. **Sorting data**\n", "6. **Summary**\n" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "## Elementwise Operations\n", "\n", "### Basic Stuff\n", "\n", "\n", "#### With Scalars\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([2, 3, 4, 5])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([1,2,3,4])\n", "a+1" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ 2, 4, 8, 16])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2**a" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### All arithmetics element-wise operation " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([-1., 0., 1., 2.])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = np.ones(4) + 1\n", "a - b" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([2., 4., 6., 8.])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a*b" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ 2, 3, 6, 13, 28, 59, 122, 249, 504, 1015])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "j = np.arange(10)\n", "2**(j + 1) - j" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### ... and lets see how fast is Numpy compared to Python's in-built operations " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.09 µs ± 154 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "a = np.arange(10000)\n", "%timeit a + 1" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "600 µs ± 55.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "p = range(10000)\n", "%timeit [i + 1 for i in p]" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Clearly, Numpy is winning hands down! \n", "\n", "Another thing about multiplication * VS matrix multiplication" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[1., 1., 1., 1.],\n", " [1., 1., 1., 1.],\n", " [1., 1., 1., 1.],\n", " [1., 1., 1., 1.]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = np.ones((4,4))\n", "c" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[1., 1., 1., 1.],\n", " [1., 1., 1., 1.],\n", " [1., 1., 1., 1.],\n", " [1., 1., 1., 1.]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c*c" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[4., 4., 4., 4.],\n", " [4., 4., 4., 4.],\n", " [4., 4., 4., 4.],\n", " [4., 4., 4., 4.]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dot(c, c)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[4., 4., 4., 4.],\n", " [4., 4., 4., 4.],\n", " [4., 4., 4., 4.],\n", " [4., 4., 4., 4.]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Or you can simply do\n", "c.dot(c)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "## More operations\n", "\n", "### Comparisions" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "a = np.array([1,2,3,4])\n", "b = np.array([5,2,6,4])\n", "c = np.array([1,2,3,4])" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([False, True, False, True])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a == b" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### array-wise comparisons " ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.array_equal(a, b)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.array_equal(a, c)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Logical operations " ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, True, True, False])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([1,1,0,0], dtype=bool)\n", "b = np.array([1,0,1,0], dtype=bool)\n", "np.logical_or(a,b)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, False, False, False])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.logical_and(a,b)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([False, False, True, True])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.logical_not(a,b)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, True, True, True])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.logical_xor(a,b)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Trancendentaal functions: " ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.arange(5)\n", "np.sin(a)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":1: RuntimeWarning: divide by zero encountered in log\n", " np.log(a)\n" ] }, { "data": { "text/plain": [ "array([ -inf, 0. , 0.69314718, 1.09861229, 1.38629436])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.log(a)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ 1. , 2.71828183, 7.3890561 , 20.08553692, 54.59815003])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(a)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Shape mismatches: \n", "\n", "Throws a broacasting error. We'll get to that soon enough..." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "ename": "ValueError", "evalue": "operands could not be broadcast together with shapes (5,) (2,) ", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0ma\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (5,) (2,) " ] } ], "source": [ "a = np.arange(5)\n", "a + np.array([2,4])" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Transposition. what is it?" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[0., 1., 1., 1.],\n", " [0., 0., 1., 1.],\n", " [0., 0., 0., 1.],\n", " [0., 0., 0., 0.]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.triu(np.ones((4,4)), 1)\n", "a" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0., 0.],\n", " [1., 0., 0., 0.],\n", " [1., 1., 0., 0.],\n", " [1., 1., 1., 0.]])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.T" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "##### See, what happened above? \n", "\n", "It flipped upside-down and then right-to-left" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "hideCode": false, "hidePrompt": false }, "source": [ "### Food for thought\n", "\n", "- Look at np.allclose versus np.isclose. What is the difference between the two according to you?\n", "- What is the difference between np.trui and np.tril?\n", "- Also play with Ranges: How is np.linspace different than np.logspace?\n" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "## Basic reductions\n", "\n", "### Calculating sums" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "26" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([5,6,7,8])\n", "np.sum(x)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Adding by rows and by columns " ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[2, 2],\n", " [6, 6]])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([[2,2,], [6,6]])\n", "x" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([8, 8])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.sum(axis=0) #column-wise addition, first dimension" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "(8, 8)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# more complex way to doing it, but you get the idea\n", "x[:, 0].sum(), x[:, 1].sum()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ 4, 12])" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Le's do row-wise\n", "x.sum(axis=1)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "(4, 12)" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x[0, :].sum(), x[1, :].sum()" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[[0.34592414, 0.15460061],\n", " [0.8683121 , 0.68157837]],\n", "\n", " [[0.23800664, 0.93651685],\n", " [0.9379463 , 0.69130543]]])" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# In higher dimensions\n", "x = np.random.rand(2,2,2)\n", "x" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "1.5498904649649252" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.sum(axis=2)[0, 1]" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "1.5498904649649252" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# inspecting it more in detail\n", "x[0, 1, :].sum()" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "### Other forms of reductions\n", "\n", "#### Extrema" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([2,5,6])\n", "x.min()" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.max()" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.argmin() # Index of mimimum" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.argmax() # index of maximum" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Logic Functions for Truth Value Testing\n", "\n", "**Numpy.all**\n", "\n", "`all(a[, axis, out, keepdims])`\tTest whether all array elements along a given axis evaluate to True.\n", "\n", "Try help(numpy.all) for more details." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all([True, True, False])" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all([[True, True], [False, True]])" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([False, True])" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all([[True, True], [False, True]], axis=0) #axis=1 should give you the opposite result as you'll see" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all([-5, 4, 7])" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all([1.0, np.nan])\n", "# We'll get to np.nan - Not a number i.e; later" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "(140591000095120, 140591000095120, array([ True]))" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = np.array([False])\n", "d = np.all([-1, 4, 5], out=c, keepdims=True)\n", "id(d), id(d), c" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Try the above operation with keepdims\n", "\n", "What did you encounter? And why?\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "**Numpy.any**\n", "\n", "`any(a[, axis, out, keepdims])`\tTest whether any array element along a given axis evaluates to True." ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any([True, True, False])" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any([[True, False], [True, True]])" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, True])" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any([[True, False], [True, True]], axis=0)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any([-1, 0, 7])" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any(np.nan)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "(140591000888144, 140591000888144, array([ True]))" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = np.array([False])\n", "d = np.any([-1, 4, 5], out=c, keepdims=True)\n", "id(d), id(d), c" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d is c" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "(140591000888144, 140591000888144)" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "id(d), id(c)" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True])" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "# Can we try array comparisons?\n", "\n", "a = np.zeros((200,200))" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.any(a != 0)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.all( a == a)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([3,4,5,6])\n", "b = np.array([4,5,6,7])\n", "c = np.array([6,7,3,2])\n", "((a <= b) & (b <= c)).all() # np.all(((a <= b) & (b <= c)))" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, True, True, True])" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a <= b" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([ True, True, False, False])" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b <= c" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### If we would have done `np.any()`...\n", "\n", "What would then the answer be then?" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "((a <= b) & (b <= c)).any()" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "### Finite and Infinite\n", "\n", "Test element-wise for finiteness and infinity" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(1) # 1 is finite" ] }, { "cell_type": "code", "execution_count": 63, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(np.exp(100))" ] }, { "cell_type": "code", "execution_count": 64, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":1: RuntimeWarning: overflow encountered in exp\n", " np.isfinite(np.exp(1000))\n" ] }, { "data": { "text/plain": [ "False" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(np.exp(1000))" ] }, { "cell_type": "code", "execution_count": 65, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":1: RuntimeWarning: overflow encountered in exp\n", " np.exp(1000)\n" ] }, { "data": { "text/plain": [ "inf" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(1000)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Above is a limitation of my computer...\n", "\n", "The answer should actually be true...\n", "\n", "...**but** if we do `np.isinf()'" ] }, { "cell_type": "code", "execution_count": 66, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":1: RuntimeWarning: overflow encountered in exp\n", " np.isinf(np.exp(1000))\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isinf(np.exp(1000))" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(np.inf) # Kind of obvious as infinity is definitely NOT finite" ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":1: RuntimeWarning: invalid value encountered in log\n", " np.isfinite([np.log(-1.), 1., np.log(0)]) # Trying same but in a different way...\n", ":1: RuntimeWarning: divide by zero encountered in log\n", " np.isfinite([np.log(-1.), 1., np.log(0)]) # Trying same but in a different way...\n" ] }, { "data": { "text/plain": [ "array([False, True, False])" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite([np.log(-1.), 1., np.log(0)]) # Trying same but in a different way..." ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### How about `np.NINF` " ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on float object:\n", "\n", "class float(object)\n", " | float(x=0, /)\n", " | \n", " | Convert a string or number to a floating point number, if possible.\n", " | \n", " | Methods defined here:\n", " | \n", " | __abs__(self, /)\n", " | abs(self)\n", " | \n", " | __add__(self, value, /)\n", " | Return self+value.\n", " | \n", " | __bool__(self, /)\n", " | self != 0\n", " | \n", " | __divmod__(self, value, /)\n", " | Return divmod(self, value).\n", " | \n", " | __eq__(self, value, /)\n", " | Return self==value.\n", " | \n", " | __float__(self, /)\n", " | float(self)\n", " | \n", " | __floordiv__(self, value, /)\n", " | Return self//value.\n", " | \n", " | __format__(self, format_spec, /)\n", " | Formats the float according to format_spec.\n", " | \n", " | __ge__(self, value, /)\n", " | Return self>=value.\n", " | \n", " | __getattribute__(self, name, /)\n", " | Return getattr(self, name).\n", " | \n", " | __getnewargs__(self, /)\n", " | \n", " | __gt__(self, value, /)\n", " | Return self>value.\n", " | \n", " | __hash__(self, /)\n", " | Return hash(self).\n", " | \n", " | __int__(self, /)\n", " | int(self)\n", " | \n", " | __le__(self, value, /)\n", " | Return self<=value.\n", " | \n", " | __lt__(self, value, /)\n", " | Return self>> (10.0).as_integer_ratio()\n", " | (10, 1)\n", " | >>> (0.0).as_integer_ratio()\n", " | (0, 1)\n", " | >>> (-.25).as_integer_ratio()\n", " | (-1, 4)\n", " | \n", " | conjugate(self, /)\n", " | Return self, the complex conjugate of any float.\n", " | \n", " | hex(self, /)\n", " | Return a hexadecimal representation of a floating-point number.\n", " | \n", " | >>> (-0.1).hex()\n", " | '-0x1.999999999999ap-4'\n", " | >>> 3.14159.hex()\n", " | '0x1.921f9f01b866ep+1'\n", " | \n", " | is_integer(self, /)\n", " | Return True if the float is an integer.\n", " | \n", " | ----------------------------------------------------------------------\n", " | Class methods defined here:\n", " | \n", " | __getformat__(typestr, /) from builtins.type\n", " | You probably don't want to use this function.\n", " | \n", " | typestr\n", " | Must be 'double' or 'float'.\n", " | \n", " | It exists mainly to be used in Python's test suite.\n", " | \n", " | This function returns whichever of 'unknown', 'IEEE, big-endian' or 'IEEE,\n", " | little-endian' best describes the format of floating point numbers used by the\n", " | C type named by typestr.\n", " | \n", " | __set_format__(typestr, fmt, /) from builtins.type\n", " | You probably don't want to use this function.\n", " | \n", " | typestr\n", " | Must be 'double' or 'float'.\n", " | fmt\n", " | Must be one of 'unknown', 'IEEE, big-endian' or 'IEEE, little-endian',\n", " | and in addition can only be one of the latter two if it appears to\n", " | match the underlying C reality.\n", " | \n", " | It exists mainly to be used in Python's test suite.\n", " | \n", " | Override the automatic determination of C-level floating point type.\n", " | This affects how floats are converted to and from binary strings.\n", " | \n", " | fromhex(string, /) from builtins.type\n", " | Create a floating-point number from a hexadecimal string.\n", " | \n", " | >>> float.fromhex('0x1.ffffp10')\n", " | 2047.984375\n", " | >>> float.fromhex('-0x1p-1074')\n", " | -5e-324\n", " | \n", " | ----------------------------------------------------------------------\n", " | Static methods defined here:\n", " | \n", " | __new__(*args, **kwargs) from builtins.type\n", " | Create and return a new object. See help(type) for accurate signature.\n", " | \n", " | ----------------------------------------------------------------------\n", " | Data descriptors defined here:\n", " | \n", " | imag\n", " | the imaginary part of a complex number\n", " | \n", " | real\n", " | the real part of a complex number\n", "\n" ] } ], "source": [ "help(np.NINF)" ] }, { "cell_type": "code", "execution_count": 70, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(np.NINF)" ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isfinite(np.nan) # Remember, NaN = Not is number" ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 0])" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([-np.inf, 0., np.inf])\n", "y = np.array([2,2,2])\n", "np.isfinite(x, y)" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([-inf, 0., inf])" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 0])" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "## Statistics - A quick brush up" ] }, { "cell_type": "code", "execution_count": 75, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "a = np.array([5,6,7,8])\n", "b = np.array([[1,2,3], [4,5,6]])" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "6.5" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Doing a simple pythonic mean\n", "a.mean()" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "6.5" ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(a) # Numpy gives the same value" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "6.5" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.median(a) # OK, same" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "3.5" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's do for b\n", "np.median(b) " ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([2., 5.])" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# For last axis\n", "np.median(b, axis=-1)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([2.5, 3.5, 4.5])" ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# for second last axis\n", "np.median(b, axis=-2)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### Standard deviation\n" ] }, { "cell_type": "code", "execution_count": 82, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "1.118033988749895" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.std()" ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "1.707825127659933" ] }, "execution_count": 83, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b.std()" ] }, { "cell_type": "code", "execution_count": 84, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "1.118033988749895" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(a)" ] }, { "cell_type": "code", "execution_count": 85, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "2.0615528128088303" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = np.array([[4,5], [8,9]])\n", "np.std(c)" ] }, { "cell_type": "code", "execution_count": 86, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([[4, 5],\n", " [8, 9]])" ] }, "execution_count": 86, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c" ] }, { "cell_type": "code", "execution_count": 87, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([2., 2.])" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, axis=0)" ] }, { "cell_type": "code", "execution_count": 88, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([0.5, 0.5])" ] }, "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, axis=-1)" ] }, { "cell_type": "code", "execution_count": 89, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "array([0.5, 0.5])" ] }, "execution_count": 89, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, axis=1)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "#### But standard can also be inaccurate... " ] }, { "cell_type": "code", "execution_count": 90, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "c = np.zeros((2, 256*256), dtype=np.float16)" ] }, { "cell_type": "code", "execution_count": 91, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "c[0, :] = 1.0" ] }, { "cell_type": "code", "execution_count": 92, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "c[1, :] = 0.1" ] }, { "cell_type": "code", "execution_count": 93, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(d)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "##### What happened up there???\n", "\n", "Now doing same with float32..." ] }, { "cell_type": "code", "execution_count": 94, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "0.4500123" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, dtype=np.float32)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "##### How about `np.float64`? " ] }, { "cell_type": "code", "execution_count": 95, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "0.45001220703125" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, dtype=np.float64)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "##### `np.float128`?" ] }, { "cell_type": "code", "execution_count": 96, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "0.45001220703125" ] }, "execution_count": 96, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.std(c, dtype=np.float128) # as you can see the accuracy doesnt change as you increasing floating point." ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false }, "source": [ "### Working with some data\n", "\n", "Let's look at the population of hares, lynxes and carrots in Northern Canada for a period of 20 years." ] }, { "cell_type": "code", "execution_count": 97, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# year\thare\tlynx\tcarrot\r\n", "1900\t30e3\t4e3\t48300\r\n", "1901\t47.2e3\t6.1e3\t48200\r\n", "1902\t70.2e3\t9.8e3\t41500\r\n", "1903\t77.4e3\t35.2e3\t38200\r\n", "1904\t36.3e3\t59.4e3\t40600\r\n", "1905\t20.6e3\t41.7e3\t39800\r\n", "1906\t18.1e3\t19e3\t38600\r\n", "1907\t21.4e3\t13e3\t42300\r\n", "1908\t22e3\t8.3e3\t44500\r\n", "1909\t25.4e3\t9.1e3\t42100\r\n", "1910\t27.1e3\t7.4e3\t46000\r\n", "1911\t40.3e3\t8e3\t46800\r\n", "1912\t57e3\t12.3e3\t43800\r\n", "1913\t76.6e3\t19.5e3\t40900\r\n", "1914\t52.3e3\t45.7e3\t39400\r\n", "1915\t19.5e3\t51.1e3\t39000\r\n", "1916\t11.2e3\t29.7e3\t36700\r\n", "1917\t7.6e3\t15.8e3\t41800\r\n", "1918\t14.6e3\t9.7e3\t43300\r\n", "1919\t16.2e3\t10.1e3\t41300\r\n", "1920\t24.7e3\t8.6e3\t47300\r\n" ] } ], "source": [ "# Let's view the data\n", "!cat data/populations.txt" ] }, { "cell_type": "code", "execution_count": 98, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "# load the data\n", "data = np.loadtxt('data/populations.txt')" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1900., 30000., 4000., 48300.],\n", " [ 1901., 47200., 6100., 48200.],\n", " [ 1902., 70200., 9800., 41500.],\n", " [ 1903., 77400., 35200., 38200.],\n", " [ 1904., 36300., 59400., 40600.],\n", " [ 1905., 20600., 41700., 39800.],\n", " [ 1906., 18100., 19000., 38600.],\n", " [ 1907., 21400., 13000., 42300.],\n", " [ 1908., 22000., 8300., 44500.],\n", " [ 1909., 25400., 9100., 42100.],\n", " [ 1910., 27100., 7400., 46000.],\n", " [ 1911., 40300., 8000., 46800.],\n", " [ 1912., 57000., 12300., 43800.],\n", " [ 1913., 76600., 19500., 40900.],\n", " [ 1914., 52300., 45700., 39400.],\n", " [ 1915., 19500., 51100., 39000.],\n", " [ 1916., 11200., 29700., 36700.],\n", " [ 1917., 7600., 15800., 41800.],\n", " [ 1918., 14600., 9700., 43300.],\n", " [ 1919., 16200., 10100., 41300.],\n", " [ 1920., 24700., 8600., 47300.]])" ] }, "execution_count": 99, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data" ] }, { "cell_type": "code", "execution_count": 100, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [], "source": [ "year, hares, lynxes, carrots = data.T # What did I just do and why?" ] }, { "cell_type": "code", "execution_count": 101, "metadata": { "hideCode": false, "hidePrompt": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 101, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEGCAYAAAAkHV36AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXiU1fXHP3dmsu+ZhCQkgbAk7EsEZHOpIogruNC6VbRU6oa11tb6U0uLWq1arVtRtCq4b1WpaxW1dUEUBAIKaPYECNlD9mRm7u+P950whElmJrOF5H6eJ8/M3Pe+73uHkDN3zvmec4SUEoVCoVD4D0OwF6BQKBQDHWVoFQqFws8oQ6tQKBR+RhlahUKh8DPK0CoUCoWfUYZWoVAo/IzJnUlCiN8AvwQksAO4HEgDXgISgW+Bn0spO4QQYcA6YBpQA/xMSlmsX+dmYBlgBa6TUn6gjy8EHgSMwJNSyrtdrSkpKUlmZWW5/UYVCoX7bNmypVpKmRzsdQwUXBpaIUQ6cB0wXkrZKoR4BbgAOB14QEr5khDiMTQDulp/rJNSjhZCXAD8FfiZEGK8ft4EYCjwkRAiR7/No8B8oBz4RgixXkr5fW/rysrKYvPmzX14ywqFwhVCiJJgr2Eg4a7rwARECCFMQCSwHzgZeE0/vhZYrD9fpL9GPz5PCCH08ZeklO1SyiIgHzhW/8mXUhZKKTvQdsmLvHtbCoVC0X9waWillHuB+4BSNAPbAGwB6qWUFn1aOZCuP08HyvRzLfp8s+N4t3N6Gj8CIcRyIcRmIcTmqqoqd96fQqFQBB2XhlYIkYC2wxyB9pU/CjjNyVR7Lq/o4Zin40cOSrlGSjldSjk9OVm5jxQKxdGBO66DU4AiKWWVlLIT+BcwB4jXXQkAGcA+/Xk5kAmgH48Dah3Hu53T07hCoVAMCNwxtKXALCFEpO5rnQd8D3wCnK/PWQq8pT9fr79GP/6x1CrXrAcuEEKECSFGANnA18A3QLYQYoQQIhQtYLbe+7emUCgU/QOXqgMp5SYhxGtoEi4LsBVYA7wDvCSEuEMf+6d+yj+BZ4UQ+Wg72Qv063ynKxa+169zjZTSCiCEuBb4AE3e9ZSU8jvfvUWFQqEILuJoLZM4ffp0qeRdCoV/EEJskVJOD/Y6BgoqM0yhUCj8jDK0PkRKya9f2sq/vi0P9lIUCkU/QhlaH/LpD1W8tW0fn+xRGl+FQnEIZWh9yD8+yQegtrk9yCtRfJlfzSd7KoO9DIUCcLOojMI1mwpr+Ka4jhCjoKapI9jLGfTc9d5uWjutnDRmSLCXolCoHa2vePTTApKiQzlz8lBqmpWhDSZWm+SHA42U1rRgtR2dqhrFwEIZWh+wo7yB//1QxbLjRjI0Ppy65g5s6g88aJTUNNNusdFhtVFxsC3Yy1EolKH1BY9+kk9MuIlLZg0jMSoMi01ysK0z2MsatOypaOx6XlzdHMSVKBQaytB6yY8HGnn/uwoum5NFTHgISdGhAMp9EER2OxraGmVoFcFHGVovWf3fAiJCjFw+dwQAiVG6oVUBsaCxp6KR4eZIQk0GSmpagr0chUKpDryhrLaFt7bt47I5WV0G1hwVBiiJVzDZc6CRcamxhBoNynWg6BeoHa0XPP6/AoxCcMXxI7vGzLrroFrtaINCa4eV4ppmxqTGMNwcpVwHin6BMrR9pPJgG69sLue8aRmkxoV3jSdEaoa2Vvlog8KPlY1ICWNTY8gyR1JS06IUIIqgowxtH3ny8yIsVhtXnjjysPFQk4HYcJMytEHCHggbkxrD8KQo2i02DjQqiZciuChD2wfqWzp47qsSzpoylOHmqCOOm6PDqG5SPtpgsKeikfAQA8PNUYzQfzfF1SogpgguytD2gae/KKalw8rVPxnt9Lg5KlTtaIPEnopGsofEYDQIhpsjASXxUgQfZWg9pKndwjNfFjN/fApjUmOczkmMClXyriCxu6KRnBTt9zI0PoIQo1CGVhF0lKH1kOe/KqGhtZNrTnK+mwXNdaASFgJPTVM71U3tjNU/AI0GQWZiJCXKdaAIMsrQekBbp5UnPiviuNFJTM2M73GeOSqUuhZV7yDQ7DlwKBBmZ4SSeCn6AcrQesCrW8qpbmrn6pNG9TrPHB2K1SZpaFX1DgKJvcbBWAdDO9wcRUlNC0drbzzFwMCloRVCjBFCbHP4OSiEuF4IkSiE+FAI8aP+mKDPF0KIh4QQ+UKIPCHEMQ7XWqrP/1EIsdRhfJoQYod+zkN6W/N+RafVxmOfFpA7LJ7ZI829zu1Kw1Xug4Cyp6KRhMgQkmPCusaykiJp7bRS2ahUIIrg4dLQSin3SCmnSimnAtOAFuAN4A/ABillNrBBfw1wGpCt/ywHVgMIIRKBlcBM4Fhgpd0463OWO5y30Cfvzoes37aPvfWtXHvSaFx9DtjTcGuUxCug7K5oZExqzGG/n+FdEi/lPlAED09dB/OAAillCbAIWKuPrwUW688XAeukxldAvBAiDTgV+FBKWSulrAM+BBbqx2KllBul9v1uncO1+g3/2lrOqOQoTh7rumK/PQ1XSbwCh00v9j02NfawcbuWVhWXUQQTTw3tBcCL+vMUKeV+AP3RboHSgTKHc8r1sd7Gy52MH4EQYrkQYrMQYnNVVWAbIBZUNjM1M8Hlbha0YBhAtTK0AaO8rpWWDusRkruh8eGYDIIiFRBTBBG3Da0QIhQ4G3jV1VQnY7IP40cOSrlGSjldSjk9OTnZxTJ8R3O7hYqDbYxMPjILzBkJuqGtVVragLG74iDAEYbWZDRoEi9laBVBxJMd7WnAt1LKA/rrA/rXfvRHe8vRciDT4bwMYJ+L8Qwn4/2GIt2/NzLJPUMbYjQQFxFCjSqVGDDsigN7soIjw82RKg1XEVQ8MbQXcshtALAesCsHlgJvOYxfqqsPZgENumvhA2CBECJBD4ItAD7QjzUKIWbpaoNLHa7VL7Ab2hFu7mhB89Mq1UHg2H2gkczECKLDjiyxnGWOoqSmWUm8FEHDrcLfQohIYD7wK4fhu4FXhBDLgFJgiT7+LnA6kI+mULgcQEpZK4S4HfhGn7dKSlmrP78KeAaIAN7Tf/oNRdXNCKH9wbqLOSpUuQ4CyJ6KRsakxDo9lmWOpLnDSlVTO0Niwp3OUSj8iVuGVkrZApi7jdWgqRC6z5XANT1c5yngKSfjm4GJ7qwlGBRWNTE0LoLwEKPb5yRGhXbthBX+pd1ipai6mYUTUp0eH550SHmgDK0iGKjMMDcoqm52OxBmxxwdpuRdASK/sgmrTfZY5CdLaWkVQUYZWhdIKSmsbmaEm4EwO/ZSiaregf/ZU3FkjQNHMhIiMBqE0tIqgoYytC6obuqgsc3ituLAjjkqFJuEelXvwO/sqWgkxCh6/DAMMRrISIhQWlpF0FCG1gWHFAfRHp2XGK3ScAPF7opGRiVHE2Ls+b/zcF15oFAEA2VoXVBU3QS4r6G1k6QKywQMLfXWudvATpZZq0urJF6KYKAMrQsKq5oJNRkYGh/h0XmJer0D1WnBvzS0dLK/oY0xqc6lXXayzFE0tltUgFIRFJShdUFhdTNZ5kiMBs8qN9pLJdaq7DC/Yi/27XJHm6T6hymChzK0Lijqg+IAIDFSuQ4CwZ4eahx0Z7jqiKsIIsrQ9oLVJimpaWZEkmeBMNCKmcRHhijXgZ/ZXdFITLiJtLjeExEyEyIxCFRATBEUlKHthfK6Fjqt0uNkBTuq7bj/2VOhBcJcla8MNRlIT4igWGlpFUFAGdpeKPSwald3zFFhVCt5l9+QUrLnQKNLt4GdLNWoUREklKHthaIqXUPbR0ObqHa0fmVfQxuNbRaXigM7w82RFFWrKl6KwKMMbS8UVjcRFxHSpSDwFFUq0b/YA2GuFAd2ssxRNLZZqG9R2XqKwKIMbS/YFQd9bcprjgqlrqUDq6p34Bd291Ls2xldxWWU+0ARYJSh7YWiquY++2dBq+AlJdS3qF2tP9hT0cjQuHDiIkLcmq+0tIpgoQxtD7R0WNjX4H6fMGckqjRcv7KnopEcN90GABkJkQihtLSKwKMMbQ/Y/xj7oqG1Y1ZpuH6j02qjoKrJbcUBQHiIkaFxEUpLqwg4ytD2QFfVLm9cB1F6BS+VhutzCqua6bRKtwNhdrKSIpWWVhFwlKHtAXvVLrtfry8cqnegdrS+pqu9eA99wnpiuNLSKoKAW4ZWCBEvhHhNCLFbCLFLCDFbCJEohPhQCPGj/pigzxVCiIeEEPlCiDwhxDEO11mqz/9RCLHUYXyaEGKHfs5Doq9hfh9SWNXM0LhwIkPdaqvmlITIEITQiocrfMueikaMBsGoIZ5948gyR1Lf0qkClIqA4u6O9kHgfSnlWGAKsAv4A7BBSpkNbNBfA5wGZOs/y4HVAEKIRGAlMBM4FlhpN876nOUO5y307m15T2F1s0ftxZ1hMhqIjwhRFbz8wJ6KRkYmRRFmcr9hJhwqLqPa2igCiUtDK4SIBU4A/gkgpeyQUtYDi4C1+rS1wGL9+SJgndT4CogXQqQBpwIfSilrpZR1wIfAQv1YrJRyo95Bd53DtYKClJLCqiav/LN2zNFhKhjmBzxJvXXE/jtV7gNFIHFnRzsSqAKeFkJsFUI8KYSIAlKklPsB9Mch+vx0oMzh/HJ9rLfxcifjQaO2uYODbRZGeqE4sJMYpbLDfE1Tu4XyulaPA2EAwxJ1La2SeCkCiDuG1gQcA6yWUuYCzRxyEzjDmX9V9mH8yAsLsVwIsVkIsbmqqqr3VXvBoT5h3u9ok6J7qHew/SVoqfX6+oORQ11vPQuEgSbxSosLVxIvRUBxx9CWA+VSyk3669fQDO8B/Ws/+mOlw/xMh/MzgH0uxjOcjB+BlHKNlHK6lHJ6cnKyG0vvG95W7XIkMSr0yAaNdcXwxq9g8z+9vv5gxG5o+7KjBa24jHIdKAKJS0MrpawAyoQQY/ShecD3wHrArhxYCrylP18PXKqrD2YBDbpr4QNggRAiQQ+CLQA+0I81CiFm6WqDSx2uFRQKq5oJMQrSPewT5ozEqDDqWzuxWG2HBmsKtMf9272+/mBkT8VBokKNff79ZJmjVDBMEVDc1S6tAJ4XQoQChcDlaEb6FSHEMqAUWKLPfRc4HcgHWvS5SClrhRC3A9/o81ZJKe3fna8CngEigPf0n6BRVN3EcHMUpl7aV7tLUnQoUkJdSyfJMVoCA7WF2uP+PK+vPxjZrafeGjzs42YnKymKmuYODrZ1EhvuXp0EhcIb3DK0UsptwHQnh+Y5mSuBa3q4zlPAU07GNwMT3VlLIOhrnzBnOCYtHGFo60ugtQ4iEno4W9Ede7Hv0yam9vkaWWYtIFZS3cKkjDhfLU2h6BGVGdYNq01SXNPiE/8sOKThOvppawroigGqXa1HVDa2U9/SyRg3SyM6Y7gql6gIMMrQdmNffSsdFptXVbsc6Sos46g8qC2E4XO058pP6xFdNWj7GAgDLRgGqlGjInD0Pb90gFLYVUzGew0taMW/waHegc2qqQ7GngF1JVChdrSecKirgufSLjuRoSZSYsMoUlpar9myZcsQk8n0JJrrbzBv3GzATovF8stp06ZVdj+oDG03iqq0YjK+8tHGR4YihIProKEMbJ2QOBLSpqgdrYd8t+8gKbFhfW4vZGe4OUrtaH2AyWR6MjU1dVxycnKdwWAYtK1EbDabqKqqGl9RUfEkcHb344P5E8gphdXNxISbSIr27g/ZjtEgSIh0yA6zB8LMozRDW/0jtDf55F6DgbzyBiZnxHt9nSyzKpfoIyYmJycfHMxGFsBgMMjk5OQGegjqK0PbjaJqrX2NLwuImaNCD9U7sGtoE0dC2mRAwoGdPrvXQKahtZOi6mam+EApMNwcRXVTO03tFh+sbFBjGOxG1o7+7+DUpipD243CKt9Ju+wc1na8tghMERCdqu1oQSkP3GTn3gYAn+xou4rLVCv3wdFOZGRkruPrhx56yHzppZcOC9Z6nKEMrQNtnVb2NbT6LBBmJyk6jGp7qcTaQm03azBATBpEJSs/rZtsL68HYLJPdrR25YFyHwx2Ojv9335eGVoHimuakRKfSbvsHL6jLYDEEdpzIVRAzAPyyhoYlhhJfKT3/nOlpR0cvPDCC3GTJ08eO27cuPFz5szJKSsrMwHccMMNQy+88MLhc+fOzT733HNHWCwWfvWrX2VMnDhxXE5Ozvh77703yZfrUKoDB4qqvO8T5gxzdCj1LZ1YOjsx1RXDmNMOHUydDIUPgaUdTGE+ve9AY8feBnKHee82AIgOM5EUHaaUBz7kd69tz/yhorHvvZ+ckJMa03Lv+VPKepvT3t5uGDt27Hj764aGBuP8+fMbAObPn990wQUX7DYYDNx///1Jq1atSn3iiSfKAfLy8iI3bdq0Ozo6Wt53331JcXFx1p07d+5qbW0VM2bMGHvWWWcdHDt2rE9qnCpD60ChDxoyOsOupa0/UESStUNzHdhJmwI2C1R+D0Nze7iCorqpnb31rVw2J8tn1xyRFKnq0g4AwsLCbLt37/7e/vqhhx4yb968OQqgqKgodPHixRlVVVUhHR0dhszMzK4UzYULF9ZHR0dLgI8++ih29+7dkevXr08AaGxsNH7//ffhytD6gcKqZlJjw4kK8+0/S6KehttS8aM+0M3QguY+UIa2R/J0/6wvaxMMN0fxvx/8V9d4sOFq5xkMrr322mG//vWvKy6++OKGt99+O2bVqlVD7ceioqK6SupJKcXf/va30vPOO++gP9ahfLQOFFX7pn1Nd+xpuJ2VdmnXqEMHE7IgLE75aV2wvawBIWBiuu8MbZY5ksrGdlo6lMRroNLY2GgcNmxYJ8Azzzxj7mne/PnzG1avXp3c3t4uAPLy8sIOHjzoM/uoDK0DRT5oyOgMu+uA2gIwhWtqAztCaHpaZWh7ZcfeBkYnRxPtw28bqlHjwOeWW27Zd+GFF46aNm3aGLPZ3OMn6m9+85vqsWPHtk2aNGlcdnb2hCuuuGJ4Z2enz8T0ynWgU9fcQV1Lp8+qdjlijtZcByENxYekXY6kTYFvngSrBYzqV9IdKSV55fWcmDPE9WQPcNTSjkvre+0ERXBpaWnZ6vj6uuuuqwFqAC655JL6Sy65pL77Offff/9hXVyMRiOPPPLIXmCvP9aodrQ6Xe1r/LCjjY8IwSAgsqnkcP+snbQpYGmD6h98fu+BwL6GNqqbOnyin3VkmK6lVam4Cn+jDK1OkY+rdjliMAjMkSbi2vYe0tA64hgQUxxBXpnvEhUciQ0PITEqlNJaJfFS+BdlaHUKq5owGQQZCd73CXNGTmQjIbLj8ECYHfNoCIlUhrYH8vY2YDIIv3y9T4sLp6KhzefXVSgcUYZWp6i6mWHmSEJ80CfMGeNC9RKVzlwHBiOkTFSGtgfyyusZmxZDeIjR59dOiQ2n4mC764kKhRcoQ6tjr9rlL0YZdUNrdrKjBc19ULEDbDbnxwcpNpskr7yBSem+yQjrTkpsOAcOutjRtjXACz+D8i1+WYNi4OOWoRVCFAshdgghtgkhNutjiUKID4UQP+qPCfq4EEI8JITIF0LkCSGOcbjOUn3+j0KIpQ7j0/Tr5+vn+q5GoRvYbNKnDRmdkUkF7YRAzFDnE9KmQEcj1BX5bQ1HIyW1LTS2WXxSGtEZqbHh1DZ30G6x9jzpw5Xww/uw5x2/rEEx8PFkR3uSlHKqlNLeDfcPwAYpZTawQX8NcBqQrf8sB1aDZpiBlcBM4Fhgpd0463OWO5y3sM/vqA/sa2il3WJjZLLvA2F20iz7KLENobOnyp1dAbFtflvD0UheV8Uu/+xoU+M06V1lT+6Dki9hy9Pa86o9flmDwju6l0nsj3jjOlgErNWfrwUWO4yvkxpfAfFCiDTgVOBDKWWtlLIO+BBYqB+LlVJu1FuVr3O4VkAo8lONA0fMHeUUy1TqmntInU4eC4YQ5aftxvayBsJMBnJS/PMhmBIbDkCFM/eBpR3+/WuIHwYjT1LyO0WfcdfQSuA/QogtQojl+liKlHI/gP5oV5OnA445z+X6WG/j5U7Gj0AIsVwIsVkIsbmqync56nZD6zcfrc1GTItmaKubejC0plBIGa+KgHcjr7yeCUNjMfkpSJkapxlap37az/6mGdcz/67VoagtBKv/a5cqvKOurs6Qnp4+yZ5OW1tb2/X62GOPHXPVVVelT5o0aVxWVtbE999/PxrgT3/6U8qSJUuyAL7++uuI7OzsCY2NjT77T+duGtJcKeU+IcQQ4EMhxO5e5jrzr8o+jB85KOUaYA3A9OnTfdY+o7CqmegwE8kxfipT2LgPo62dYpnK+J52tKC5D3a9DVJqqbmDHIvVxnf7DvKzGZl+u0eqfUfbXeJVuQs+ux8m/wxGz4PmKq3KWm0hJI/x23qOat68JpPK731aJpEh41tY/KhHxWoSEhJss2fPbnzllVfifv7zn9c/9dRTiaeffnpdWFiYBLBYLGLHjh27Xn755bhVq1YNXbhw4Q+33XbbgZkzZ45Zt25d/D333JP26KOPFsfExPgsMu2WxZZS7tMfK4E30HysB/Sv/eiP9ha75YDjX0YGsM/FeIaT8YBRqAfC/BaD0xsyFssUapp7kRKlTYHWWmgo73nOICK/qonWTitTMv0TCAOIiwgh1GQ4fEdrs8H66yAsBk79izZmN67KT3tUsHz58ip7EZnnnnsuafny5dX2Y0uWLKkDmDNnTnN5eXkoaCm469atK7ryyitHzJ49u3HBggU+zWJxuaMVQkQBBillo/58AbAKWA8sBe7WH9/ST1kPXCuEeAkt8NUgpdwvhPgA+ItDAGwBcLOUslYI0SiEmAVsAi4FHvbdW3RNUXUTuZkJrif2Fb0hY4kt5VCTRmekOmSIxftvF3e0kFfmux5hPSGEILW7lnbzP6H8azjncYjSC+0n5WiP1crQ9oiHO09/smDBguYVK1aEvfPOO9FWq1XMmDGj65M0PDxcAphMJqxWa9fuateuXeGRkZG2ioqKEF+vxx3XQQrwhr7bMwEvSCnfF0J8A7wihFgGlAJL9PnvAqcD+UALcDmAblBvB77R562SUtbqz68CngEigPf0H6/4tOxT7t9yP2HGMEINoYQaHX4Modq4MRSjCKHCamNE0s+8vWXP1BYijWEcMCQdamnjjJQJIAxQkQfjzvTfeo4S8vbWExNmYoTZf0FK0NwHB+yug4a98NGfteDXZIf/E6FREJepdrRHERdccEHN5ZdfPvK3v/3tfldza2pqjDfeeGPmxx9/vPuqq64a9vTTTydcfvnldb5ai0tDK6UsBKY4Ga8B5jkZl8A1PVzrKeApJ+Ob6aEfel+JDokmOz6bDmsHHbYO2q3tNHU0dT3vsHbQYe2gubOF8KHNhETOAXJ8uYRD1BYiErKIF+G9uw5CIyFpjFIe6OSVNzAxPQ6Dwb/+6pS4cLaX1Wu+8Xdv1HyxZz5wpJ88KUcZ2n5IW1ubISUlZbL99VVXXXXgT3/604Fly5bV/PWvf01ftmxZbW/nA1x55ZWZy5Ytq5o8eXL72rVri08++eQxCxYsaExPT/dJseIBW5Nveup0pqdOdznvze2F3LL5Ir6sfYFrj/zc8A21hWAehbkztHfXAWh+2qL/+mcdRxHtFiu79h/kF8c5KcLjY1Jjw/jgYBvy+7cQe96F+bc7L/6TPEbT1dpsR5a6PAqoa6sjIdyPLrIgYbPZnKbsbdiwIWbhwoV1SUlJXdkoX3/9ddcnZVpammXv3r07AF599dVi+/jo0aM7S0tLd/pyjUff/xYfs6/ORkfNieyo/ZotB/yQYmmzQW0RJI7EHB1KTW+uA9CKgDfuh6bK3ucNcHbvb6TTKpniR/+snZTYcMItB5Hv/k77oJt1tfOJSTlgaYWGfuOKdJumjibOeesc/rHtH8FeSkBYunRp5sqVK9NXrVoV0MB6Twx6Q7t7fyNJtp+QFJHEI1sfQfN8+JDG/dofZ+KIw9uO90RXhtjg0NNKKSk9WMqb+W/yzM5naOnUasPm7dUCYZN82LqmJ1LjwvmD6UVESw2c9VDPxdeTx2qPR2HiwhM7nqCmrYYTMk4I9lICwtq1a8tKS0t3Tp48uV9UDBqwrgN32bm3gUlDkzlx0i+5++u72VSxiVlps3x3A13aReIozFGhVDe5+L2nTtIe92+D7FN8t45+Qqetk901u9laubXrp6atpuv4+8Xv88i8R8grqycxKtRvZSsdGd2yjbGmTygds4xhQ6f2PNFR4pU93+/r8hVlB8t49vtnOXvU2UxM8mkoROEmg9rQNrZ1UljdzDm56SzJWcIz3z3Dw1sfZmbqTN9parsM7UjM0R00tlnosNgINfXwZSI8TiulOEACYlJKNh/YzFf7v2Jr5VZ2VO2gzapF+NOj05k9dDa5Q3LJHZJLeWM5N312Exe/czGd+5YxOSPLf9pmO51tjPrqVkpsQ/gm61cM621uZCJEJkFVb/k6/Y+/bfkbJoOJXx/z62AvZdAyqA3td/u0zsITM+IINYayfPJyVm1cxWd7P/PdV6zaAjCGQlwGiVFaIkJdS0dXjr1T0qbA3m99c/8gs+77ddy3+T6MwsiYxDGcl3Nel2EdEnl4D7DshGyePvVprt5wDTUx9zI9+f/8v8D/3UtIfQG3WG5mRpMbnrTkMUeV6+Dr/V+zoXQDK3JXHPHvrQgcg9pHu7ObH3Dx6MVkRGf41ldbW6i1FDcYSdLbjrt2H0yG+hJo9ZmMLyi8V/Qe922+j/nD5/PlhV/y8pkv84dj/8CpWaf2+Ec/IWkCf5jyKDZLLP+pXcXbhW/7d5Gb/wnjzub7iGnOC8t0xy7x8rUv3w9YbVbu+eYehkYN5dLxlwZ7OYOaQW1o88obGBoXTpK9S60hhCunXMmu2l18XPqxb25SU9jVVSExSruP2wGxih2+WUMQ+KbiG275/BaOGXIMdx1/F5Eh7qfA762KoKX4SiaaJ3PzZzezJm+N74OUoH2QtdZB5rHuFQAHbUfbVq/VPujn/Cv/X+yp28MN028g3NTLN6gBQGlpqenMM88cmZmZOXHUqFETTjzxxNF5eXk+L16yatWqIX0pNjOoDe3OvQEfbn0AACAASURBVJog3pEzRp5BVmwWj2x7BJv0sqaElNqOVu8TZtZ3tG5paeGo9dMW1Bfw609+TUZMBg+d/BBhRs/+v+eVN5AWk8hTC5/g9BGn8/DWh/nzxj/TafNx5ay6Eu0xfjipsWHuG1ro94kLjR2NPLL1EY4ZcgwLhi8I9nL8is1m4+yzzx59wgknNJaVle0sKCj47q677tq7b98+l6m0NpsNq/Xwou8WS885Co8//nhKU5M7PqbDGbSG1h4I6y4fMhlMXD31avLr8/mg+AMvb1LRJe0CMEfphtbVjjYqCWIzgmJot1Vu47z157Hsg2XsbfK8xX1lSyVXfXQVYcYwVp+ymrgwz+VZO/Y2MFn3m999/N1cMekKXv/xdVZsWEFzpw9rfdTrhjYhi9Q4N3e0Sbqh7ec1Dx7f/jh1bXXcdOxN/g8oBpm33347xmQyyd///vddXzPmzJnTOnv27JbZs2fnjB8/flxOTs745557Lh5gz549oSNHjpxwySWXDJswYcL4goKC0MjIyNzrr79+6OTJk8du2LAh+q233ooZN27c+JycnPFLlizJam1tFXfccceQysrKkBNPPDFn5syZHqWRDtpgmGMgrDunZp3Kmrw1/GPbP5g/fD4mQx//mWq1YjL2PmGx4SEYDYIaVz5a0BIXAqilbepo4sFvH+TlPS+TEpXCvqZ9nLf+PP5w7B9YNGqRW3+sTR1NXP3R1dS31/PMwmdIj3ZaVrhXGlo7Kapu5vxpWkE3IQTXHXMd6dHp3P7V7Sx9bymPznuUlKgUj699BHXF2mPCcFJiD1Dd1NG7IgQgdiiExvTrHW3JwRKe3/08i0cvZrx5fMDue9sXt2Xm1+X7tEzi6ITRLbfPvb3XDJG8vLyIKVOmtHQfj4yMtL3zzjv5iYmJtv3795tmzpw59qKLLqoHKC4uDn/iiSeKn3vuuVKA1tZWw8SJE1v//ve/72tpaREjR46c9J///GfP5MmT288555yse++9N/mPf/xj5erVq1P++9///pCWluZRau6g3dF2D4Q5YhAGrp16LcUHi3mn0Is+UQ7SLgCDQbiXtACa+6D6B+jwabU2p3xa9imL31rMy3te5uJxF/PWord4/ezXGZc4jtu+uI0bPr2BurbeA3Odtk5u+PQG8uvzuf8n9/f5D3xHub1i1+G/l/NyzuPReY9S1ljGRe9exO5aH0is6oohIgHC47pUIJWNLna1QkBSdr82tPdtvo8wYxjXHXNdsJcSVGw2m7j++uszcnJyxp900kk5lZWVoeXl5SaAtLS0jnnz5nX9cRmNRi677LI6gO3bt4dnZGS025MdLrvssprPP/88xpu1DNodbfdAWHdOHnYy4xLHsXr7ak4feTohhj5UTqsp0NrTxB4qt2uOciMNF3Q/rYSKnTBspuf3doPq1mru/vpuPij+gNHxo7n/J/czOVmrzREZEsmTC57k2e+f5aGtD3Hu+nNZNWcVx2ccf8R1pJT86cs/sXH/RlbNWcVx6cf1eU3b7T3CnHS9nZs+l3WnreOaDddw6XuXcvfxd3PysJP7fC/qSjRFCIcKgB842EZGgotNWfIYKPy07/d1gw5rBx8Uf0CrpZXFoxcTagx167yN+zbyadmnXH/M9SRFJPl1jd1xtfP0F5MmTWp98803jyji8PjjjyfW1NSYduzYsSssLEymp6dPam1tNYC223WcGxoaajOZNHPoj8DroN7Rdg+EOSKE4Nrca9nbtJc389/s203s0i6HlE5zdKibrgP/BcSklLzx4xssenMRH5d+zIrcFbxy5itdRtaO0WDksomX8eIZLxIfFs/VG67mjq/uoNXSeti8R7c9yvqC9Vw95WrOyT7Hq7XtKG8gyxxJXKTzD7YxiWN48YwXGRU3ius/uZ6ndz7d9z+MumKIHw449A5rcON3k5SjpVa3NfTtvr1Q3Vqtuaxem8//ff5/3P7V7Sx+azGflH7i8n1abBbu+eYe0qPTuWT8JT5fW3/lrLPOauzo6BB/+9vfuj5Z/vvf/0aWlJSEJiUldYaFhcl///vfMfv27XPr02rq1Klte/fuDd25c2cYwLp168zHH398I0BUVJS1oaFBBcPcoadAWHeOTz+eycmTeXz747Rb+5AyrReTcSQxKsw910FMmpaF5GNDW3qwlCv+cwV//PKPZCdk8/rZr7N88nJCjD3v2MckjuGlM1/i0vGX8vKel/npv3/KzmqtuNFrP7zG43mPc87oc7hyypVery+vvJ5JLgrJJEcm8/TCp1mQtYD7t9zPyi9X0ulpLy+bFepLD+1o43pp0njEAuw1D3707J698H3N99zy+S0seG0Bq7evZmLSRNbMX8PjpzxOiCGE6z65jis/upKC+oIer/HaD6+RX5/PjdNv9FjpcTRjMBhYv359wYYNG2IzMzMnjh49esLKlSuHnn322Q3bt2+Pmjhx4rjnnnsuccSIEW78ciEyMlI+9thjxUuWLBmVk5Mz3mAwcOONN1YBLF26tPq0007LVsEwN+gtEOaIEIIVuSu44j9X8NoPr3HxuIvdv4ld2jXi8K/a5ig3SiVqN9d2tRW+M7Rv/PgGd266kxBDCH+c/UfOyz4Pg3DvszbMGMbvZvyOEzJO4NYvbuWSdy/h7FFns75gPXPT53Lb7Nu8jm5XNbazr6GNX7j4vQCEm8K554R7yIrN4vG8xylrLOOBnzxAfLib1b4a94Ots8vQJkRqLW0qPZV4ZbguxdkTVpuVT8o+4bldz7HlwBYiTBGcn3M+F429iKy4rK55r6W9xsu7X+Yf2/7BeevP48KxF3LV1KuIDY3tmtPQ3sCj2x5lRuoM5g3zU7nPfkxWVlbnu+++W9h9fNu2bU6d+T/++ON3jq9bWlq2Or5etGhR46JFi77vft4tt9xSecstt3hcWm9Q7mh7C4R1Z2bqTKanTOeJvCeO+MrcK00HoLP5iB2tOSqUxnYL7RZrDyc6kDZFaxJo8b4A0TuF77Dyy5XkDsnlrcVvsSRnidtG1pGZaTN5/ezXWThiIW/kv0FOQg73n3h/33zY3ciz+2fdLI1oEAauzb2Wu46/i+1V27n43Yspaihy72YOigPQPlRTYsPc29HGD9fSqvso8Wq3trP2u7Wc8cYZ/ObT37C/aT83Tr+Rj5Z8xP/N/L/DjCxoiTSXjL+Et899m3Ozz+X5Xc9z5r/O5NUfXsVq0/4fPbb9MRraG7hpxsCXcx2NDEpD6yoQ5ojdV1vTVsOLu190/yZ6n7DuBaQT9aQFt5UHNgtUHvHB6hH/K/8ft35+K9NTp/PIvEe8znmPDY3l7uPvZt1p63hiwRMeZX31Rl55AwYBE4bGup7swJkjz+SpU5+iqbOJi9+9mK/2f+X6pC5Dm9U1lBobfmQ3XGcYTWAe3WflwZ+//DP3bb6P1KhUHvjJA7xz7jssnbD0sB2qMxLDE/nj7D/yylmvMCJuBKs2ruKCdy7grfy3eGn3S5ybfS5jElWH3v7IoDS0rgJh3ZmWMo0TMk7g4W8f5j/F/3HvJIfyiI6Y9TRct9wHPvAFbjmwhRs+vYGcxBweOsnzLK3eyB2S26eEhJ7IK69n9JBoosI892hNHTKVF854gZTIFK768Cpe/eHV3k+oK9H6s8UdaoI5xN00XOhzW5uPSz/m34X/Zvnk5Tyz8BlOGX6KxzrtsYljeWbhM9x74r3Ut9dz6xe3Em4KZ0XuCo/XowgMbhtaIYRRCLFVCPG2/nqEEGKTEOJHIcTLQohQfTxMf52vH89yuMbN+vgeIcSpDuML9bF8IcQffPf2jsTdQFh37j7+biYlT+J3//sd/y74t+sTagvBYDrsDxkOpeG6taNNHKEZgz4a2l01u7h2w7UMjR7K6lNWEx0a3afrBAIpJXnlDV51vE2PTufZ055l1tBZrNq4ils+v4XChiPcdhp1xZrsziEIqHXDbXNPxZA8Rsss63TTMAP1bfWs2riKsYljuXKyd4FDIQQLsxayfvF6bph2A3cdfxfmCLNX1+wjNpvNpnwVaLpdwGnevic72l8Duxxe/xV4QEqZDdQBy/TxZUCdlHI08IA+DyHEeOACYAKwEPiHbryNwKPAacB44EJ9rl9wNxDWnZjQGB475TFmpMzgls9vcb1jqi04QtoFmo/WGFHMX7evYMXHK3hq51Nsq9xGh9WJ4TWFaf7AGs8NbXFDMVd+dCUxoTGsmb+GxPBEj68RSPY1tFHT3MEUD38v3YkOjebhkx/m8omX837R+yx6cxFXfnQlX+z94nADWl/S5Z+1kxobTlunjYOtbiT9JI8BaYOafLfXduemO2noaOCOuXf0qvLwhAhTBJdPvJyfZP7EJ9frAzurqqriBruxtdlsoqqqKg5w2mvMre8sQogM4AzgTuAGoXnbTwYu0qesBf4ErAYW6c8BXgMe0ecvAl6SUrYDRUKIfOBYfV6+3m0XIcRL+lzvHJM94EkgrDuRIZE8Mu8Rbvj0BlZtXEWbpY2fj/+588m1hUcEwjqtnbxWuIaI4c9Q227G0NDMp2WfAhBqCGVi0kSmDpnKMUOOYeqQqdrX8qRsqHb/jxmgormC5R8uB2DN/DWkRqV6/F4DzfYyLRDmStrlDiaDiRum3cBlEy7j1T2v8tKel7jyoysZGTeSi8ddzJkjzySyrviILgkpusTrQGNbjzreLhxrHqS67lrwQfEHvF/8PityVwwoP6rFYvllRUXFkxUVFRMZpK5IHRuw02Kx/NLZQXedQ38Hfg/Y09DMQL2U0v7RXw7YE9vTgTIAKaVFCNGgz08HHKMUjueUdRv3TyoUngXCnBFuCufBkx7kps9u4p5v7qHN0sYVk684fJKUmoZ2+NyuoYL6Am7+7GZ21e7C2jCDM7JXcOvpudS01rCtahtbD2htXdZ9t46ndmod2UfFjeKYEAunNJdyrKUDk8m13rq2rZblHy6nsaORp0596ogIdn9lY0ENkaFGxqd5FgjrjcTwRH415Vf8YuIveL/4fZ7b9Ry3f3U7D377IOeFtHNhTBJpDvNTu5IW2shJcZFxaR6tuXWqXBcBr2mt4c6v7mSCeQK/mPgLL95R/2PatGmVwNnBXkd/x6WhFUKcCVRKKbcIIX5iH3YyVbo41tO4s09Bp04yIcRyYDnAsGG9Nh3pEU8DYc4IMYZwzwn3cNsXt/HQ1odotbSyInfFIVlNUyV0NEHiKGzSxou7X+SBLQ8QaYrk7yf9nZuflTS2GAEwR5iZN2xel/ax1dLKzuqdbKvcxreV3/Lu/q95NTmOxNfmsSBrIWeMPIMpyVOcSniaOpq46qOr2Ne0j8fnP8448ziv3mcg+aKgmmNHJPZe0KWPhBhDOGvUWZw58ky2VW3j2a3/YG17A+vK3mDep02ckHECaVFp2IxxIDrdk3iFhGtuHRdtbaSU3P7V7TR3NnPncXf2vUCR4qjGnd/6XOBsIcTpQDgQi7bDjRdCmPRdbQZgb+tbDmQC5UIIExAH1DqM23E8p6fxw5BSrgHWAEyfPt3jvEvHHmHeYjKYuPO4OwkzhvHEDk1j+/sZv9cMoK44OBCVyG0fXsnG/Rs5IeME/jznzyRFJJEY9T9qmp1rYyNMEcxIncGM1BkAtBds4PPXL+ad8VN5I/8NXtrzEkOjhrJwxEJOH3E6OQk5CCFos7Sx4uMV/FD7Aw+e/CDTUqZ5/R4DRUVDG4VVzVw4o28fnu4ihNDa6GSdx75vXueluct4bf8X/KfkkJIkZizctyueNyoySItOIyUyhbSoNNKj05mTPocIk0OzSDfa2rxb9C4bSjdww7QbGBU/qte5ioGLS0MrpbwZuBlA39HeKKW8WAjxKnA+8BKwFHhLP2W9/nqjfvxjKaUUQqwHXhBC3A8MBbKBr9F2utlCiBHAXrSAmd3361P6GgjrCYMwsHL2SiJMETy36znare3cOutWDLUFvB8Vye15f6NTWvnj7D9yfvb5XbvQpOgw9wrLAGHJ45nX0sq8ISfQtGA1n5R9wjtF77D2u7U8tfMpRsWN4rQRp7GjegdbDmzh7uPvPupaSn9ZUA3AnNEBiprXFTPUYuWGGb9nRUQ8+5v3U9Fcwf7m/dz2zuckJ7cTG9ZOQX0Bn+/9vCtRZVjMMG6fezvHpByjXScpBwo+BqvFaYvyypZK/rLpL0xJnqJayQxyvPkecxPwkhDiDmAr8E99/J/As3qwqxbNcCKl/E4I8QpakMsCXCOltAIIIa4FPgCMwFNSysPS43yFN4GwnhBC8PsZvyfCFNG1s5UVebwzJInJcSP4y/F3MTz28Oh2YlQopbVHlM90TkwqhEZD9Y9Eh0Zz1qizOGvUWdS21fJh8Ye8W/Quj2x7BIBbZ97K6SNP99l7CxRf5NeQEBnCuFTf+Wd7pb4EQqIgKokQIRgWO4xhsdpuerVMIK0zgsfna98opJQc7DjI9qrt/GXTX7js/cu4eNzFXHfMdUQkjwVrh3Y98+G7VSklf974ZzqsHdwx9w6MBmNg3puiX+KRoZVSfgp8qj8v5JBqwHFOG7Ckh/PvRFMudB9/F3jXk7X0BW8DYT1hL04dbgrn4a0PYwSubjNwxWnrnPrkzNFu1qTVLq4FXrrJiBLDE/nZ2J/xs7E/Y3/TfvY37z+00zqKkFLyZUE1s0eZMRgCpBCqK9akXU783Cm6ltaOEIK4sDhOyDiB6SnTuX/L/Ty36zk+2/sZt2dfTC5oiQvdDO2b+W/yv/L/cdOMm46agKTCfwwqOYYvAmG9sXzych486UFeaIvmqqjsHgMf5qhQmtottHW6Ue8AdEPbs5Y2LTrtqDSyAMU1LexvaGPOqADWTnWoQ9ud1NhwDhx07j+PDInk1lm38uSCJ7HYLCzdchf3JsbTVnn4F7CK5gru+eYepqVM46JxfvGCKY4yBo2h7WtGmKecnHkS42vKjtDQOmKOdrMbrp2kbKgvg04PitocJXyRr/ln544OkKGV8rA6tN1JiQunuqmdTmvPjTnthXV+OuanrIuLZUnJq2yr3KZfXrLyy5VYpZXb597ep8I9ioHHoPlf4OtAWI80V0FH4xFfJR1JjPIgDRe0HS3yUP2EAcSXBdWkxYWTZfZpq6meaa7Wqqr1sqOVUivZ2BtRIVHcOutWniCVDms7l753Kfd9cx/P73qeL/d9yW+n/ZbMmMxer6EYPAwaQ+uPQJhTuvUJc4a9G261O50WQNvRgk8LTfcHbDbJxoIa5oxKClxpP4fOt85IjdO+bbilpQVmJU/lX/trWJJzPmu/X8tfv/krs9Jm8dMxP/XFahUDhEGjnvZXIOwI3DG0nroO7BXAPMirPxrYVXGQupZO5gZK1gVH1KHtzpAYPQ3XnXKJAMljiGpv5Lbxy5iftYBX97zK72b8TtWEVRzGoDG0/g6EdVFTAMII8T2L71NjwzEIKKp2s8NtWDTEDB1whvbL/BqAAAfC9MLgPfhoPWppA4fVPJg16mRmpc3ydoWKAcigcB0EKhAGaDva+GGHld/rTkSokbGpsWzTC6m4RdLoAec6+KKgmpHJUV3GLSDUlUDUEAh17hNOjAwlxCjcN7RdbW1c1zxQDF4GhaENWCAMtPKIvQTC7OQOi2dbaT02m5uZxOZsTeLlh1bIwaDTauProlrmBnI3C7qGNqvHwwaDYEhMOJU9SLyOICoZwuP73NZGMTgYFIY2YIEwe9WuhBEup+YOS6Cx3UJBVZN7107K1tpbN1d7ucj+wfayelo6rIH1z4Jehzar1ympcW62tAEt6SF5TJ/b2igGB4PC0AYsENZaB+0Hj+gT5ozcYVrd1W9L69y7tllXHgwQP+0X+TUIAbNGBtDQWjuhobzHQJidVE9a2kCf29ooBg+DwtAGLBBWqwda3NjRjjBHERcRwtZSN/20dndEH7ot9Ee+KKhmwtBY4iNd19j1GQ1lWlcEFzvaFE9a2oDW262lGlpqvV+jYkAy4A1tQANh9oi2Gztag0EwNTPefUMbPwyMYQMiINbaYWVraV0Q/LO9a2jtpMaF0dJhpbHdjZY24BAQU7tahXMGvKENaCDMhXSoO7nD4vmhspHGtk7Xkw1GTZs7AFwH3xTX0mmVzAlU2q0du4bWxe8nJdZDLW1SjvaoAmKKHhjwhjZggTCA2mKITu1ROtSd3GEJSKn5kN1igEi8viioJsQomJGVENgb1xWDIQRih/Y6zW5o3ZZ4xWVCSKSSeCl6ZMAb2oAFwkDb0brhNrAzNVMLiG31JCBWV6QVmj6K2VhQQ25mApGhAc6XqS+B+Ezt20Ev2HuH9VTF6wgMBq0ehYu2NorBy4A3tAELhIFLjWZ34iJCGD0k2oOA2GiwWQ7l6x+FNLR0smNvQ+C6KTji5u/HnkDhkfLAjbY2isHLgDa0AQ2EdbbBwX1uKQ4cyc2MZ2tZvXsR7gFQXGZjYQ1SBrAsoiO91KF1JDzESFxEiPtaWtAMbUMZtLupi1YMKga0oQ1oIKy+BJAeuQ5A89PWNne419rGPFp7PIolXl8WVBMRYmRKRnxgb9zWAK21bgcqU7t1WnCJvebBUfy7UfiPAW1oAxoI66oKleXRafbEBbfcB5GJEGk+qne0XxbU+K2teK+4Ke2ykxLnYdKCqnmg6IUBbWgDGgjzIFnBkZyUGCJDjZ5liB2lEq8DB9vIr2wKfNotuKxD252UmDDPXAeJI8FgUhIvhVNcGlohRLgQ4mshxHYhxHdCiD/r4yOEEJuEED8KIV4WQoTq42H663z9eJbDtW7Wx/cIIU51GF+oj+ULIf7gqzcX2EBYkdatNsoz36PRIJiS4UHigpNGjUcLXW3FA52oAC7r0HYnVW9pY+mlpc1hGEM0Y6uSFhROcGdH2w6cLKWcAkwFFgohZgF/BR6QUmYDdcAyff4yoE5KORp4QJ+HEGI8WuvxCcBC4B9CCKMQwgg8CpwGjAcu1Od6RUADYXAoot2Hgs+5w+LZtf8grR1uNGtMGg1NB6DtoMf3CTZf5NcQHxnC+LQAtRV3pK4YwuMgwj3tbkpsODYJ1U1uFmcHVfNA0SMuDa3UsIdSQ/QfCZwMvKaPrwUW688X6a/Rj88TWrn5RcBLUsp2KWURkI/WrvxYIF9KWSil7ABe0ud6RUADYaBX7crq06m5wxKw2CQ797mRuNBVXObo8tNKqbWtmT0ygG3FHakrcTsQBoe0tB4FxJLHavWILR4YZ8WgwC0frb7z3AZUAh8CBUC9lNKunC8H0vXn6UAZgH68ATA7jnc7p6dxZ+tYLoTYLITYXFVV1euaa5s7SIgMCcyO1mbTdkweKg7sHAqIueGn7ZJ4HV3ug5KaFvbWtwY+7daOhxrnrk4Lnkq8pFWrSaxQOOCWoZVSWqWUU4EMtB3oOGfT9Edn2xXZh3Fn61gjpZwupZyenJzc65pPn5TGt7fND0wgrHE/WNv7vKNNig5jWGKke37ahBEgDEfdjvaLLv9sEAJhNhvUl3r0++mqd+DJjjZ1svZYvtmDxSkGAx6pDqSU9cCnwCwgXghhz6HMAPbpz8uBTAD9eBxQ6zje7Zyexr0mYA3yugItfdvRgrardcvQmkK1r8BHWUDsy4IaUmPDGZkUFfibN1XoH4Tuuw7MUaGYDB60tAFtRxuZBCVf9GGRioGMO6qDZCFEvP48AjgF2AV8ApyvT1sKvKU/X6+/Rj/+sdTSntYDF+iqhBFANvA18A2QrasYQtECZut98eYChgflEXsiNzOeioNt7G9odT05Kfuoch10tRUfbQ5Od9g+aJy1ljZh7lfwAi0QmjUXipWhVRyOOzvaNOATIUQemlH8UEr5NnATcIMQIh/NB/tPff4/AbM+fgPwBwAp5XfAK8D3wPvANbpLwgJcC3yAZsBf0ecePdQWaZ1v4zJdz+2B3GFaNNytXa1dS2tzU3oUZHZXNFLb3BH4+rN27MkK8VkenZYS52F2GMDw46Ch9NA9FQrcaDcupcwDcp2MF6L5a7uPtwFLerjWncCdTsbfBd51Y739k7piiMvotfOtK8alxRJqMrC1tI7TJ6X1PjlpNFha4eBerRpVP6dLPxuMRAXQd7TC43+r1NhwfjjQ6Nm9suZqjyVfeOSqUAxsBnRmWMDwsDyiM0JNBialx/GtuztaOGoCYl8W1DAyKYq0uIjgLKCuGGLTweRZYDQlNtz9Uol2ksdBRKJyHygOQxlaX+CFhtaR3Mx4duxtoMPiwiVgl3jV9H8ZUVunlU2FNcHbzYLe+dbz3WVqXDhN7Raa3G1pA1pt2uFzoPgzj++nGLgoQ+st9qpQXigO7OQOS6DDYmPXfhdZX9EpWrpvPy4u026x8txXJcz7239p7rByyriU4C3GQw2tna6kBU8CYgBZx2nGvb7M9VzFoEAZWm+xR7S9dB2AB4kLQug1D/qfoW3rtLJuYzE/ufdTbn1zJ0Niw3jm8hn8ZMyQ4Cyos03TOffB0PZJSwuaoQUl81J0EeBeIgOQrqpdWV5fKi0unJTYMLaW1XOZq8lJ2VC6yet7+oq2Tisvfl3KY/8t4MDBdqYPT+Ce8ydz3Oik4Ei67NSXao8epN/a6VN2GMCQCRAeD8Wfw5QLPL6vYuChDK23+CBZwY4QgtzMBPclXjteg85WCAlSkAmtdfjzm0p4/H+FVDW2c+yIRB746VRmjwqSZrY7fawTDJASqwXPPJZ4GQwwfK7a0Sq6UIbWW+qKtGLc4b6pSHXM8Hje/66C6qb23tOHk0YDUitikjLBJ/f2BCkla78s5pFP8qlu6mD2SDMPX5jLrJFBDHo5o6sOrec72shQEzHhJio9NbSgybz2vKO1N3LRdVcx8FE+Wm/xkeLAjj1xYZurXa29rU0QAmJSSu54Zxd/+vf35KTE8MqvZvPi8ln9z8iCtqM1hWsBxD7gcUsbO8N1Pa2SeSlQhtZ76op94jawM3FoHCaDYGuZi4BYkPqHSSlZ9fb3/PPzIi6bk8Xzv5zJV9l6JgAAGuBJREFUsSMSA7oGj/CiTjBoftoKT7W0AKmTICwOSj7v0339wYMf/ciWEjc7eSh8ijK03mDthIZynygO7ESEGhmXFuvaTxsapYnwA1jzQErJn9Z/x9NfFLPsuBGsPGt8//DD9oaHdWi7kxIb7lm9AzsGIwyfrQXE+gH5lU088NEPbCqqCfZSBiXK0HpDfalWf9SHrgPQZF7by+qx2ly0IA+gxMtmk9z21k7Wbixh+QkjufWMcf3fyErZZw2tndTYcKqa2l3/LpwxfK5Wk6Kxos/39xUvbColxChYMq3/p2wPRJSh9QYfKg4cyR0WT3OH1XWefZJeXEb2wQh4gM0mueXNnTz3VSlXnjiKm08b2/+NLEBrHXQ0emVoU+LCsdok1U19cB/0Ez1tW6eV178t59QJqSTHBKA+s+IIlKH1Bh+UR3RGbqablbzMo7XMtOZqn97fEZtNcvO/dvDi16Vcc9Ioblo45ugwsnDo9+NFcZcU3TB5rKUFrRB4aEzQ3Qfv5O2nobWTi2YOC+o6BjPK0HpDbZEe0U716WWHmyNJiAxxnSHm5+IyVpvk96/n8fLmMq47eTQ3LjiKjCx4paG1Y09a8Dg7DMBogmGzgq48eOHrUkYmRTG7P6pCBgnK0HpDXbEWaDH49p9RCEHusAS2lrnY0Sb5T+JltUl+9+p2XttSzvWnZHPD0WZkwaEObd93tKl9TcO1k3UcVO+Bpso+r8Ebdu0/yJaSOi6aOezo+/0NIJSh9QYvGjK6IjcznvzKJhpaO3ueFJcJxjCf72gtVhu/fWUb/9q6l9/Oz+H6U3J8ev2AUVestZYJi+7zJczRYRg9bWnjSJD9tC9sKiXUZOD8aRlBub9CQ2WG9RUpNddB1vF+ubw9cWF7WT0n5PTQiNJgBPMor8sltnZYya9sYnfFQfZUNPJNSR3by+r53aljuOak0V5dO6h4qTgAMOotbSoa+hAMA0ibAiFRmvtgwjlercVTmtstvLF1L2dOSiM+MjSg91YcjjK0faW5CjqbfS7tsjMlMw4htIBYj4YWNENbuduta9pskqKaZvZUNLK7opEfKhrZc6CR4prmLuFCmMlAdko0ty+eyM9nHeUdAupLIH2a15fRCoD3cUdrDNH9tIEPiP17+z6a2i1cPEsFwYKNMrR9xYflEZ0REx5CzpAYNzLEsmHPe1ryRC+tdJrbLfzq2S18nq8pFAwCssxRjE2NYdHUoYxJiWFMagzDzVEYDQPAl2e1aPVgJ57n9aVSY8PJr2rq+wWy5sKGVdBcA1GBC0g9v6mUMSkxHKN/O1IED5eGVgiRCawDUgEbsEZK+aAQIhF4GcgCioGfSinrhOZxfxA4HWgBLpNSfqtfaylwq37pO6SUa/XxacAzQARa77Bf651z+y9d5RH9Y2hB09O+t7MCKWXPgYykbLBZtMBPkvOv+QfbOrn86W/YVlbPzaeNZe7oJEYPiSY8xOi3tQedg3u1ZBIvAmF2UmLD+KLACwmd3b1U8gWMP9vr9bhDXnk9O/Y2sGrRBBUE6we4EwyzAL+VUo4DZgHXCCHGo3W33SClzAY26K8BTkNrJZ4NLAdWA+iGeSUwE62p40ohhP2jdrU+137eQu/fmp+pK0Jr+Oe/r2W5w+JpaO3k0x+qep7kQuJV19zBxU9sIq+8nkcuzOVXJ45iYnrcwDay4BNpl52UuHAa2yy0dHjQ0saRobkQEhlQ98ELm0qJCDGyODc9YPdU9IxLQyul3G/fkUopG9FagqcDi4C1+rS1wGL9+SJgndT4CogXQqQBp6K1Kq+VUtYBHwIL9WOxUsqN+i52ncO1+i91xVr5u5Bwv93ilHEpjEiKYtkz33D/hz9gsTrpJWYepT3WHFnzoLKxjQvWfMWeA42s+fl0TnPVXXcg4UND2+eWNnaMIZB5bMCUBwfbOlm/fR9nTxlKbHjfOzMrfIdH8i4hRBZa6/FNQIqUcj9oxhiw9ypJBxybJZXrY72NlzsZ79/UFvnVbQCatOjtFcdx7jEZPLThRy5Y8xXldS2HT4pM1OrhdtPS7m9o5YLHv6K0toWnL5vBSWOD1EomWNSXgDBqhXe8pMvQ9jUgBprM68B30FLr9Xpc8dbWvbR0WFUQrB/htqEVQkQDrwPXSyl76x7ozCEk+zDubA3LhRCbhRCbq6p6+TodCOp8W4e2J6LCTNy3ZAoPXjCV3RWNnP7gZ7y3Y//hk8zZh+1oy2pbWPLYRqoa23l22bHMHZ3k93X2O2qLID5Ty87ykpS+trRxZPhxgISSL71eT29IKXl+UykT02OZnBHv13sp3MctQyuECEEzss9LKf+lDx/Qv/ajP9pTX8oBxxJBGcA+F+MZTsaPQEq5Rko5XUo5PTm5F8mTv+logaYDkJgVsFsumprOu9cdz4jkaK56/ltu/tcOWjus2sGk0V072oKqJpY8tpGmdgvPXzGT6Vn9uFasv/j/9u47vIoqb+D495cEQgnNJCA1IIRmA4wkSFDK0hQF28oiRcFFXdu66qNge23L+r52EQuKqIAgKiq6UgQR6UmQEooQeggkQOgl9bx/nIleMO1e72Tm3pzP88yTy7lzJ+c45ndnzpzzO6ePQNoP0DjOL4drUq86tcLD+Km0vvKyNO6kp2vb3H2wevcRNu8/zq3xAT40L8iUGWitUQQfAJuUUq94vPUNMMJ6PQL42qN8uGgJwFGra2Eu0EdE6lkPwfoAc633jotIgvW7hnscy51sytpVlmaRNfj8ri7c3b0l05N2c+34JXpp8uh2cDKL7RtXccu7y8kvLGT66ITKe0WT9D7kHIOu9/vlcOFhodx4WRP+u36fb1m8AMLCdT+tzQ/Epq7cRUR4GNddapbPcZPyXNF2BYYBPUVkjbVdDfwH6C0iW4He1r9BD8/aDqQBE4F/ACilsoHngCRre9YqA7gbeN/6zDbgez+0zT6H7R/aVZIqoSE82q8tn4yM5+jpPAa+tZQZeYnkh9ch87OHCBNhxp1daHu+f9YwCzi5p2DFBGjVW8/K8pOhCTHkFShmJO0pe+eSxCTC/vU6faMNjpzK5dt1+xjUsRE1w80QeTcp82wopZZQfD8qQK9i9lfAPSUcaxIwqZjyZOCisuriGjZPViiPxNgo5jzQjYdnruXR7/fya9ggngr7iNl9TxAd7fvc/oC3+iM4dQi6PeTXw7aqH8EVLSOZtlLn5PVpUkfzroCC3SugTX+/1g/gi9V7yc0vZEhn023gNiapjC+yd+j1oKo7O+MmMiKcSbddzlMD2rO9+S3knxdL9LJnIT/X0Xo5Jj8Xlr6hVzaI6eL3ww9LiGHvkdMs3OxjJq7GcToJkA3dB/oh2C46NatL+0aV9G7GxUyg9cXhHTqZtAtm3IgIIxNbMPmOroT1HwfZ22DVe05XyxnrpsPxDOj2L1sO37t9AxrUDueTFbt8O0CVatDkclsC7cod2Ww/cJIh5iGYK5lA6wsb0yP+KbG9dd/kT/9r66oLrlSQD0tehYYdoOUferT8Iiw0hCGdY1i85QA7D5707SDNu8L+dXplDD+aunI3tauFMeCSSjQpJYCYQOutwgKdV8CBB2Hl0vffkHsCFj7vdE0q1savIHu77pu18U5jcOemhIUIU1f6eFXbPBFUoe6n9ZODJ3KYk7qPmy5rGvxTqwOUCbTeOrYXCvMqZLKCT6JbQ+e/64dC+1Odrk3FUAp+fgWi2kDbAbb+qga1q9H3wvP5LDmdM3kF3h+gyeUQWtWv3Qczk9PJK1AMiTcr3LqVCbTeKsra5caugyJXPQrV6sCcx2xfIdcVtsyBrA26b9bPywoVZ2hCDEdP5zF7bbHzakpXpbp+KOaHiQsbMo5y/6e/8NK8X+lyQSSt6tf608c07GECrbccmqzglRrnQY/HYefPsPk7p2tjL6Vg8Us6i5ofcs+WR8IF59GqfgRTfH0o1jwRMtbA0fSy9z2HUopl2w4yfNIqrnljCQs2ZTIqsQVvDunoW12MCmECrbcO74CQML8kK7HVZbdDdFuY9wTk+zibKRDs/Bn2JkPXB0pNfO5PIsKwhBjWph9lbVkLaBan0zBd1wXPlvsjBYWK79fvY9BbSxkycSUbM47ySN82LBvTi7FXtyMqItz7ehgVxgRab2Xv0FdPfkhWYqvQMP1g7PAOWPG207Wxz+KXIKIBdBhaob/2+k6NqVE11LehXnWbQZd7YN0MSE8pddec/AI+XbWbv7zyE3dPXc2R03k8P+giljzak3t6tKJOdZMGMRCYQOutwzvd3W3gqVUvaN1PByOHlru2VXoy7PgJutxra17g4tSuVoXrOzZm9toMDp/0YYJI4oNQMxrmji2xH33BpkwSX/yRMV+up2Z4KOOHdGThQ90ZmhBjRhcEGBNovVVB6RH9ps8LkH/Gq9vUgPHzK1CtLsTd7sivH5oQQ05+IZ+neN/XSngt6PkE7FkBG/+YQ+lMXgFjZ62nbvUqTBkVz+x7ExlwSaPgWM+tEjKB1hunsvVAczePODhXVCuIvxN+mQL71jpdG//J3AC/fgcJd+ug5YB2DWtzefN6TFm5i8JCH0Z3dBwG9S+E+U9B3tm5bmcm7yHzWA7PDLyQxNgos+5XgDOB1huBMOKgOFc+okcizBkTPMO9lrwKVWpC59GOVmNoQgy7Dp3i5zQfZuKFhELfF/RqEKve/a04J7+ACYu2ERdTjy4XVNyquYZ9TKD1xm/pEZs7Wg2vVa+rh3vtWlrsbWrAyd4OqV/A5SP1F4iD+l10PlERVflkuY9DvVr28OhH14nFv0jZy76jZ7i/V6y5kg0SJtB6IztAAy1ApxH6NnXuWNvyoVaYJa9BSBX9EMxh4WGh3HJ5UxZuzvzjem7l1fs5yDsFi8aRV1DIhEVpdGhal26xlXAJoiBlAq03Du+EmvUhPADzvYaGwcA39RI839wXuF0IxzJgzTToOBRqne90bQB+y5g1beVu3w4Q3RriRkHKhyxY/BPph0/zgLmaDSom0Hrj8M7AvJot0vgy6PU0bJoNyX/Iv+5+SsG8J3VSFj8tU+MPjetWp1e7BsxI2kNOvg/5DwC6P4YKr0W9Jc9yceM6dG/j4Jp4ht+ZQOuN7B2BNeKgOF3u1WkE54zRT+4DyYoJkPo5dB/jui+8YQkxHDqZy5zU/b4doMZ5rG85mviC1fxP+wxzNRtkTKAtr/wcnbkr0EYcnCskBK5/Vz8gm3m7XmMrEGz/SV/Nth3g92Vq/CGxVRTNI2v4/FCsoFDx8M549oY0pNPml3V+XSNomEBbXkd2A8p1V1I+iYjWwfbgFp3hy+2O7IbPb4fIVnD9OxWSoctbISHC0IQYkncdZuHmTK8//+26DLYcyiWz81jkwGad5tIIGuVZbnySiGSJSKpH2XkiMl9Etlo/61nlIiJviEiaiKwTkU4enxlh7b9VREZ4lF8mIuutz7whbr1nCoT0iN5o2QMS/6n/oFO/cLo2Jcs9BdNv1Vd4g6c5NjmhPAZ3bkb7hrW5e8pqlm87VO7PFRYq3lyYRusGEXToPVSvlvvjv/2+CoPhnPJcGkwG+p1T9hiwQCkVCyyw/g3QH4i1ttHA26ADM/A0EA90Bp4uCs7WPqM9Pnfu73IHB5cYt02Px3Ui6tn//H0yhpsoBbMf0Et03zhRz3JzsYjwMD4Z1Zmm59Vg1EdJrN5dvmF036fuJy3rBPf2jCUkNERPYjh1CH5+2eYaGxWlzECrlFoMZJ9TPBAourf5CBjkUf6x0lYAdUWkIdAXmK+UylZKHQbmA/2s92orpZZby5R/7HEsdzm8E6rUgIj6TtfEf0KrwI0fAAKfj4KCPKdrdLYVE2D9Z/oLoXVfp2tTLpER4Uy9I57oWuHcNmkVGzJKvyrVV7NbuSC6JtdcbK331agDdBiis64VfQEWFurEQHtTYOM3sHwCzBkLM4bBxF7wxR2wJ8nexhk+8zXXXwOl1D4ApdQ+ESmKPo2BPR77pVtlpZWnF1NeLBEZjb76pVmzZj5W3Qf5ubBrme6fdWnPhs/qxcB1r8PM2/Q6Y72fcbpGmssffpWmQe1qTL0jnr++s5xhH6ziszsTSlz9YP6mTDbvP86rt1x6dsKYnk/Ahlkw+Vo9VfdYBhSck1e4Sg2o0wRqNYQt82D9TL16Q8Ld0H5gheXnNcrm76SqxUUh5UN5sZRS7wHvAcTFxVXMiPvCApg1Gvat0Q+QgtGF18P2RbD0NWhxpU6v6KQju3Xgd/HDr7I0qVeDqX9P4OZ3ljNk4kpm3tWFmMiaZ+2jlOKNBVtpHlmDay9pdPYBajfS+YTXfqqTzLe7Fuo01YG1TmP9unq937/4c07ofVe8DV+M0l9Sne/QCeAdnqZs+D7qINO67cf6WZTsNB3wXCGuCZBRRnmTYsrdQSn47l/6yqLPC3DpYKdrZJ++4yC6Hcy6q+TctXmnYedSPS9/6s3wf61g2mB9O+svRQ+/Cgtc//CrLC2iajL1jnjyCgoZMnElGUdOn/X+ws1ZbMg4xj09WhEWWsyfYtztMGoe3Pwh9HkO4kdD26uh4aU6eHreXYVH6EU5702GITMhuo1OjflKO/jmfsjaZHNrjdL4Gmi/AYpGDowAvvYoH26NPkgAjlpdDHOBPiJSz3oI1geYa713XEQSrNEGwz2O5bwFz0DKZOj2MFzh/Lx6W1WtATdNgpxjMOtO3Sd48iBs+lYvh/P+X2BcU5h8NSx8Ti+53uIq2L0cJvaEKTfpRNx/RoA9/CqPNufX4uOR8Rw7ncet768k67hOh1h0NdukXnUGdfTjskghIdC6Dwz/Cv6xQl8crJsBExLg44GwY7H/fpdRbqLKmPMuIp8C3YEoIBM9euAr4DOgGbAbuFkplW0Fy/HokQOngNuVUsnWcUYCY63DvqCU+tAqj0OPbKgOfA/cp8qqFLrrIDn5T/5hl2bJa/DD03oO+jUvB1/fbEmSJ8G3D0KtRnDcurkIrQqNOkGzeGjWBZrG/347euYYJE2EZePhdDa07AlXPab39dbyt3TSmx5PwFWP+K9NLpCyK5thH6yiab0aTB+dwNr0I9z2YRLjbriYv3W2+XnDqWx9wbBqor5g6HJPmR8RkRSlVJy9Fas8ygy0bmVroE2ZrK+sLroJbpgYkH2EPlNKB7vsHb8H1oYdyl4qJucEJL0Py96EUwf11W73xyDmiuL3z8+FQ2mQtVFPBc7aBFvnQZv+8NdPgvK/+bK0g9w2OYk2DWoRInDgeA6LHulB1bAKamtBnu6SKceyPybQ+pcJtOfaMEtPTY3trfsIzZNb7+Se1FfFS1+HkwegeTe9Qm1h/u8BNWsjHNwKhdZwspAwiIyFJnHQb1xA98uWZeHmTO78JIW8AsVzgy5iWEKM01Uqlgm0/mUCrae0BTDtFv0HP/RL3W9p+Cb3lL4zWPqaTs1YpE4zaNAe6rfT+XHrt4OoWAirPMtlz9+YybfrMnjxxktcu8iiCbT+ZQJtkT2r9MOCyJYw4luddMX48/JO6y+wiPoQ3Raq1Xa6RkY5mEDrX/4eRxuY9qfC1Jv0wO+hX5og609VqkO7AU7XwjAcFXxPHLx1YAtMuUEv9Df8q+CaYmsYhitU3kCbdwZ+HAfvJOoHNcO/groVOK3XMIxKo3J2HWz9Af77sM7IdfHN0Od516w/ZRhG8KlcgfboXpg7Ri+5HRkLw7+GC7o7XSvDMIJc5Qi0BXmw8h3dVaAKoOeTcMV9lWpIkWEYzgn+QLtruU4Mk7URWveD/i8Gx3I0hmEEjOANtCcPwvynYc0UnVJu8DRoc3XlyVlgGIZrBG+g3TIX1k2HxAfhykegas2yP2MYhmGD4A20l/4NmiXomV6GYRgOCt5xtCEhJsgahuEKwRtoDcMwXMIEWsMwDJuZQGsYhmEzE2gNwzBsZgKtYRiGzUygNQzDsJkJtIZhGDYzgdYwDMNmAbtmmIgcAHaVsVsUcLACqmOnYGgDBEc7KlMbYpRS0XZXprII2EBbHiKSHOgLzAVDGyA42mHaYPjKdB0YhmHYzARawzAMmwV7oH3P6Qr4QTC0AYKjHaYNhk+Cuo/WMAzDDYL9itYwDMNxJtAahmHYLOACrYhMEpEsEUn1KLtURJaLyHoRmS0itT3eGyMiaSLyq4j09SjvZ5WlichjAdqGndb+a0Qk2a1tEJFIEflRRE6IyPhzjnOZtX+aiLwhUnGLuvmxDYusc7PG2uq7tA29RSTFKk8RkZ4en3HsPFQKSqmA2oArgU5AqkdZEnCV9Xok8Jz1uj2wFggHWgDbgFBr2wZcAFS19mkfSG2w3tsJRAXAeagJJAJ3AePPOc4qoAsgwPdA/wBswyIgLgDOQ0egkfX6ImCvG85DZdgC7opWKbUYyD6nuA2w2Ho9H7jRej0QmK6UylFK7QDSgM7WlqaU2q6UygWmW/tWCD+1wVHetEEpdVIptQQ447mziDQEaiulliv91/4xMMjWinvwRxuc5mUbflFKZVjlG4BqIhLu9HmoDAIu0JYgFbjOen0z0NR63RjY47FfulVWUrmTvG0DgALmWbeBoyuklqUrqQ0laYxuTxE3n4eyfGh1Gzzpgtvu8rThRuAXpVQO7jwPQSVYAu1I4B4RSQFqAblWeXH/w6tSyp3kbRsAuiqlOgH9rc9eaX81S1VSG0oSSOehNLcqpS4GulnbMBvrVx6ltkFELgReBO4sKirmGE6fh6ASFMuNK6U2A30ARKQ1cI31Vjpnf5s3AYpunUoqd4QvbSi6DVRKZYnILHSXwmIcUkobSpKObk8RN5+H0j6z1/p5XESmoc/Dx3bWs4z6lNgGEWkCzAKGK6W2WcWuOw/BJiiuaIue8opICPAE8I711jfAYKsfqgUQi+70TwJiRaSFiFQFBlv7OsbbNohITRGpZX2mJvoPK/WPR644pbShWEqpfcBxEUmwbreHA1/bXtFSeNsGEQkTkSjrdRVgAC49DyJSF/gOGKOUWlq0vxvPQ9Bx+mmctxvwKbAPyEN/E48CHgC2WNt/sGa8Wfs/jn5S/yseT1KBq639twGPB1ob0CMm1lrbhgBow070Q5sT1v7trfI4dGDaBoz3/EwgtAE9GiEFWGedh9exRoW4rQ3ooHsSWOOx1Xf6PFSGzUzBNQzDsFlQdB0YhmG4mQm0hmEYNjOB1jAMw2Ym0BqGYdjMBFrDMAybmUBrGIZhMxNoDcMwbPb/av6oYeeEm40AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.axes([0.2, 0.1, 0.5, 0.8])\n", "plt.plot(year, hares, year, lynxes, year, carrots)\n", "plt.legend(('Hare', 'Lynx', 'Carrot'), loc=(1.05, 0.5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### So what was mean (median) population over time?" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([34080.95238095, 20166.66666667, 42400. ])" ] }, "execution_count": 102, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pop = data[:, 1:]\n", "pop.mean(axis=0)" ] }, { "cell_type": "code", "execution_count": 103, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([20897.90645809, 16254.59153691, 3322.50622558])" ] }, "execution_count": 103, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pop.std(axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Which animal had the highest population each year? " ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 0, 0, 1, 1, 2, 2, 2, 2, 2, 2, 0, 0, 0, 1, 2, 2, 2, 2, 2])" ] }, "execution_count": 104, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.argmax(pop, axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Diffusion using a random walk algorithm\n", "\n", "\n", "\n", "Let us consider a simple 1D random walk process: at each time step a walker jumps right or left with equal probability.\n", "\n", "\n", "We are interested in finding the typical distance from the origin of a random walker after t left or right jumps? We are going to simulate many “walkers” to find this law, and we are going to do so using array computing tricks: we are going to create a 2D array with the “stories” (each walker has a story) in one direction, and the time in the other:\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 1 : Let's initialize number of stories and max duration where we follow the walkers " ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [], "source": [ "n_stories = 1000 # number of walkers\n", "t_max = 200 # duration of time in which we follow the walkers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 2: Let's randomly choose all steps 1 or -1 of the walk " ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: DeprecationWarning: This function is deprecated. Please call randint(0, 1 + 1) instead\n", " steps = 2 * np.random.random_integers(0, 1, (n_stories, t_max)) - 1\n" ] } ], "source": [ "t = np.arange(t_max)\n", "steps = 2 * np.random.random_integers(0, 1, (n_stories, t_max)) - 1\n", "# use instead np.random.randint(0, 1 + 1) logic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 3: Lets verify if all steps are 1 or -1" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1, 1])" ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.unique(steps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 4: OK, we build the walks by summing steps along the time" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [], "source": [ "positions = np.cumsum(steps, axis=1) # this is axis = 1, dimension of time\n", "square_dist = np.square(positions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 5: Lets get the mean in the axis of the stories" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [], "source": [ "mean_square_dist = np.mean(square_dist, axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Step 6: Finally lets plot the results " ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$\\\\sqrt{\\\\langle (\\\\delta x)^2 \\\\rangle}$')" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAF1CAYAAAD1DaP0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUVfrH8c8zExKKSle6oa2IKAKxREUjICKiCCqiKy0IwoorYsUGrq6i7rp2XFzBCupPVGyAiMa2oy5VRcAaCUoziiCQOuf3x2QwhCSEtDsz+b5fr7yS3IzJc51Mvpxzz32OOecQERGR2OHzugARERGpXAp3ERGRGKNwFxERiTEKdxERkRijcBcREYkxCncREZEYE+d1AZWlSZMmLjEx0esyREREqs3SpUt/ds41LXo8ZsI9MTGRJUuWeF2GiIhItTGzH4o7rml5ERGRGKNwFxERiTEKdxERkRijcBcREYkxCncREZEYo3AXERGJMQp3ERGRGKNwFxERiTEKdxERkRijcBcREYkxCncREZEqFsgIcOcHdxLICFTLz4uZ3vIiIiKRJpAR4KmVTzFrxSzygnnE++NZPHwxya2Tq/TnKtxFREQqWeFQz8nPweGolQc55JCWnhb74W5mM4EBwGbnXJciX7sauAdo6pz72Yv6RERE9iWQESAtPY3GdRuzfMPy3aHuz3ec8S0MXwkp6XD41bVISUyp8no8D3fgCeAh4KnCB82sNXAasM6DmkRERMokkBGg91O9yc7LJkgQc3D0Bhj2GVz0ORyyA36uC5/36sLCQXdzTBWP2iECwt05976ZJRbzpX8B1wLzqrUgERGRMgpkBJiaNpXs/GyabQvy589Co/QuWyDbD28cZmwafDrdRk7m1PYnV1tdnod7cczsbOBH59xKMyvtcWOBsQBt2rSppupERKSmKjr9/sInMxmwKpcFKx29vwvdgvZRa7jsbD+1LxzOeSeOYXA1jNSLirhwN7O6wI1A33091jk3A5gBkJSU5Kq4NBERqWGKu5ael5vDST84RqyEe1bBAbnwXQN49px21Boxiu8b+7k4MaXKF82VJuLCHWgPtAXCo/ZWwDIzO9Y5t9HTykREpEYovNo9Nz+XIEHa/wKTV4am3dtuhW3xMOdIeKorLG1Xm8UjnvE00AuLuHB3zn0OHBz+3MzSgSStlhcRkapS0mr3etmOi1bByBVw8joIAm+3gxt7wauHG7m1a5F6dCp3dx0eMcEOERDuZjYHSAGamNl6YIpz7nFvqxIRkZqguBG6Lxi6bW3ECjh3NdTLhbWN4YbexvPd4uh7ymhSmnfjyJ2ZpHg8/V4Sz8PdOXfhPr6eWE2liIhIDVFck5n2mTCiYNr90N9gawLMPtrHr0MG0iDldA7c9QvPRGiYF+V5uIuIiFSn8H3pWXlZHJjluLhg2v2kDMg3WNQebjjdT+Mho7jw2NSoCPOiFO4iIhIzwtfOw9Plha+lZ+7MpHHdxrz8xYucvDaLYSscg1dDnTxY3RReGnkcOy8YTEa9fCZEyQi9JAp3ERGJSoWDHNhjgxa/z0//Dv2Z/8383dfSD/sZhq+Ax1ZCq+3wSx14sruPrRecwynnXcXgNid4e0KVSOEuIiJRJzy1npOfg9/nx7Dd184B8vPzeWXtKxyYBcNWQepySF4PeQYLOsCkfsbOfr258bS/RfUIvSQKdxERiSqFW74GXZBgfhBgd7DjoOcPoUA//8vQavdVTeGa0+CZo2DzQT4S/AksjtFgB4W7iIhEkaKbtPjMR5wvDsNoujWXESuNv3xRmxYbd7AtAeZ09ZE5dCANTzmdRrt+4daCa++RegtbZVG4i4hIVNhjxE4QHz76te7FvXm9Ofi516ifFsAXDMLJPfj62lOYd4SfEw/vG9MhXhKFu4iIRKTiusaFF8cdscUYs8LHX9YsoVbm29CiBVw/GUaOhI4d6Qhc7fUJeEjhLiIiEaW4rnGGcUCWY9gqGL0Mjv/REYxz+M7uBaNHQ9++EKdIC9P/CRER8VRJfd0drtDiOLd7cdyXBxvpN08g8fKboWlTr8uPSAp3ERHxREkjdIejxbZQG9jU5dDxF/ZYHHfy+VfTOYbuSa8KCncREakWpY7QAX8+9P/aMWYZ9P8a/A6+6tKcd/96OkuTE2vs4rjyULiLiEiVKm2EDnDorzB6eWiU3nI7bDgAFp3XjUMun0y3nufzJ+BUb08h6ijcRUSkShS381qYP99x9loYsxT6fhc6tj65Cy+eeTStho6lX7ueHlUdGxTuIiJS6QrvvFY41NtnwiXLYORKaPY7/NKkHj9ecQGtr5xCmzZtaONhzbFE4S4iIpWi8DX1uV/OJTs/G4cjPg8Gr4axy4xTv3cE/T6+Of4wfr0klcOHXUkjv9/r0mOOwl1ERCqkpGvqh21xjF0Gw1ZCk52Q1boZ3H4ZvlGj+FOLFl6XHdMU7iIiUm5Fp99r58L5q2DMMkfPdZDnN347PQWuuJ7affqAz+d1yTWCwl1ERPZLcdPvXTaGbmEbthIaZMPXjeDG0+MYdMcrJHU/0+uSaxyFu4iIlFnhXdnic4NcsAo++h8c/yNk++GzEzuw8aKz+eLwxgxoeypJui/dEwp3EZEaLjwSD2+DWnhknrkzc4/3c7+cy6GbshjzP8fIFdAoC9Y0hkcvOowe1/6LY7qeAcBZHp9TTadwFxGpoQovhMsL5uH3+enfoT/zv5m/V7OZWvkwcA1cuwR6fw+5PnjpcJhxjBFon8DiEbM4RqP0iKFwFxGpYUpqLpOfn88ra1/Z47GttjrGLg11kGv+O6TXhxt7Gd+eczK9TriIPjszub1gxC+RQ+EuIlIDFJ5qn7hg4l7NZQrzBeH0b2D8klCPd3Pwxp/g0SR4q6NRq1ZtFp9zpwI9gincRURiWNGpdzMj6IK7g90wavlr0b9Df5aueJNhS3IZs9SRuBV2Nj6Ijy/uzg9D+vJDfRhYtzEn7czcfW1eIpfCXUQkBpU09e5zPvw+P4bh9/lJ7TqKCbuO5IiXPiD4UhBfruO3E7rDFddR95xzODE+nhM9PhfZfwp3EZEYUtpmLYaREJfAff3uY8em9Zzz6TbaXrsQ1vwbGjTAd9kEuPRS6nfq5OEZSGVQuIuIRLF97ZEOf0y9px6dyng7lqMe/gCeew527YLjjoNZs+CCC6BOHQ/PRCqTwl1EJAoUd+95OMyL2yMd/gj1sUeM4K8ZLej494Xw8aNQty4MGwbjxkG3bh6elVQVhbuISAQraVMWh9srzIsukpvUagiXLfPT6sF5sHkzdOwI998PI0ZA/fpenZJUA4W7iEgEKu3aefjjorey+cxHnPmZ5u/H8A+30/itORAMwoABMGECaOOWGkPhLiISYYrutFZUeMTuMx9xvjhSj07lmPqH0+KVxZz82ufU/fo1aNQIrroqNPXetq0HZyFe8jzczWwmMADY7JzrUnDsHkKtiXOAb4FRzrmt3lUpIlK1ittpreg0e+rRqXRr3m2Pfu/98tvS7eUAPHEzbNsG3bvDrFu0QK6G8zzcgSeAh4CnCh1bBEx2zuWZ2V3AZOA6D2oTEalSpV1TLzwyH951+B+NY/Lz4Y034KEXYdENUKsWDBkSmno/7jgw8/akxHOeh7tz7n0zSyxy7K1Cn34MnFedNYmIVIeSpt8dDh8++rTtw9SUqX+EemYmPP44PPII/PADtGwJt98Ol1wChxzi0VlIJPI83MsgFXje6yJERCpTICPA1LSpe0y/h/nMR4I/4Y9gX7kytMp99mzIzoaUFPjnP2HgQIiLhj/jUt0i+rfCzG4E8oBnS/j6WGAsQJs2baqxMhGR8iluGr7w9Hv4mnpK654kL90Ew1LgvfdC96aPGgWXXQZdunh9GhLhIjbczWwEoYV2vZ1zxW5d5JybAcwASEpKKn57IxGRCFBir/ei0+9bt8LMmTBiOHz/PbRpA3ffHZp6b9jQ47OQaBGR4W5m/QgtoDvFObfT63pERPZXWdvCJsQVTL/vahxaEPfEE7BjB/TsCffco6l3KRfPf2PMbA6QAjQxs/XAFEKr4xOARRZa9fmxc26cZ0WKiJRRaavfw3bf2tZ1FJdv60TncbfDm29CfDwMHQpXXBG6pU2knDwPd+fchcUcfrzaCxERqaDSVr/DH6E+rtMwrviqEe2uex1W/zu00n3q1FDDGa16l0rgebiLiMSCfa1+j/PFcVWL87nsU2h570vw66+h0flTT4XuUU9I8KhyiUUKdxGRCtjn6vdmR1Pnk2X0f/NrGi+YE/qPBg8OTb2feKIazkiVULiLiJTDvla/33riTRz/cQZMvReWLoUGDUK93i+7DA491OPqJdYp3EVE9kNpu7UZxsG58cxY05FDb7gI1q+Hww6D6dND+6fXq+dh5VKTKNxFRMqopAVzhtFhWxzTv+lEyuJv8e94ONRFbvp06N9f26xKtVO4i4gUUfge9cydmaQkpgDstWDOME74yc/9q9rQPZCO+VaHdmObNEm3somnFO4iUuMV13Cm6OI4w3Yfi3PGoLV+7ljZhA6rN0L9zND19Msvh9atvT4dEYW7iNRcZWk4E3RBcvNzAaiT40hdbkxeWpsWm3dBYm247z5ITYUDD/TqNET2onAXkZhXdJq9tJawxd2j3vp3P5d9HGT0knwa7XJs79YOHpoCgwapNaxEJP1WikjMKm1kXnSEHlb4HvVeOw6m06zXOOKdz7G8IGt7HsGPE//KkYPGenA2ImWncBeRmFGWzVrCHxdtCbt7u9UdP3P2pgYcMf0NePPR0Farl46DiRPp1L69J+clsr8U7iIStUpbCFfSyDx8vPAIfXjX4SS3OBbmzYO7HodPP4WmTeG222D8eGjc2IOzEyk/hbuIRKXwPefZednFhnmJI/OCa+7hW9ySm3aDp5+Gf4yEr76Cdu3gkUdg5EioU8ebkxOpIIW7iESdwpu0BAkCxS+E22Nk3jp5z2+ydSs8+ijcPxg2bgzdl/7886G+71okJ1FOv8EiElVKGrHvsVlLwQg9JTFl71D/8cfQ7Wv//jds3w59+8Izz0CvXtrERWKGwl1EokbREbsPH33a9eHczueWHOZhX34J99wDzz4L+fmhTnLXXAPdulXvSYhUA4W7iES8krZVTfAnMDVlasmBDvDRR3DXXfDaa6Fr6JdeGmoP27Zt9Z2ASDVTuItIRCtus5bwtqolBrtzsGAB3HEHfPhhaLX7lCkwYQI0aVLNZyBS/RTuIhKxCk/DF179nhBXwog9Px9efjkU6suXh/q8338/jB6t7ValRlG4i0hEKrpwrtTV77m5MHs2TJsGa9ZAx44wcyb8+c8QH+/dSYh4ROEuIhGn2IVzxU3D79oFs2bB3XfDDz9A166h29nOPRf8fu9OQMRjCncRiRhlXji3fXvoHvV//hM2bYLkZHj4YejfX7eziaBwF5EIUaaFc5mZ8MAD8OCD8Ouv0KcPPPccnHKKQl2kEIW7iHhunwvn4hLh6qtDo/UdO2DgQLjhBjj2WG8LF4lQCncR8VRpC+fGNDqN7nc+HVocl5sLQ4fC5MnQpYvXZYtENIW7iHimpIVzdyZeQvcnFsBTF4Sm20eOhOuuA225KlImCncRqTalbdHqMx+df63F06/GcfDLF4Y2bxk/Hq69Flq18rp0kaiicBeRKlPW/dY7/Ww8uOwQen+8CYt/By6/PBTqzZt7fAYi0UnhLiKVqnCgT1wwsdT91jttgZveh6FfOEj4BZs4MbSZS7NmHp6BSPRTuItIuRQO8cydmXuMzvOCeZgZQRcsdr/1zpvhlvfh/FWQGx/HxrEX0PJv98LBB3t1OiIxReEuIvuluEYz4VF54dG5z/nw+/zg2H1NvetmH4+taEP3/35PTu1abBh/AS2n/pOWTZt6fFYisUXhLiJlUjjUc/Jz9ppiL/w+fI/6ff3uI3NnJh3X76TLoy/T6f1VcOAWuOEGEq68kpaNG3tyLiKxTuEuIqUqLdTDwiP2vTZ32XYQ3DkF5s6Fgw6Cm2+GiROhUSMPzkSk5vA83M1sJjAA2Oyc61JwrBHwPJAIpANDnHO/elWjSE1UWqgbRi1/LVKPTqVb8267r7ln7swkJTGF5J2N4LpbQ61hDzgAbrklFOoNG3p4RiI1h+fhDjwBPAQ8VejY9cBi59w0M7u+4PPrPKhNpEYqrs877Bnqe227CvDtt3DzbfD001C7dqjxzNVXg6bfRaqV5+HunHvfzBKLHB4IpBR8/CSQhsJdpFqU1Oe91FBftw5uuw2eeCLUfGbixFCwa/W7iCc8D/cSHOKc2wDgnNtgZvoLIVINSuvzXmyo//QT3HEHPPZY6PNx40K931u0qP7iRWS3SA33MjGzscBYgDZt2nhcjUh0K6nP+x77qIdt3gzTpsH06ZCXB6mpcNNN0Lq1N8WLyB4iNdw3mVnzglF7c2BzcQ9yzs0AZgAkJSXtvYRXRMqkuBF7gj9h72DPzIR77gntp56dDcOHh1bAt23rXfEishef1wWU4FVgRMHHI4B5HtYiEtNKGrEvHr74j2Dftg2mTg2F+N13w6BBsHp1aCtWBbtIxPF85G5mcwgtnmtiZuuBKcA04AUzGw2sA873rkKR2FRcp7m9RuxZWfDoo/D3v8PPP8O558Lf/gadO3tdvoiUwvNwd85dWMKXeldrISI1REn3r+9xjb35MTBrVmi0vm4dnHZaaOFcUpK3xYtImXge7iJSPfbVlCYhLoGpp0wheclGuPGo0LT7MceEpt5769/aItFE4S4S48raaW7CjiM44vwr4dNPoVOnUMvYQYPAzMPqRaQ8FO4iMawsnebG27Ecdf8cWPRo6Fa2mTNh2LBQMxoRiUp69YrEmML7rM/9cm6JnebGHnAK3R55Cf4vNdQe9t57Yfz4UNtYEYlqCneRGFL0fvXidmu75JAz6PHY6zDz4lCQ33ILXHVVaNc2EYkJCneRGFH0fnUI7a8eXgV/W49rOPbZd+FfQ0Nd5S67DG68Uf3fRWKQwl0kBpQ2Yq9HPI+md6HtpAtD96pfdBHcfruaz4jEMIW7SJQrtsNcuz6ce/hgGi94n/4zP6DO9/dCSkqodazuVReJeQp3kShVWoe5fxx4Lkde8ST897+hbnKvvw79++u2NpEaQuEuEmVK6zA3vE4y/3i3Fo1vuhSaNYMZM2DUKN3WJlLD6BUvEiVKa0bTZAfc+oExbsnH+BJqw623wqRJcMABHlYsIl5RuItEuNJCvXYuTPrUz00f+amdnY9dckmoH3yzZt4VLCKeU7iLRKhS28Y6uGBNHA++W4cmm7fDWf3hrrvg8MM9rFhEIoXCXSQCldY29tiNfp7+oCkdV22Ao9rCnH9Br14eVisikUbhLuKhcKvYlMQUgFLbxrbeGcecZe1IXvwV1iQvtFguNRX8fg/PQEQikcJdxCPh0XlOfg5+nx/Ddt/SVrgJzQH5fp767mgGvPQF/rzv4Zpr4IYboH59r09BRCKUwl3EA3s0nnFBgvl/tIsNv/c547bNXbhq3iYSfvwfDB4Md98N7dt7WbqIRAGFu0g1K9oqNrypS+GR+7E/Gf9aYJyw7jPo2hWeeS7UYU5EpAwU7iLVqNhWsW37MDVlKgD/W/Iq/R5P409vfExOkwbwn3tg5EhdVxeR/aJwF6kmxY3YE/wJTE2ZSnKzJHjgAZJvfRiys+G664i/4QZtwyoi5aJwF6kGpY3Yk9f8Dn2PgjVr4Mwz4b77oEMHr0sWkSimcBepQqVt7nJHh0vpceU/Ye7c0CK5116DAQO8LllEYoDCXaSKFNeIxoePM1qeyiOr29Gm98WhB95+O1x1FdSu7WG1IhJLFO4iVaDwNPzuRjQOBn8Tx9MzV1N73WI4/3z4xz+gTRuPqxWRWKNwF6lEJU3Dd/rVz/MfNKPL0gzo3AAWP62WsSJSZRTuIpWgpE1eDsgxHvvsUIa8tR5fwla4916YMAFq1fK4YhGJZQp3kQoobee2gWvggfnQ5rfvYcQImDZNW7GKSLVQuIvsh/BGL43rNmb5huXFhvqhW+HBBcZZaxw7D2sLbzwFJ57oYdUiUtMo3EXKqGgTmvDmLmG18uHqT/xMec+I88XBPbdR94orNAUvItVO4S5SBkWb0AB7bMeakuHnmUUH0mLdr3DOOXD//VoFLyKeUbiLlKK41e+Ft2M9ZJefF5e044S318KhB8GrT8JZZ3ldtojUcAp3kRKU1ISmT7s+nNtpEC1fXETfGW9T6/dv4brr4OaboV49j6sWEVG4i+yh8IK5uV/O3bMJDUZCXALTWo2k21XT4YMPoGdPmD4djjjC48pFRP4Q0eFuZlcClwAO+BwY5ZzL8rYqiVUlLZgL77c+5ogRXPsRtOk3MjRCf/xxGDUKzLwuXURkD6WGu5nlO+c82UjazFoCfwU6O+d2mdkLwFDgCS/qkdiXlp5GTn7OHgvmwru3/ePAcznypgdg1SoYOjS0c9shh3hcsYhI8fY1cvd6SBIH1DGzXKAu8JPH9UiMCmQEWPfbOuJ8cbh8t7ttbOO8eJ545yCaPzkOWrbUzm0iEhX2Fe5uH1+vMs65H83sH8A6YBfwlnPuLa/qkdhSXDOavGAefp+fsT3G0q15Nxq+/RFnP7qAhI1zQy1j//53OPBAr0sXEdmniL3mbmYNgYFAW2Ar8H9mdrFz7plCjxkLjAVoo3uKpYxKbUYThE75DRl792J44YXQQrmX5sHxx3tbtIjIfvCV9z80s4ZmNtzMXjazVWb2upmNMbPKuhDZB/jeObfFOZcLvAScUPgBzrkZzrkk51xS06ZNK+nHSiwrtRmNg9SVPi4b8RC88grcdhssW6ZgF5GoU66Ru5m9BDQE3gCuc859ZWZtCI20nzazeOdcSgVrWwccb2Z1CU3L9waWVPB7Sg21r2Y0h27389o7zThiWQacdBw89hh06uR12SIi5VLeaflU59zWwgecc+uAB4EHzaxBRQtzzn1iZi8Cy4A8YDkwo6LfV2qeUpvRHD6Yli++xen/XkhcfiY88ABcdhn4yj2pJSLiuXKFe9Fg39+v78fPmQJMqYzvJTVT4Wn4os1o7vjTeHrc9CgsXAinnBK6b719e48rFhGpuP0OdzM7DRgCPOycW2FmY51zGlFLxCm6cC7cjCa16yiuXt2Q9qcNh/x8eOghGD9eo3URiRnlGbn/BRgF3GRmjYCjK7ckkYorunAu3Izmjo7j6DHlUXjr35CSEhqtt2vndbkiIpWqPOG+pWDa/WozmwYcU8k1iVRIcSP2BF88j/zUjfaXjYBgEB5+GMaN02hdRGJSecL9jfAHzrnrzezySqxHpEKKG7EPrX8SD7+URYP37wqN1mfOhLZtvS5VRKTK7He4O+fmFfn8QTOrB2Q55/IrrTKR/VDcrW4+8/HnVX5mLlhOXF6+RusiUmOU9z53H6FNXP4MJAE5QIKZbQHeBGY4576utCpFSlHcrW6Ndhkvvncwp368EY7vAU89BR07elypiEj1KO8Q5l2gPTAZaO6ca+2cOxjoCXwMTDOziyupRpESFXerW59v4bPpjlOWbAl1mfvgAwW7iNQo5W1i06egJewenHO/AHOBuWZWq0KViexD0YVzdfOMaYuNywNBdrVvg++5uZCU5HWZIiLVrrxNbHIBzOw+4Ern3F67xxUX/iKVIbyj27rf1u3ef73HT8bc1+ty6E874PLLqTNtGtSt63WpIiKeqOiucL8Dr5rZUOfcDjPrC0xxzp1YCbWJ7CU8Ws/Jz8Hv85Pg/Fz5QZAp7zpc0zqw8CXo29frMkVEPFWhcHfO3WRmFwFpZpYN7ACur5TKRIrY4zY3FyQx07FwYVPar97EzwN60+TJF6BRI6/LFBHxXIXC3cx6A2MIhXpzYLRzbm1lFCYSVtxtbhd9YTzyWpB68b/Ds8/S5MILwczrUkVEIkJFp+VvBG5xzn1gZkcCz5vZJOfcO5VQm9RwhUM9Jz8Hh6NeNjw0H0aucGzv0YW4F1+DxESvSxURiSgVnZbvVejjz83sDEKr5U+oaGFScxUX6gDdfoLnXoQOv8L6v46k1T8fg7iK/vtURCT2lLeJTSrwrHMu28wGAs2Az51z/y2Yqhcpl+Ia0lgQJn4C096GHQ3qsfq5uzhiyGUeVyoiErnKO+y5wjk308ymAqcCAWCwmR0IDAJ2VVJ9UoMU15Dm4N/hyXlGv68dv/TtSaPZL9OwcWOPKxURiWzlDfecgvf9geRwT3kzOxOYDgyuhNqkBiluJ7fTvvfx/Cu1OHBXPjz8LxqNH69FcyIiZVDe9rMZZvYEcDBQJ3zQOfcGoO22ZL+lpaftbkgTFzRmrWjL/Cfzqd88Ed//lsBf/qJgFxEpo/KO3EcC5wL/ItRqdgGwCujGH6N6kTIJZARY99s64nxxHLzN8fRL0Pvbb2HYMJg+HerV87pEEZGoUt72s9uAWQBmdj5wKaHA/xW4oLKKk9hWeFV8XjCPU9YZc1+uzUE78uCxh2H0aI3WRUTKobyr5S3cT74g6O8p7TEiRRVeFU/Qce1H8Pd34LeW9fG9sxi6dvW6RBGRqFXuLV/N7HIza1P4oJnFm1kvM3sSGFHx8iRWha+xN9zpeG0OTFsML3Xx883C5xTsIiIVVN5r7v2AVGCOmbUFthJaWOcD3gL+5ZxbUTklSqwJX2NP/tHHs8/nc8jvMGfcSSROvotj26j/kYhIRZX3mnsW8AjwSMG+7U2AXc65rZVZnMSG8Batjes2ZvmG5cxaMYuRn+aw+E3HtiYHsvbZ+7nwzFFelykiEjPKe83d55wLwu592zeU9HWp2Yrev147Fx6cD2OWwcIO8OV9V3Clgl1EpFKVd1p+kZn9DLwCvOGc22ZmdQlN1w8CuhC6LU5quML3r7f6DeY+D8f+BH/vCXeeVptFR/X3ukQRkZhT3mn53mbWGRgIvFEwNe+AhYSuty+rxBolShW+f/3kb4I893+O2vlw/oVxNLnoEhZ1HU5y62SvyxQRiTnl3lLLOfcl8CVwp5nVLrgOLwL8MR2fk5fNlZ8Y0xbCtjbNmHPbUCadMkShLiJShfYV7mZm5wELnXPbS3qQgl0KC28A49+VxVPzHBd9AWt7drRiZpUAABa2SURBVOaw1wNcetBBXpcnIhLzSg1355zPzFoBQ8zsIGA78JZzbl21VCdRJzxiP/jnLD6Y4zhqM0zpE0e//8wABbuISLXY57S8c2498DiAmdUDTi9oOZsL/Nc5t6RqS5RoER6xd0vPYu4cR508uHlSNwZc8bCm4UVEqtF+XXN3zu0AXoLQ7W7ACWY2ETBgLbDYOZdd6VVKxAuP2M9bmsWrrzrWHwRnjE7gEQW7iEi1q8iCuiDwYcEbZtYRGGVmzYF04MXSrtNL7AhkBLj1nSncvHAXkz+AdxPh4WtSeOSsOxTsIiIeKHe4F2Zm7YCzgDOBIJAJ+Cvh+zYA/kPovnkHpDrnAhX9vlJ5AhkBzn6sF4/9XxbnrIEZPeCas2qzQMEuIuKZcoe7mRlwG9AbWAfMA86t5Ba09wMLnHPnmVk8ULcSv7dUUCAjwMMvXseix7I4chNc0c9Yc2EfFpx6q4JdRMRDFZmWd2aWA/Stiun3gtX5JxPaJx7nXA6QU9k/R/ZfeB/2L199nBdm55KQBwMuNt47rDaLFewiIp6r6LT8bVW4Z3s7YAswy8y6AkuBKwoW9QFgZmOBsQBt2rQp9ptI5QmH+qwVszj7s2wWvgTrD4JTRxqtjzuNxSlTFewiIhHAqi6bK8bMkoCPgROdc5+Y2f3ANufczcU9PikpyS1Zorvyqkp4NXxW7i6uDMA/34IPW8OgobCjfh0WD1+sYBcRqWZmttQ5l1T0uM+LYspoPbDeOfdJwecvAt09rKfGCt+/npubxf3zQ8H+QmfoP6oW5508TsEuIhJhKmW1fFVwzm00swwzO8w5t5bQwr0vva6rpgmP2P07s3hxrmPgWrj3RB/fXDuGhd1GKNRFRCJQxIZ7gcuBZwtWyn8HaOPvahQesdf/LYt5zzp6bICHhnci+faZTFKoi4hErIgOd+fcCmCvawlS9cIj9labsvjoaUezHTD0ongm3T5To3URkQgXydfcxSPhEfvh67P4YKbjoGy4/objmHRnmoJdRCQKRPTIXapfeMR+/NdZvDvH8WsdGDgygemj/6VgFxGJEhq5yx7S0tPo/3k2859xZNSHG249melXvatgFxGJIgp32S2QEaD1Cwt4/oUgy1pA3zG1mTB4moJdRCTKaFpeAAis+y8LRp/CrW/nsaCjMf+OVF48brSCXUQkCmnkLuAcXD+ZW9/O45kj4ZwLjWaHtFewi4hEKYV7TRcMwl//SvKc95lxrJ9Rg3344hNISUzxujIRESknTcvXZPn5cOml8PjjfDLkJLjmYv626xdSElM0ahcRiWIK95oqLw9GjIDZs5mWEsdNnf9L/MKl6hMvIhIDNC1fE+XkwAUXwOzZzDy/AzeeGiSfIDn5OaSlp3ldnYiIVJBG7jVNVha/ntmLhu8EuKq/n/uO+I6gC+IzH/H+eF1rFxGJARq51yTZ2fzaPxTs4wbAvcfmEySIDx992vbRlLyISIzQyL2myMmB886j4bsBLj3LmNHDAWAYCXEJTE2ZqmAXEYkRGrnXBLm5MGQIvP46s//SkyePjcdvfuL98Vza41KN2EVEYoxG7rEuNxcuvBDmzWPSgFo8cMh/8eNnTPcxDO86XKEuIhKDNHKPZXl5cPHFMHcub084kweOCZLv8skP5tOmfhsFu4hIjFK4x6r8/NB97C+8wOK/nMF3I88m3v/HdLxWxYuIxC5Ny8ci52DcOJg9m5tPi+POQ94ifkEa9/W7j8ydmepAJyIS4xTusWjyZPjPf5h9dlvu6PEDQRdqUJO5M5PJPSd7XZ2IiFQxTcvHmB9unAB33cWMY3wM65auBjUiIjWQwj2GfHvPDRx6x8M8dwSMPyNI0Jwa1IiI1EAK9xixdsadJF53Jws7wPBBEPSpQY2ISE2la+5RKpARIC09jZTEFA76aAkd/nIDn7SEwUMgv5aPeF8cqUen6l52EZEaSOEehQIZAXo/1Zuc/BySNhiLZwX5qjEMuAiy4kPT8Bqti4jUXAr3KJSWnkZOfg6tfsnnladhS204/WL4rZ6PBL+m4UVEajpdc48ygYwA635bR5NsP28+Cwn5cMbFsOkgLZwTEZEQjdyjSHg6nuxs5s92dNzqo/9wP18fHCTBH68Ru4iIAAr3qBHICDA1bSrZeVk8+YrjlHSYd/P5/G3MxN0L6xTsIiICCveoEB6xZ+dlc8u7jos/hyl94ug3ZiLJrZMV6iIisgddc48C4QV0530RZMp7sLBnC/o9/p5CXUREiqVwj3DhBXTHbPDxxCvw0aE+6s+cTXKbE7wuTUREIpSm5SNYeDq+0dZsPp4dZGfDA4h/+TmO6XCK16WJiEgE08g9gqWlp2FZ2cydE6RBFrxy1yiO6Xam12WJiEiEi+hwNzO/mS03s9e9rqW6BTICrNv6A4++Acf9CKPPi6dznwu9LktERKJARIc7cAWw2usiqlt4Ot5mzGDY8iCvX5TExGlpWkAnIiJlErHhbmatgDOB/3hdS3VLS0/j6B+yue9Nx5sd4fOx5yjYRUSkzCI23IH7gGuBYEkPMLOxZrbEzJZs2bKl+iqrQoGMAL+sW8sLzwf58SAYc35tUtr18rosERGJIhEZ7mY2ANjsnFta2uOcczOcc0nOuaSmTZtWU3VVJ5AR4LQnetHvlidpvBNmTz2PF8e9o1G7iIjsl4gMd+BE4GwzSweeA3qZ2TPellS1wu1lb1iURe/vYMIAw9e9u4JdRET2W0SGu3NusnOulXMuERgKvOOcu9jjsqpMeAGdW7SI6z+Ax7vDnKTapCSmeF2aiIhEoYgM95omLT2N+r9l8+TLjjVNYN743tq6VUREyi3iO9Q559KANI/LqFIpbU6m+zxomAUDRyZw/+m3KdhFRKTcNHL3WCAjwI5/3MHpXwV5b8JZ3H/tuwp2ERGpkIgfuceyQEaAq6adStq/s3ntcB9NrrxewS4iIhWmkbtHAhkB7njrZma+kM3mujD6bEj74T2vyxIRkRigkbsHwqvjb39jF51+hr7Djd8PTNDqeBERqRQauXsgLT2NY7/LZmIApieB9TlNq+NFRKTSaOTugV4HH8eQVxzfN4RbzqjNqylTFewiIlJpNHKvZoGMAHGTb6T9L47ArWN4dYzay4qISOXSyL0aBTIC/G1KCvNfzuGhZD89zhmlYBcRkUqnkXs1CWQEuHPhzTzycg5fNYLJvRxp6WlelyUiIjFII/dqEF4df/OCXbTdCqeOMvLraHW8iIhUDY3cq0FaehrtN2Rz9X/hya4Qf6pWx4uISNXRyL2KBTICrNv6A4+84dieADf3r83zWh0vIiJVSCP3KhSejs+aNYOe6Y43R5/M8+O1Ol5ERKqWwr0KpaWnccD2bO5e6Ai0hozz+irYRUSkymlavgqlJKZwyGKjYRZcMTCB+9v18rokERGpATRyryKBjABfvj2HUUuDLDvvRG3lKiIi1UYj9yoQvtY+b9Yufq0N3HSzgl1ERKqNRu5VIC09jZS12Zz2Lfz9ZGPx1mVelyQiIjWIRu5VIKV1T85c5PiuITyenMB8NasREZFqpJF7JQtkBMiccR9HbXR89tehzE/VrW8iIlK9NHKvRIGMAAP+04uVD2fxv1bGIamXK9hFRKTaaeReidLS0xj/UTattsE1fY20H97zuiQREamBFO6VJJAR4Of1X3Hth455neDTdtoYRkREvKFwrwThW98aPfokB2XDZxOGaGMYERHxjMK9EqSlp1H392wmfOJ4sTPEHXW0gl1ERDyjcK+gQEaAdb+t48qPjfrZcFcvTceLiIi3tFq+AsLT8Qk7svn+4yDLk9vywLXPatQuIiKe0si9AtLS08jJz2H0kiANsmBFan8Fu4iIeE7hXgEpiSnUc7WY+AmktfXR6fQ/e12SiIiIwr28AhkB0tLTeIWhtNoGh0y9R6N2ERGJCLrmXg7ha+05edksn+7YeVg7Dh92pddliYiIABq5l0v4WnvP74McucnxzqBuYOZ1WSIiIoDCvVxSElOI98fzlyVGZh1omjrB65JERER2i9hwN7PWZvauma02s1VmdoXXNYUlt07m/b7PMXiNkT38Io7rmOJ1SSIiIrtF8jX3POAq59wyMzsQWGpmi5xzX3pdGEDS/BWQH6TF1bd6XYqIiMgeIjbcnXMbgA0FH283s9VAS8DTcA9kBHjvm8VcOf0hEk4/HTp08LIcERGRvUTstHxhZpYIdAM+KXJ8rJktMbMlW7ZsqfI6wqvkP3p8Cgkbt7D23JQq/5kiIiL7K+LD3cwOAOYCE51z2wp/zTk3wzmX5JxLatq0aZXXEl4lf9HKIJl14JX2uVX+M0VERPZXRIe7mdUiFOzPOude8rqelMQUGufW4pw1MPdIPyd37ON1SSIiInuJ2HA3MwMeB1Y75+71uh4IrZL/qMEk6uRB8uRH1JFOREQiUsSGO3AiMAzoZWYrCt76e11UhzcC0KEDRw4c43UpIiIixYrk1fIfAhHT9i2QEWDJ/+YxIS0NmzpVHelERCRiRWy4R5LwKvlRgSzMOVac1JGjvS5KRESkBJE8LR8xwqvkz1rj+KoxzI/73uuSRERESqRwL4PwKvle6fBmpzhS2p7qdUkiIiIlUriXQXLrZN5vcRPx+dBn4gNaJS8iIhFN19z3IZARIC09jUsWvQ+NG9Nl0FivSxIRESmVwr0U4YV0wZxsxs0PsnlAfw72+70uS0REpFSali9FeCFdj/VBGmbBh0c18LokERGRfVK4lyIlMYV4fzwnrwvd0976rD97XJGIiMi+KdxLkdw6mcXDFzN6x5/Y2eFQjjna8wZ5IiIi+6Rr7vuQ3PI4WL0JzjvP61JERETKRCP3ffnyS9i6FU46yetKREREykThvi8ffhh6r3AXEZEooXDflw8/hGbNoF07rysREREpE11zL0a4cU1KYgrJH34IPXtqFzgREYkaCvciwo1rcvJzaPt7Lb7+IQsmTfK6LBERkTJTuBcRblyT7/Jpv8mFDnbt6m1RIiIi+0HX3IsIN67xm59WOwpazbZo4W1RIiIi+0Ej9yLCjWvS0tO4aFcGMB2aN/e6LBERkTJTuBcjuXVyaFvXuRPhgANCbyIiIlFC0/Kl2bBBo3YREYk6CvfSKNxFRCQKKdxLo3AXEZEopHAvjcJdRESikMK9JNu3w44dCncREYk6CveSbNgQeq9wFxGRKKNwL4nCXUREopTCvSQKdxERiVIK95Io3EVEJEop3EuyYQMkJEDDhl5XIiIisl8U7iXZsAGaNdM+7iIiEnUU7iXZuFFT8iIiEpUU7iVRAxsREYlSCveSKNxFRCRKRXS4m1k/M1trZt+Y2fXV9oOzs+GXXxTuIiISlSI23M3MDzwMnAF0Bi40s87V8sM3bgy9V7iLiEgUithwB44FvnHOfeecywGeAwZWy0/WPe4iIhLFIjncWwIZhT5fX3BsNzMba2ZLzGzJli1bKu8n16oFp50GbdtW3vcUERGpJpEc7sXdYO72+MS5Gc65JOdcUtOmTSvvJ/foAW+9BYcfXnnfU0REpJpEcrivB1oX+rwV8JNHtYiIiESNSA73/wEdzaytmcUDQ4FXPa5JREQk4sV5XUBJnHN5ZjYBWAj4gZnOuVUelyUiIhLxIjbcAZxzbwJvel2HiIhINInkaXkREREpB4W7iIhIjFG4i4iIxBiFu4iISIxRuIuIiMQYhbuIiEiMUbiLiIjEGIW7iIhIjFG4i4iIxBhzzu37UVHAzLYAP1Tit2wC/FyJ389LOpfIpHOJTDqXyKRzKd6hzrm9tkWNmXCvbGa2xDmX5HUdlUHnEpl0LpFJ5xKZdC77R9PyIiIiMUbhLiIiEmMU7iWb4XUBlUjnEpl0LpFJ5xKZdC77QdfcRUREYoxG7iIiIjFG4V4MM+tnZmvN7Bszu97revaHmbU2s3fNbLWZrTKzKwqOTzWzH81sRcFbf69rLQszSzezzwtqXlJwrJGZLTKzrwveN/S6zn0xs8MK/b9fYWbbzGxitDwvZjbTzDab2ReFjhX7PFjIAwWvn8/MrLt3le+thHO5x8zWFNT7spk1KDieaGa7Cj0/j3pX+d5KOJcSf6fMbHLB87LWzE73purilXAuzxc6j3QzW1FwPNKfl5L+Dlffa8Y5p7dCb4Af+BZoB8QDK4HOXte1H/U3B7oXfHwg8BXQGZgKXO11feU4n3SgSZFjdwPXF3x8PXCX13Xu5zn5gY3AodHyvAAnA92BL/b1PAD9gfmAAccDn3hdfxnOpS8QV/DxXYXOJbHw4yLtrYRzKfZ3quDvwEogAWhb8HfO7/U5lHYuRb7+T+CWKHleSvo7XG2vGY3c93Ys8I1z7jvnXA7wHDDQ45rKzDm3wTm3rODj7cBqoKW3VVW6gcCTBR8/CZzjYS3l0Rv41jlXmU2XqpRz7n3glyKHS3oeBgJPuZCPgQZm1rx6Kt234s7FOfeWcy6v4NOPgVbVXlg5lPC8lGQg8JxzLts59z3wDaG/dxGhtHMxMwOGAHOqtahyKuXvcLW9ZhTue2sJZBT6fD1RGo5mlgh0Az4pODShYMpnZjRMZRdwwFtmttTMxhYcO8Q5twFCLyLgYM+qK5+h7PlHKhqfFyj5eYj211AqoVFUWFszW25m75lZT6+K2k/F/U5F8/PSE9jknPu60LGoeF6K/B2utteMwn1vVsyxqLulwMwOAOYCE51z24DpQHvgaGADoSmuaHCic647cAZwmZmd7HVBFWFm8cDZwP8VHIrW56U0UfsaMrMbgTzg2YJDG4A2zrluwCRgtpkd5FV9ZVTS71TUPi/Ahez5D+KoeF6K+Ttc4kOLOVah50bhvrf1QOtCn7cCfvKolnIxs1qEfqGedc69BOCc2+Scy3fOBYHHiKDpuNI4534qeL8ZeJlQ3ZvCU1YF7zd7V+F+OwNY5pzbBNH7vBQo6XmIyteQmY0ABgB/dgUXQgumsDMLPl5K6Dr1n7yrct9K+Z2K1uclDhgMPB8+Fg3PS3F/h6nG14zCfW//AzqaWduCUdZQ4FWPayqzgmtTjwOrnXP3Fjpe+PrNIOCLov9tpDGzemZ2YPhjQoueviD0fIwoeNgIYJ43FZbLHiOQaHxeCinpeXgVGF6wAvh44LfwVGSkMrN+wHXA2c65nYWONzUzf8HH7YCOwHfeVFk2pfxOvQoMNbMEM2tL6Fw+re76yqEPsMY5tz58INKfl5L+DlOdrxmvVxVG4huhlYtfEfrX4I1e17OftZ9EaDrnM2BFwVt/4Gng84LjrwLNva61DOfSjtDq3pXAqvBzATQGFgNfF7xv5HWtZTyfukAmUL/Qsah4Xgj9g2QDkEtolDG6pOeB0BTjwwWvn8+BJK/rL8O5fEPommf4NfNowWPPLfjdWwksA87yuv4ynEuJv1PAjQXPy1rgDK/r39e5FBx/AhhX5LGR/ryU9He42l4z6lAnIiISYzQtLyIiEmMU7iIiIjFG4S4iIhJjFO4iIiIxRuEuIiISYxTuIiIiMUbhLiIiEmMU7iJSbmbWyswu8LoOEdmTwl1EKqI3oT24RSSCqEOdiJSLmZ1EqDf2VmA7MMiF9gkXEY8p3EWk3MxsAXC1cy6aNrwRiXmalheRijiM0CYkIhJBFO4iUi5m1pjQ1pS5XtciIntSuItIebUFfvK6CBHZm8JdRMprDdDEzL4wsxO8LkZE/qAFdSIiIjFGI3cREZEYo3AXERGJMQp3ERGRGKNwFxERiTEKdxERkRijcBcREYkxCncREZEYo3AXERGJMf8P2mmu+A1T+GoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(8, 6))\n", "plt.plot(t, np.sqrt(mean_square_dist), 'g.', t, np.sqrt(t), 'r-')\n", "plt.xlabel(r\"$t$\")\n", "plt.ylabel(r\"$\\sqrt{\\langle (\\delta x)^2 \\rangle}$\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "hide_code_all_hidden": false, "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.8.3" } }, "nbformat": 4, "nbformat_minor": 2 }