{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Some python techniques useful for your assignment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "IHE Transient groundwater\n", "@Theo Olsthoorn\n", "2019-01-17\n", "2019-12-21" ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2019-01-17T20:55:11.471174Z", "start_time": "2019-01-17T20:55:11.468247Z" } }, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To help you get up and going with your assignment, I explain in this notebook some features of Python and some techniques in Python that will prove useful." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.684154Z", "start_time": "2019-01-17T22:53:50.849637Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.special import exp1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Strings (unchangable series of characters)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.691967Z", "start_time": "2019-01-17T22:53:51.687207Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is a string\n" ] } ], "source": [ "a = 'This is a string'\n", "print(a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function `dir(a)` gives all objects associated with `a`, a string. Show only the public ones, i.e. the ones that do not start with '_'.\n", "\n", "Here are the functions that can be applied on strings. But there is none with which you can change an existing string; you can, however, replace your string it with a new one." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.700309Z", "start_time": "2019-01-17T22:53:51.695284Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['capitalize', 'casefold', 'center', 'count', 'encode', 'endswith', 'expandtabs', 'find', 'format', 'format_map', 'index', 'isalnum', 'isalpha', 'isascii', 'isdecimal', 'isdigit', 'isidentifier', 'islower', 'isnumeric', 'isprintable', 'isspace', 'istitle', 'isupper', 'join', 'ljust', 'lower', 'lstrip', 'maketrans', 'partition', 'removeprefix', 'removesuffix', 'replace', 'rfind', 'rindex', 'rjust', 'rpartition', 'rsplit', 'rstrip', 'split', 'splitlines', 'startswith', 'strip', 'swapcase', 'title', 'translate', 'upper', 'zfill']\n" ] } ], "source": [ "# I use a list comprehension to generate a list from an iterable (see further down\n", "# for explanation)\n", "methods = [x for x in dir(a) if not x.startswith('_')]\n", "print(methods)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now apply some of these so-called string methods on the string `a`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.819628Z", "start_time": "2019-01-17T22:53:51.815285Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This Is A String , result is a with capitalized first letter.\n", "This is a string , a is still unchanged.\n" ] } ], "source": [ "print(a.title(), \", result is a with capitalized first letter.\")\n", "# Note you get a new string, a is still the same\n", "print(a, \", a is still unchanged.\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.887541Z", "start_time": "2019-01-17T22:53:51.879991Z" } }, "outputs": [ { "data": { "text/plain": [ "'THIS IS A STRING'" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.upper()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:51.944613Z", "start_time": "2019-01-17T22:53:51.938955Z" } }, "outputs": [ { "data": { "text/plain": [ "['This', 'is', 'a', 'string']" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.split() # split the string on whitespace" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.008925Z", "start_time": "2019-01-17T22:53:52.003477Z" } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.startswith('John')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.074796Z", "start_time": "2019-01-17T22:53:52.069390Z" } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.endswith('ing')" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.139030Z", "start_time": "2019-01-17T22:53:52.133599Z" } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a.upper().endswith('ING')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Insert values into a string on successive locations indicated by {}. This is done by applying the function `format` on the string. The aguments of `format(..)` are the values that will be formatted and placed in the successive placeholders `{}`.\n", "\n", "We can put code inside these placeholders to specify how each value should be formatted. Read the documentation to find out. There is a wealth of possibilities." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.429986Z", "start_time": "2019-01-17T22:53:52.423729Z" } }, "outputs": [ { "data": { "text/plain": [ "'Today is Monday 3 Feb in the year 2019.'" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "'Today is {} {} {} in the year {}.'.format('Monday', 3, 'Feb', 2019)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tuples (unchangeble lists)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.691858Z", "start_time": "2019-01-17T22:53:52.686853Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1, 2, 'a', 'string', array([9, 3, 4]))\n" ] } ], "source": [ "# a tuple is a list with parentheses ()\n", "a = (1, 2, 'a', 'string', np.array([9, 3, 4]))\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:52.759535Z", "start_time": "2019-01-17T22:53:52.753886Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 , first\n", "2 , second\n", "[9 3 4] , last\n", "string , last but one\n" ] } ], "source": [ "# Items can be accessed by referring to them using a numerix index\n", "print(a[0], ', first')\n", "print(a[1], ', second')\n", "print(a[-1], ', last')\n", "print(a[-2], ', last but one')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Selection of more than one item at once to get a sub-tuple, is done by by\n", "slicing like so: a[:] = 'all' of tuple `a`, `a[:5]` means `a` from start up to but not including item 5, `a[4:]` means `a` from item 4 to the last item.\n", "\n", "Examples" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:53.052604Z", "start_time": "2019-01-17T22:53:53.044365Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "('a', 'string', array([9, 3, 4])) , third to last\n", "(1, 2) , frist up to last but 3\n", "1 , first item\n", "[9 3 4] , last item\n", "('a', 'string', array([9, 3, 4])) , last 3 items\n" ] } ], "source": [ "print(a[2:], ', third to last')\n", "print(a[:-3], ', frist up to last but 3')\n", "print(a[0], ', first item')\n", "print(a[-1], ', last item')\n", "print(a[-3:], ', last 3 items')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By slicing (with the colon :) you always get a tuple back, unless the result is only a single item, then you get the item itself a[0:1] would yield the first item (0 up to but not including 1)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A tuple consisting of only one item is indicated with a comma:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:53.543553Z", "start_time": "2019-01-17T22:53:53.538013Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(2,) b is a tuple, due to the comma\n", "2 c is just aa number, not a tuple\n" ] } ], "source": [ "b = (2,)\n", "c = (2)\n", "print(b, ' b is a tuple, due to the comma')\n", "print(c, ' c is just aa number, not a tuple')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A series of items separated by comma's is also a tuple, like in the arguments of a function." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:53.906456Z", "start_time": "2019-01-17T22:53:53.900409Z" } }, "outputs": [ { "data": { "text/plain": [ "(1, 2, 3, 4)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1, 2, 3, 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A single number or any single object ended with a comma indicates a tuple of one item" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:54.215790Z", "start_time": "2019-01-17T22:53:54.209775Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(3,) \n", "('sheep',) \n", "sheep \n" ] } ], "source": [ "a = 3,\n", "b = 'sheep',\n", "c = 'sheep' # no comma here\n", "print(a, type(a))\n", "print(b, type(b))\n", "print(c, type(c))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lists are just like tuples, but can be changed in place" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lists are demarcated by straight brackets `[ ]` instead of parentheses `( )`. But they work further exactly like tuples:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:54.724509Z", "start_time": "2019-01-17T22:53:54.718709Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 'a', 'string', array([9, 3, 4])]\n" ] } ], "source": [ "a = [1, 2, 'a', 'string', np.array([9, 3, 4])] # a list, not a tuple, because it has [..]\n", "print(a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Don't use a comma to indicate a list a list of only one item.\n", "\n", "`[3]` means a list of one item (item = 3)\n", "`[3,]` means turn the tuple `(3,)` into a list.\n", "It is namely the same as `[(3,)]` because `3,` is same as `(3,)` as was explained above." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:55.006894Z", "start_time": "2019-01-17T22:53:55.002282Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3]\n", "[3]\n" ] } ], "source": [ "b = [3,] # turn tuple (3,) into a list\n", "c = [3] # put item 3 into a list\n", "print(b)\n", "print(c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A list is changable in place. So you can insert items and extend it or sort it, all in place, which you cannot with a tuple, which you have to replace." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:55.313698Z", "start_time": "2019-01-17T22:53:55.307957Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['The', 'quick', 'fox', [['cow', 'chick'], 'bull', 'John', array([4, 3, 2])]]\n" ] } ], "source": [ "a = ['The', 'quick', 'fox']\n", "a.append([['cow', 'chick'], 'bull', 'John', np.array([4, 3, 2])])\n", "print(a) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note we have appended an entire list to the list `a` and we did so in-place, meaning that `a` itself has now changed. \n", "\n", "Notice that a now consists of 4 items. First three strings and then an entire list as the fourth item, which itself also contains items and even a list. But `a` has only four items after the `append()`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is what you can do with a list (dir(a) gives the methods associated with list a, so that the can be applied on `a` with the dot `.`, like `a.count()`, `a.append(..)` etc.):" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:55.835863Z", "start_time": "2019-01-17T22:53:55.830916Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['append', 'clear', 'copy', 'count', 'extend', 'index', 'insert', 'pop', 'remove', 'reverse', 'sort']\n" ] } ], "source": [ "c = [x for x in dir(a) if not x.startswith('_')]\n", "print(c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is what you can do with a tuple" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:56.135057Z", "start_time": "2019-01-17T22:53:56.128635Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['The', 'quick', 'fox', [['cow', 'chick'], 'bull', 'John', array([4, 3, 2])]]\n", "('The', 'quick', 'fox', [['cow', 'chick'], 'bull', 'John', array([4, 3, 2])])\n", "\n", "Methods associated with tuples:\n", "['count', 'index']\n" ] } ], "source": [ "b = tuple(a) # assign to b the tuple version of a (but that does not change a)\n", "print(a)\n", "print(b)\n", "\n", "# Methods associated with a tuple of which `c` is an instantiation\n", "print()\n", "print('Methods associated with tuples:')\n", "c = [x for x in dir(b) if not x.startswith('_')]\n", "print(c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Only count and index, no append insert or pop or sort, i.e. anything that could alter the tuple.\n", "\n", "Hence, trying to append something tuple `b` won't work. But we can replace `b`." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:56.437283Z", "start_time": "2019-01-17T22:53:56.433354Z" } }, "outputs": [], "source": [ "#('one', 'two').append('cow') # this gives an error, because a tuple has no attribute 'append'\n", "['one', 'two'].append('three') # works" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, because append changes the list in place, it does not produce output so we have no result unless we assign the list to a variable." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:56.748847Z", "start_time": "2019-01-17T22:53:56.743010Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "None , This is \"b\", it is a \"None\" because a.append() changes \"a\" in place and yields no output.\n", "['one', 'two', 'cow'] , This is \"a\", it was changed in place.\n" ] } ], "source": [ "a = ['one', 'two']\n", "b = a.append('cow')\n", "print(b, ', This is \"b\", it is a \"None\" because a.append()',\n", " 'changes \"a\" in place and yields no output.')\n", "print(a, ', This is \"a\", it was changed in place.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numpy arrays, so-called ndarrays (n-dimensional arrays)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Arrays may look like lists, but they are not. In a lists each item can be of a different type, as was shown above, but in an array all items must be of the same type like floating point numbers. An array forms one contiguous block in memory, in which each item occupies the same number of bytes (a floating point number occupies 8 bytes, an integer 4 bytes). This means that items can be accessed at full speed using indexing and slicing. You can also add arrays of the same size together, take the power on a whole array or other functions etc. You cannot do this with a list, because items in the list may be of completely different types. In fact, the items in the list are stored on different locations in memory, while the list is just a number of pointers that indicate these places in memory so that they can be accessed indirectly, yet still quite fast." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you multiply a list, you copy the list and append it like so:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:57.537210Z", "start_time": "2019-01-17T22:53:57.532082Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5 times list a = [3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse']\n", "\n", "a + a + a = [3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse', 3, 2, 5, 'a', 'horse']\n" ] } ], "source": [ "a = [3, 2, 5, 'a', 'horse']\n", "print('5 times list a =', 5 * a)\n", "print()\n", "print('a + a + a = ', a + a + a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But if you multiply an array, then each item in the array is multiplied as expected.\n", "\n", "First generate an array (a numpy ndarray):" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:57.824540Z", "start_time": "2019-01-17T22:53:57.820132Z" } }, "outputs": [], "source": [ "b = np.array([3, 4, 5, 2, 9, 8])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`b` is not a list, but an array, because `np.array(..)` function turns the given list into an array." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.129407Z", "start_time": "2019-01-17T22:53:58.124854Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n" ] } ], "source": [ "print(type(a))\n", "print(type(b))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.195968Z", "start_time": "2019-01-17T22:53:58.191562Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5 times b = [15 20 25 10 45 40]\n" ] } ], "source": [ "print('5 times b =', 5 * b)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.286931Z", "start_time": "2019-01-17T22:53:58.251816Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b + b + b = [ 9 12 15 6 27 24]\n" ] } ], "source": [ "print('b + b + b =', b + b + b)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.319385Z", "start_time": "2019-01-17T22:53:58.314490Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b * 2 * b**3 = [ 162 512 1250 32 13122 8192]\n" ] } ], "source": [ "print('b * 2 * b**3 = ', b * 2 * b ** 3)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.381788Z", "start_time": "2019-01-17T22:53:58.376841Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.73205081 2. 2.23606798 1.41421356 3. 2.82842712]\n" ] } ], "source": [ "print(np.sqrt(b))" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.450186Z", "start_time": "2019-01-17T22:53:58.444270Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.02021246 0.02694853 0.03368336 0.01347549 0.06060436 0.05387748]\n" ] } ], "source": [ "print(np.sin(b * np.exp(-5)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some ways to generate arrays" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.768760Z", "start_time": "2019-01-17T22:53:58.764048Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3 2 4 5]\n" ] } ], "source": [ "# Tell numpy to convert a list or a tuple into an array\n", "a = np.array([3, 2, 4, 5])\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.829885Z", "start_time": "2019-01-17T22:53:58.824533Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[3. 3.3 3.6 3.9 4.2 4.5 4.8 5.1 5.4 5.7 6. 6.3 6.6 6.9]\n" ] } ], "source": [ "a = np.arange(3, 7, 0.3) # start at 3 upto (not including 7, with steps of 0.3)\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.895511Z", "start_time": "2019-01-17T22:53:58.890559Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 1.3 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7 4. ]\n" ] } ], "source": [ "a = np.linspace(1, 4, 11) # give 11 points starting at 1 and ending at 4\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:58.961649Z", "start_time": "2019-01-17T22:53:58.955194Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.01 0.01258925 0.01584893 0.01995262 0.02511886 0.03162278\n", " 0.03981072 0.05011872 0.06309573 0.07943282 0.1 0.12589254\n", " 0.15848932 0.19952623 0.25118864 0.31622777 0.39810717 0.50118723\n", " 0.63095734 0.79432823 1. 1.25892541 1.58489319 1.99526231\n", " 2.51188643 3.16227766 3.98107171 5.01187234 6.30957344 7.94328235\n", " 10. ]\n" ] } ], "source": [ "a = np.logspace(-2, 1, 31) # 10**(-2) to 10**1 in 30 steps (31 values)\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:59.024368Z", "start_time": "2019-01-17T22:53:59.018686Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0. 0.]\n" ] } ], "source": [ "a = np.zeros_like(a)\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:59.097335Z", "start_time": "2019-01-17T22:53:59.091946Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[5. 5.]\n", " [5. 5.]\n", " [5. 5.]]\n" ] } ], "source": [ "a = 5 * np.ones((3, 2)) # 2D array with 3 rows and 2 columns\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:59.164534Z", "start_time": "2019-01-17T22:53:59.160010Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[5. 3.14159265]\n", " [5. 3.14159265]\n", " [5. 3.14159265]]\n" ] } ], "source": [ "a[:, 1] = np.pi # replace all rows of column 2 (with index 1) by np.pi\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:59.235131Z", "start_time": "2019-01-17T22:53:59.230051Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-1.2246468e-15 7.7685322e-01]\n", " [-1.2246468e-15 7.7685322e-01]\n", " [-1.2246468e-15 7.7685322e-01]]\n" ] } ], "source": [ "b = np.sin(2 * np.pi * a) # apply function on all values in array\n", "print(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## For-loops with zipping and enumerate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A for loop cycles over an iterable, which is something you can iterate over like a `list`, a `tuple`, an `array`, a `string`. Here are some examples." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:53:59.885666Z", "start_time": "2019-01-17T22:53:59.876943Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 3 5 4 6 2 2 " ] } ], "source": [ "a = [1, 3, 5, 4, 6, 2, 2] # a list, this is a so-called iterable\n", "b = np.array(a) # turn list a into an array, an array is also an iterable\n", "c = ['a', 'bear', 'c', 'd', 'echo', 'f', 'gulf'] # list of strings\n", "for x in a:\n", " print(x, end=' ')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In each new cycle of the loop, the next item of the iterble is taken and can be used inside the loop. But very often we have a set of iterable like the array of x-coordinates of the wells, the array of y-coordinates of the wells and the list of extractions from the wells. In each cycle of the loop we need the next `x`, the next `y` and the next `Q` at the same time. To deal with such a set of iterables in a loop just put them in a `zip(..)` function." ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:00.290895Z", "start_time": "2019-01-17T22:54:00.284776Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 1, 'a']\n", "[3, 3, 'bear']\n", "[5, 5, 'c']\n", "[4, 4, 'd']\n", "[6, 6, 'echo']\n", "[2, 2, 'f']\n", "[2, 2, 'gulf']\n" ] } ], "source": [ "for aa, bb, cc in zip(a, b, c): # a, b, c are iterables, aa, bb and cc are items\n", " print([aa, bb, cc])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oftentimes you not only want the next times from the iterable but also have a counter, an index while looping. This counter can be obtained by packing the `iterable` or the entire `zip(..)` with all its iterables in an `enumerate(..)` function. Note that counting starts with 0 in python, not 1." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:00.648674Z", "start_time": "2019-01-17T22:54:00.642826Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0, 1, 1, 'a']\n", "[1, 3, 3, 'bear']\n", "[2, 5, 5, 'c']\n", "[3, 4, 4, 'd']\n", "[4, 6, 6, 'echo']\n", "[5, 2, 2, 'f']\n", "[6, 2, 2, 'gulf']\n" ] } ], "source": [ "for i, (aa, bb, cc) in enumerate(zip(a, b, c)):\n", " print([i, aa, bb, cc])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This type of looping is efficient and extremely useful in practice. For instance we have 5 wells so 5 x vales, 5 y values 5 flows for which we want the `drawdown s` at a given point `x0`, `y0`" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:01.285385Z", "start_time": "2019-01-17T22:54:01.010271Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa6ElEQVR4nO3dfZQV9Z3n8fcnPAgDCCoKgg8Qgw+sg08djDNqQHYWNcchcUyGjKOJicNxcnyIO3HVzYnRzLhJJptVE3VY1jiZZDJhYnQNKDNklqURBx9QUAwqEVFjC6yCgjbhqeG7f1Q13Nu53V1Nd93bt+vzOueervpV3arv/R2oz62HW6WIwMzMiutDtS7AzMxqy0FgZlZwDgIzs4JzEJiZFZyDwMys4BwEZmYF5yCwQpD0Q0l/09fXaXYgHARmZgXnILBCk9S/1jWY1ZqDwPokSadJWiHpA0n/DAxK26dIapJ0o6SNwN9LOkTSI5LekfReOnxUOv9USS+ULPf/SHq6ZPxxSZ/saJ0l8/6FpLWS3pU0T9KYtP02Sd9PhwdI2ibpb9PxwZJ2pDWOkxSSPifpN5I2Sfpqnv1oxeAgsD5H0kDgYeDHwKHAA8CflMwyOm0/FphF8v/g79PxY4DtwN3pvE8AH5E0Mt17OBk4StIwSYOBM4Clna1T0nnAN4HPAEcCbwBz08lLgCnp8EeBjcDH0/GzgDUR8V5J/WcDJwDTgFskndTFLjIr4yCwvuhjwADgzojYHRE/B5aXTN8LfD0idkbE9ojYHBEPRsRvI+ID4HbSDXFE7ACeAc4FGoBVwOPAH6breSUiNmdY56XA/RGxIiJ2AjcDZ0kaRxI2EyQdlq7nB8BYSUPTOpa0+Xy3pXU/DzwPnNL9LrMi8/FR64vGAG9F+R0V3ygZfifdwAMg6feAO4DzgUPS5mGS+kXEHvZ/Y29Kh98j2UDvZP9GurN1jgFWtI5ERLOkzcDYiHhd0jPpMs8lCaJTScLm48D323y+jSXDvwWGttsTZhl4j8D6og0k36hV0nZMyXDbW+7+FcmhljMj4mCSjTFA6/tbg+DcdHgJyQa69Nt6Z+tcT3LoKVmwNAQ4DHirZB3nAaeR7EksAaYDk4HHOvvAZt3hILC+6AmgBbhWUn9JF5NsUNszjOS8wBZJhwJfbzN9GUlQTAaejojVJBv1M9m/ke5snf8EXCHpVEkHAf8NeCoiXk+nLwEuB16MiF1AI3Al8FpEvNPVDjDrCgeB9TnphvRi4PMkh3H+FHiog7fcCQwGNgFPAv/aZnnbSA7rrE6XDcmG/42IeDvLOiNiEfA14EGSvYfjgJklq1mW1tAaLC8CO/DegFWB/GAaM7Ni8x6BmVnBOQjMzArOQWBmVnAOAjOzgqu7H5SNHDkyxo0bV+syumXbtm0MGTKk1mX0Gu6Pcu6P/dwX5brTH88+++ymiDi80rS6C4Jx48bxzDPP1LqMbmlsbGTKlCm1LqPXcH+Uc3/s574o153+kPRGe9N8aMjMrOAcBGZmBecgMDMrOAeBmVnBOQjMzAqu7q4aMjMrmlWrVrFo0SK2bt3KypUrmTZtGpMmTeqx5TsIzMx6sVWrVjF//nx2794NwNatW5k/fz5Aj4WBDw2ZmfViixYt2hcCrXbv3s2iRYt6bB0OAjOzXmzr1q1daj8QDgIzs15s+PDhXWo/EA4CM7NebNq0aQwYMKCsbcCAAUybNq3H1uGTxWZmvVjrCeHWq4aGDx/uq4bMzIpm0qRJTJo0Kbeb8PnQkJlZwTkIzMwKzkFgZlZwDgIzs4JzEJiZFZyDwMys4Hz5qJlZO7atfJv3F77Oni076TfiIA6ePo4hpx1R67J6nIPAzKyCbSvfZstDrxC79wKwZ8tOtjz0CkCfCwMfGjIzq+D9ha/vC4FWsXsv7y98vTYF5chBYGZWwZ4tO7vUXs9yDQJJ50taI2mtpJsqTJ8iaauk59LXLXnWY2aWVb8RB3WpvZ7lFgSS+gH3ABcAE4HPSppYYdalEXFq+vpGXvWYmXXFwdPHoQHlm0gN+BAHTx9Xm4JylOcewWRgbUSsi4hdwFxgRo7rMzPrMUNOO4IRF0/YtwfQb8RBjLh4Qp87UQz5XjU0FnizZLwJOLPCfGdJeh5YD3wlIlbnWJOZWWZDTjuiT27428ozCFShLdqMrwCOjYhmSRcCDwMTfmdB0ixgFsCoUaNobGzs2UqrrLm5ue4/Q09yf5Rzf+znviiXV3/kGQRNwNEl40eRfOvfJyLeLxleIOleSSMjYlOb+eYAcwAaGhoij/txV1Ne9xSvV+6Pcu6P/dwX5erxeQTLgQmSxksaCMwE5pXOIGm0JKXDk9N6NudYk5mZtZHbHkFEtEi6GlgI9APuj4jVkq5Kp88GLgH+UlILsB2YGRFtDx+ZmVmOcr3FREQsABa0aZtdMnw3cHeeNZiZWcd8ryEzqysvLV3M0rk/4oPNmxh22EjOmXk5J50ztdZl1TUHgZnVjZeWLuaXc+6mZVdym4cPNr3DL+ckBxUcBgfO9xoys7qxdO6P9oVAq5ZdO1k690c1qqhv8B6BmfU6v35qI0/84lWa393LG7/8d86acRzHnzmaDzZvqjh/e+2WjYPAzHqVXz+1kcU/eZmWXcktoJvf3cnin7wMwLDDRvLBpnd+5z3DDhtZ1Rr7Gh8aMrNe5YlfvLovBFq17NrLE794lXNmXk7/geV3/+w/8CDOmXl5NUvsc7xHYGa9SvO7le/33/zuTk465zwAXzXUwxwEZtarDD30oIphMPTQZE/gpHOmesPfw3xoyMx6lbNmHEf/geWbpv4DP8RZM46rUUV9n/cIzKxXOf7M0QDpVUM7GXroQfuuGrJ8OAjMrNc5/szRHH/m6PRum39Y63L6PB8aMjMrOAeBmVnBOQjMzArOQWBmVnAOAjOzgnMQmJkVnIPAzKzgHARmZgXnIDAzKzgHgZlZwTkIzMwKzkFgZlZwDgIzs4JzEJiZFZyDwMys4BwEZmYF5yAwMys4B4GZWcE5CMzMCs5BYGZWcLkGgaTzJa2RtFbSTR3M91FJeyRdkmc9Zmb2u3ILAkn9gHuAC4CJwGclTWxnvm8DC/OqxczM2pfnHsFkYG1ErIuIXcBcYEaF+a4BHgTezrEWMzNrR55BMBZ4s2S8KW3bR9JY4FPA7BzrMDOzDvTPcdmq0BZtxu8EboyIPVKl2dMFSbOAWQCjRo2isbGxh0qsjebm5rr/DD3J/VHO/bGf+6JcXv2RZxA0AUeXjB8FrG8zTwMwNw2BkcCFkloi4uHSmSJiDjAHoKGhIaZMmZJTydXR2NhIvX+GnuT+KOf+2M99US6v/sgzCJYDEySNB94CZgJ/VjpDRIxvHZb0Q+CRtiFgZmb5yi0IIqJF0tUkVwP1A+6PiNWSrkqn+7yAmVkv0GEQSJqXYRnvRsTnK02IiAXAgjZtFQOgvWWYmVm+OtsjOAm4soPpIvmtgJmZ1anOguCrEbGkoxkk3daD9ZiZWZV1+DuCiPhZZwvIMo+ZmfVemU4WS2oAvgocm75HQETEpBxrMzOzKsh61dBPgBuAF4C9+ZVjZmbVljUI3omILFcQmZlZnckaBF+XdB+wCNjZ2hgRD+VSlZmZVU3WILgCOBEYwP5DQwE4CMzM6lzWIDglIn4/10rMzKwmst6G+slKD5UxM7P6l3WP4Gzgc5JeIzlH4MtHzcz6iKxBcH6uVZiZWc1kCoKIeCPvQszMrDY6PEcgaUVnC8gyj5mZ9V6d3n1U0qoOpgsY3oP1mJlZlXUWBCdmWMaenijEzMxqo7O7j76Rnh/oD2xMh8cDM4Ct6fSmKtRpZmY5yfo7ggeBPZI+AvyAJAz+KbeqzMysarIGwd6IaAEuBu6MiOuBI/Mry8zMqiVrEOyW9FngcuCRtG1APiWZVdfW+fN55bxpvHTSRF45bxpb58+vdUlmVZU1CK4AzgJuj4jXJI0H/jG/ssyqY+v8+Wz42i20rF8PEbSsX8+Gr93iMLBCyRQEEfFiRFwbET9Nx1+LiG/lW5pZ/t6+405ix46yttixg7fvuLM2BZnVQNY9ArM+qWXDhi61m/VFDgIrtP5HVr7mob12s77IQWCFdsT1X0aDBpW1adAgjrj+y7UpyKwGMt10TtLxJA+vP7b0PRFxXk51mVXF8IsuApJzBS0bNtD/yCM54vov72s3K4Kst6F+AJgN/C98SwnrY4ZfdJE3/FZoWYOgJSL+LtdKzMysJrKeI5gv6UuSjpR0aOsr18rMzKwqsu4RfC79e0NJWwAf7tlyzMys2rI+oWx83oWYmVltZDo0JGmppNslnS9pWNaFp/OvkbRW0k0Vps+QtErSc5KekXR2V4o3M7Puy3qO4HPAGuBPgGXpRvuOjt4gqR9wD3ABMBH4rKSJbWZbBJwSEacCXwDu60LtZmbWA7IeGlonaTuwK31NBU7q5G2TgbURsQ5A0lySB9q8WLLc5pL5h5CcdzAzsyrKemjoVeBhYBTJg2lOjojzO3nbWODNkvGmtK3tsj8l6WXgUZK9AjMzqyJFdP4lXNJ1wNnA0cDLwBLgsYh4tYP3fBqYHhFXpuOXAZMj4pp25j8XuCUi/mOFabOAWQCjRo06Y+7cuZ3W3Js1NzczdOjQWpfRa7g/yrk/9nNflOtOf0ydOvXZiGioODEiMr+AocA1wBvAnk7mPQtYWDJ+M3BzJ+95DRjZ0TxnnHFG1LvFixfXuoRexf1Rzv2xn/uiXHf6A3gm2tmuZj009F1JTwFPAacCtwATOnnbcmCCpPGSBgIzgXltlvsRSUqHTwcGApuz1GRmZj0j6w/KngT+NiL+X9YFR0SLpKuBhUA/4P6IWC3pqnT6bJKrkC6XtBvYDvxpmlxmZlYlWa8aekDSH6fH8QGWRESnz/KLiAXAgjZts0uGvw18uwv1mplZD8t6aOibwHUkl36+CFybtpmZWZ3LemjoE8CpEbEXQNI/ACtJTgCbmVkd68oTykaUDA/v4TrMzKxGsu4RfBNYKWkxIOBcvDdgZtYnZD1Z/FNJjcBHSYLgxojYmGdhZmZWHR0GQXptf6mm9O8YSWMiYkU+ZZmZWbV0tkfw3fTvIKABeJ5kj2ASyY/LfNtoM7M61+HJ4oiYGhFTSW4pcXpENETEGcBpwNpqFGhmZvnKetXQiRHxQutIRPyK5FYTZmZW57JeNfSSpPuAfyR5ZsCfAy/lVpWZmVVN1iC4AvhLkl8XAzwG/F0uFZmZWVVlDYI/AP5nRHT4eEozM6s/WYPg88BsSZuBpenr8Yh4L6/CzMysOrL+oOxyAEljgEtIHko/Juv7zcys98q0IZf058A5wO8Dm4C7SfYKzMyszmX9Rn8n8CowG1gcEa/nVZCZmVVXpt8RRMRI4AskvzC+XdLTkn6ca2VmZlYVWR9MczBwDHAsMI7kNtR78yvLzMyqJeuhocdLXndHRFMn85uZWZ3IetXQpLwLMTOz2sh61dDhwH8B/gPJeQIAIuK8nOoyM7MqyXrTuZ8ALwPjgduA14HlOdVkZmZVlDUIDouIHwC7I2JJRHwB+FiOdZmZWZVkPVm8O/27QdIngPXAUfmUZGZm1ZQ1CP5G0nDgr4DvAwcD1+dWlZmZVU2nQSCpHzAhIh4BtgJTc6/KzMyqptNzBBGxB/jjKtRiZmY1kPXQ0DJJdwP/DGxrbYyIFblUZWZmVdOVB9MAfKOkLQD/jqAOPbruUe5acRcbt21k9JDRXHf6dXziw5+odVlmViNZf1ns8wJ9xKPrHuXWZbeyY88OADZs28Cty24FcBiYFVSHQSDpP3c0PSL+R8+WY3m7a8Vd+0Kg1Y49O7hrxV0OArOC6uxk8bD01UDy8Pqx6esqYGJnC5d0vqQ1ktZKuqnC9EslrUpfyySd0vWPYF2xcdvGLrWbWd/X4R5BRNwGIOmXwOkR8UE6fivwQEfvTS87vQf4I6AJWC5pXkS8WDLba8DHI+I9SRcAc4AzD/CzWAajh4xmw7YNFdvNrJiy3mLiGGBXyfgukucSdGQysDYi1kXELmAuMKN0hohYFhHvpaNP4l8r5+66069jUL9BZW2D+g3iutOvq1FFZlZrWa8a+jHwtKT/TXK10KeAf+jkPWOBN0vGm+j42/4XgX/JWI8doNbzAL5qyMxaKSKyzSidTvIAe4DHImJlJ/N/GpgeEVem45cBkyPimgrzTgXuBc6OiM0Vps8CZgGMGjXqjLlz52aqubdqbm5m6NChtS6j13B/lHN/7Oe+KNed/pg6deqzEdFQaVrWPYLWH4915QdkTcDRJeNHkdysroykScB9wAWVQiBd9xyS8wc0NDTElClTulBG79PY2Ei9f4ae5P4o5/7Yz31RLq/+yHqO4EAsByZIGi9pIDATmFc6g6RjgIeAyyLi1znWYmZm7ci8R9BVEdEi6WpgIdAPuD8iVku6Kp0+G7gFOAy4VxJAS3u7LmZmlo/cggAgIhYAC9q0zS4ZvhK4Ms8azMysY3keGjIzszrgIDAzKzgHQb1Z9TO442S4dUTyd9XPal2RmdW5XM8RWA9b9TOYfy3s3p6Mb30zGQeY9Jna1WVmdc17BPVk0Tf2h0Cr3duTdjOzA+QgqCdbm7rWbmaWgYOgngxv55587bWbmWXgIKgn026BAYPL2wYMTtrNzA6Qg6CeTPoMXPQ9GH40oOTvRd/ziWIz6xZfNVRvJn3GG34z61HeIzAzKzgHgZlZwTkIzMwKzkFgZlZwDgIzs4JzEJiZFZyDwMys4BwEZmYF5yAwMys4B4GZWcE5CMzMCs5BYGZWcA4CM7OCcxCYmRWcg8DMrOAcBGZmBecgMDMrOAeBmVnBOQjMzArOQWBmVnAOAjOzgss1CCSdL2mNpLWSbqow/URJT0jaKekredZiZmaV9c9rwZL6AfcAfwQ0AcslzYuIF0tmexe4FvhkXnWYmVnH8twjmAysjYh1EbELmAvMKJ0hIt6OiOXA7hzrMDOzDuQZBGOBN0vGm9I2MzPrRXI7NASoQlsc0IKkWcAsgFGjRtHY2NiNsmqvubm57j9DT3J/lHN/7Oe+KJdXf+QZBE3A0SXjRwHrD2RBETEHmAPQ0NAQU6ZM6XZxtdTY2Ei9f4ae5P4o5/7Yz31RLq/+yPPQ0HJggqTxkgYCM4F5Oa7PzMwOQG57BBHRIulqYCHQD7g/IlZLuiqdPlvSaOAZ4GBgr6QvAxMj4v286jIzs3J5HhoiIhYAC9q0zS4Z3khyyMjMzGrEvyw2Mys4B4GZWcE5CMzMCs5BYGZWcA4CM7OCcxCYmRWcg8DMrOAcBGZmBecgMDMrOAeBmVnBOQjMzAou13sN9RYPr3yL7yxcw/ot2xkzYjA3TD+BT57mZ+SYmUEBguDhlW9x80MvsH33HgDe2rKdmx96AcBhYGZGAQ4NfWfhmn0h0Gr77j18Z+GaGlVkZta79PkgWL9le5fazcyKps8HwZgRg7vUbmZWNH0+CG6YfgKDB/Qraxs8oB83TD+hRhWZmfUuff5kcesJYV81ZGZWWZ8PAkjCwBt+M7PK+vyhITMz65iDwMys4BwEZmYF5yAwMys4B4GZWcEpImpdQ5dIegd4o9Z1dNNIYFOti+hF3B/l3B/7uS/Kdac/jo2IwytNqLsg6AskPRMRDbWuo7dwf5Rzf+znviiXV3/40JCZWcE5CMzMCs5BUBtzal1AL+P+KOf+2M99US6X/vA5AjOzgvMegZlZwTkIzMwKzkGQI0nnS1ojaa2kmypMv1TSqvS1TNIptaizGjrri5L5Pippj6RLqllftWXpD0lTJD0nabWkJdWusZoy/F8ZLmm+pOfT/riiFnVWg6T7Jb0t6VftTJek76V9tUrS6d1eaUT4lcML6Ae8CnwYGAg8D0xsM88fAIekwxcAT9W67lr1Rcl8/xdYAFxS67pr/G9jBPAicEw6fkSt665xf/xX4Nvp8OHAu8DAWteeU3+cC5wO/Kqd6RcC/wII+FhPbDe8R5CfycDaiFgXEbuAucCM0hkiYllEvJeOPgkcVeUaq6XTvkhdAzwIvF3N4mogS3/8GfBQRPwGICL6cp9k6Y8AhkkSMJQkCFqqW2Z1RMRjJJ+vPTOAH0XiSWCEpCO7s04HQX7GAm+WjDelbe35IknK90Wd9oWkscCngNlVrKtWsvzbOB44RFKjpGclXV616qovS3/cDZwErAdeAK6LiL3VKa/X6eq2pVOFeEJZjahCW8VrdSVNJQmCs3OtqHay9MWdwI0RsSf50tenZemP/sAZwDRgMPCEpCcj4td5F1cDWfpjOvAccB5wHPBvkpZGxPs519YbZd62ZOUgyE8TcHTJ+FEk32bKSJoE3AdcEBGbq1RbtWXpiwZgbhoCI4ELJbVExMNVqbC6svRHE7ApIrYB2yQ9BpwC9MUgyNIfVwDfiuQg+VpJrwEnAk9Xp8ReJdO2pSt8aCg/y4EJksZLGgjMBOaVziDpGOAh4LI++k2vVad9ERHjI2JcRIwDfg58qY+GAGToD+AXwDmS+kv6PeBM4KUq11ktWfrjNyR7R0gaBZwArKtqlb3HPODy9OqhjwFbI2JDdxboPYKcRESLpKuBhSRXRdwfEaslXZVOnw3cAhwG3Jt+E26JPninxYx9URhZ+iMiXpL0r8AqYC9wX0RUvJyw3mX89/HXwA8lvUByaOTGiOiTt6eW9FNgCjBSUhPwdWAA7OuLBSRXDq0Ffkuyt9S9daaXI5mZWUH50JCZWcE5CMzMCs5BYGZWcA4CM7OCcxCYmRWcg8DMrOAcBFZokkZI+lLJ+BhJP89hPbdKekvSN9qZ/rqkkZIGp7ee3iVpZE/XYVaJg8CKbgSwLwgiYn1E5PUshDsi4paOZoiI7RFxKt28ZYBZV/iXxVZ03wKOk/Qc8G/APcAjEXGypM8DnyT5tevJwHdJ7pd/GbATuDAi3pV0XPq+w0l+6fkXEfFyRyuVdBjw0/Q9T1P5RmJmVeE9Aiu6m4BXI+LUiLihwvSTSZ4NMBm4HfhtRJwGPAG03hp6DnBNRJwBfAW4N8N6vw48ni5rHnBM9z6G2YHzHoFZxxZHxAfAB5K2AvPT9heASZKGkjxp7oGS22cflGG55wIXA0TEo5Le62R+s9w4CMw6trNkeG/J+F6S/z8fArakx/W7yjf6sl7Bh4as6D4Ahh3om9MHo7wm6dOw78Hip2R462PApel7LgAOOdAazLrLQWCFlj4M6N8l/UrSdw5wMZcCX5T0PLCays9jbus24FxJK4D/RHK/fbOa8G2ozapA0q1Ac0T894zzvw409NV77lvv4j0Cs+poBma194OyVq0/KCN5EElRH85uVeY9AjOzgvMegZlZwTkIzMwKzkFgZlZwDgIzs4L7/5SWMdW6tinOAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "kD = 600 # m2/d\n", "S = 0.1\n", "t = 1.3 # d\n", "xWells = [34, 123, 45. -70, -80.]\n", "yWells = [-23., 45., -50., -30., 20.]\n", "Qwells = [1200, 600., 350., 1200., 300]\n", "\n", "x0, y0 = 0.3, 6.7 # multiple assignment, handy\n", "\n", "times = [0.1, 0.25, 0.3, 0.5, 0.8, 0.82, 0.9, 1.0] # days\n", "\n", "# plot individual times as a dot\n", "plt.title('drawdown')\n", "plt.xlabel('time [d]')\n", "plt.ylabel('drawdown s[m]')\n", "plt.grid()\n", "\n", "for t in times:\n", " s = 0\n", " for xw, yw, Qw in zip(xWells, yWells, Qwells):\n", " r = np.sqrt((x0 - xw)**2 + (y0 - yw)**2)\n", " s += Qw / (4 * np.pi * kD) * exp1(r**2 * S / (4 * kD * t))\n", " plt.plot(t, s, 'o') # plot the point as a dot, color is automatic\n", " \n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Superposition, using logical indexing to handle switch times" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When dealing with superposition, we have a set of so called `switch times` at which the `boundary conditions` change and we have an `array of times` with values at small intervals to obtain sufficient values to produce a nice, detailed graph.\n", "\n", "During the superposition loop, we need to obtain the results of for instance the head change `s` at times larger than the the switch time, e.g. `st3` when the well was swithed on. We compute its results for times `t - st3`, but only for those times for which `t > st3`.\n", "\n", "We can select those times by logical indexing. That is, if `t` is an array of times and `ts3` is a single time, then `t > ts3` is an array of booleans, an array of True/False values of the same size as the array `t`, such that the elements that are True correspond with the condition `t > ts3`. We can then address the correct values of an array `s`, which has the same size as the array `t`. We do this like so `s[t > ts3]`. Or, in two steps: `I = t > ts3` and then select the right values like so `s[I]`." ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:01.678018Z", "start_time": "2019-01-17T22:54:01.670954Z" } }, "outputs": [], "source": [ "times = np.linspace(0, 360, 361) # from 0 to 360, 361 points (1 d interval)\n", "swt = [0, 30, 34, 60, 90, 96, 110, 150, 200, 220, 300, 305, 330]\n", "Q = np.array([300, 0, 200, 50, -30, 120, 50, 60, 150, 35, -50, 100, 0])\n", "\n", "# changes of flow. hstack((a, b)) glues two arrays or a value and array together\n", "dQ = np.hstack((Q[0], np.diff(Q)))\n", "#print('dQ = ', dQ)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:02.162151Z", "start_time": "2019-01-17T22:54:01.888135Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEWCAYAAAAgpUMxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6dklEQVR4nO29eZwU1bn//34YBoZ9HWAYkEURg+y4RjFj1G/ASNDk672amLgkIZr4zWa8EpNvfhrjV2LM1eTGGy4xC5oYNG6gF6/iMi5xA3TYRUAQBgYGhnVkWOf5/VE1bU1P93R1T/VUd83zfr361adOnXPqeaqrz6eeU6eqRFUxDMMwjNamXdgGGIZhGG0TEyDDMAwjFEyADMMwjFAwATIMwzBCwQTIMAzDCAUTIMMwDCMUTICMQBGRW0XkgWbWf0VEns/CdieLyNoA25stIv83qPaiiIgMFREVkfbucrmIfCNsu8KkJceNuy9PCtqmXEbsPqDWRUTOBe4GTgWOA2uA76vq4lANywIiMhTYCBSq6rGQzUmKiFwDfENVzw3blnwi/vcVkXLgr6qa9ASkLSEiZTj7Y5DP8gqMUNX12bQrl2gftgFtCRHpDjwD3AA8CnQAJgOHW9kOwTn5qG/N7RqGYTRCVe3TSh/gNGBvM+tvwzljalgeCijQ3l0uB+4C3gH2AfOB3p7yZwFvAHuBZUCZZ105cCfwT6AOOMlHe18AVrntlQOf8qy7BdgKHADWAhfE+wBsdu2vdT9nA9cAr3va+TSw2N3+YuDTcTbf4dp8AHge6Jtk35UBlZ7lTcCPgOVu248ARQnqfQo4hBON1jb8PsBfgF942wb+DagGqoBLgYuBD4DdwK2eNtsBM4ENQA3OyUbvRHYnsKcjcB+wzf3cB3SMs+Mmjx3XxtW9x93vO4DZQKck2/kImOSmr3J/p1Hu8jeAp1L5QuLj8xs+/SwHfoFzvNYCTwN9gL8B+91jYain/G+ALe66pcDkuP/No8CD7nGyCjjNXXcz8Hjctv8DuC+JXU2Oa6AI5z/T1y3zU+AY0N1d/kVDew3HDdDFrVPPJ8f/QKAAuNXdnwdcXwa7dRW4HlgH7AHuxx2lctdfhzNisgd4Dhji5gtwr3tM7MM55keH3d/5Og7CNqAtfYDu7p94LjAV6BW3/jZSC9BWYLR7gD/OJ519qdv2xW6ncZG7XOypuxln6K89UJiivZOBj912CnE63/U4UdtItzMY6LHzxHgf4u13867BFSCgt/tn+qpr05Xuch+PzRtcWzq5y7OS7NsymgrQO+6fvrf7x70+Sd2YTZ68v9BYgI4BP3P3xTeBncDDQDd3nx4Chrvlvw+8BQzCEYX/Av7u8xj5uVu3H1CM00HfEWfHz107LgYO4h5HOGK1wPW3G06nfleS7TwI3OSm57j7+QbPuh+k8iX+9yV9AVoPnAj0AFbjiPmF7rHwIPBnT/mrcASqPY4Ab8c9ocA55g65+6MA56TqLXddCc5x3NNdbo/TUU9KYFNzx/WrwJfc9PPu/prqWXdZkuOmMm4bNwMr3G0JMI5PjnfFGSHpCZyAc4xNcddd6u6vT7k+/BR4w133ORwh6+m2+SmgJOz+ztdxELYBbe3jHhx/wTmTPYbTYfR3191GagGa5Vk/Cjji/uluAR6K29ZzwNWeuj+PW99ce/8XeNSzrh2OWJXhRE/VOJ1FYVybMR/i7XfzruETAfoq8E5c/TeBazz2/dSz7tvA/yTZr43+7DgCdJVn+W5gdpK6MZs8eX+hcUdSBxS4y91cv870lF8KXOqm1+BGhO5yCXDUux+aOT42ABd7lj8HbIqzw7s/q3EiX8HpaE/0rDsb2JhkO18HFnjs/QYwz13+CJiYypf435f0BegnnuVfA896lqcBFc3U3wOM8xxzL8Qdx3We5WeBb7rpS4DVSdps7ri+A/it6/d24HvALJpGR/HHTbwArQWmJ9m+Aud6lh8FZnp8+Hrc//EgMAT4LI54nwW087P/c+Vjs+BaGVVdo6rXqHNhcjTOGfp9aTSxxZP+COdMuC/OgXi5iOxt+ADn4nQYieqmam+gu9xgd71btlSdi6Tfx/njV4vIPBEZmIYPDTTahseGUs/ydk/6INA1jfZbUjeeGlU97qbr3O8dnvV1nvaHAE96foc1OEN8/X1sJ36ffOTmee3wTuho8KsY6Aws9Wz3f9z8RLwCTBaRATgnHI8A57gTC3oAFQH4kor4/ZdsfyIiN4nIGhHZ59rRA+c4bSD+ty5qmJ2HM+JwlZu+CngokTEpjutXcARlIk4Eswj4DE6nv15Vd6V2F4DBOCcZyUh2zA4BfuP5HXbjnHSUqupLwO9whux2iMgc93pzzmMCFCKq+j7OGdNoN+tjnE6kgQEJqg32pE/AORvdhSMOD6lqT8+ni6rO8m4yjfa24Rz0QGziwmCcKAhVfVidWWND3HZ/mcjFBHleGm3DY8PWFPWCJpWd6bIFZ3jG+1sUqaofv+L3yQluXip24XTap3q22UNVE4qu29keBL4LvKqqB3A6vxk40WDDBJWW+BIIIjIZJ8L/F5zhxp441zrEZxNPAWNFZDROBPS3ZAWbOa7fwBk2uwx4RVVX4/w2n8cRp4TNJcjbgjPsmC5bgG/F/Q6dVPUN1+7fquoknOHgk3GG+nIeE6BWREROcc/kBrnLg3Gue7zlFqkAzhORE0SkB/DjBM1cJSKjRKQzzrWAx9wz878C00TkcyJSICJFIlLWsK1mSNbeo8DnReQCESnEGXc/DLwhIiNF5LMi0hFn7L0O56w4np04F2GHJ9n2QuBkEfmyiLQXkX/FGT55JoXNQbMDGCQiHQJqbzZwp4gMARCRYhGZ3rBSRDa5U78T8Xfgp26dvjjXnf6aaoOuYPwBuFdE+rnbKRWRzzVT7RXgRj7pQMvjllP6kgzPPUJDU5X1QTec4eqdQHsR+RnO9VRfqOoh4DGca3bvqOrmROWaO65V9SDOMOt3+GT/vAF8i+QCtAPo4/6XG3gAuENERojDWBHp48ON2cCPReRU19YeInK5mz5dRM50/6cf88mkmpzHBKh1OQCcCbwtIh/jCM9KnM4dVV2EMxSyHOdgT9QRP4QTNW3HGX/+rlt3CzAdZ4bNTpwzpptJ/Rsna28tznDFf+CcXU8DpqnqEZyL0bPc/O04F8xvjW/Y/dPeCfzTHTo4K259Dc4Z6U04Eyb+DbgkjeGMoHgJZ+bUdhEJYtu/wbm297yIHMD5nc8EcEWuD5+cdMTzC2AJzjGwAnjXzfPDLTgXqt8Skf3ACzhn7cl4BadzfzXJcrO+pGAwzvBhEJHSczjXQD5w2zxE4uHk5pgLjCHJ8JtLquP6FZwh6nc8y/H7K4Y7wvF34EP3+B8I/DvOyd3zODP6/ogzwaZZVPVJnGhsnvvbrsSZyASOGP8B57rYRzj/pXtStZkL2I2oeUTQN/rZjYOtj3sj8ndU9cqwbckmIvJTYKeq/lfYtgCIyAnA+8AAVd0ftj2Gg92IahitiKq+Drweth3ZRlX9Rm1ZR0TaAT/EmeVn4pNDmAAZhhFZRKQLzrWYj4ApIZtjxGFDcIZhGEYo2CQEwzAMIxRsCA7o2bOnnnRSdJ+C/vHHH9OlS5ewzcga5l9+E2X/ouwbwNKlS3eparKbnVNiAgT079+fJUuWhG1G1igvL6esrCxsM7KG+ZffRNm/KPsGICLxTzJJCxuCMwzDMELBBMgwDMMIBRMgwzAMIxTsGpBhGEaGHD16lMrKSg4dOpRwfY8ePVizZk0rWxU8RUVFDBo0iMLCwkDbNQEyDMPIkMrKSrp168bQoUNxHhjfmAMHDtCtW7cQLAsOVaWmpobKykqGDRsWaNs2BGcYhpEhhw4dok+fPgnFJyqICH369Eka5bUEE6AAqK6uZvPmhE94Nwwj4kRZfBrIlo8mQC1k+fLlDB06lBNPPJEXX3wxbHMMwzDyBhOgFjJ58mTq6uo4duwY999/f9jmGIbRxti7dy//+Z//2WyZTZs28fDDD6dsa9OmTYwePTpluaAwAWoBBw4cYP/+T57uXl1dHaI1hmG0RYIUoNbGBKgFPPHEE42WJ0yYEJIlhmG0VWbOnMmGDRsYP348N998MzfffDOjR49mzJgxPPLII7Eyr732GuPHj+fee+9l06ZNTJ48mYkTJzJx4kTeeOONUGy3adgt4O9//3ujZXu1hWG0XbI5GaG5vmXWrFmsXLmSiooKHn/8cWbPns2yZcvYtWsXp59+Oueddx6zZs3innvu4ZlnngHg4MGDLFq0iKKiItatW8eVV14ZyvMwTYAyZNeuXbzwwguN8kyADMMIk9dff50rr7ySgoIC+vfvz2c+8xkWL15M9+7dG5U7evQoN954IxUVFRQUFPDBBx+EYq8JUIY8+eSTHD9+vFGeCZBhGGHitw+699576d+/P8uWLaO+vp6ioqIsW5YYuwaUIQ1jq15MgAyj7aKqTT779+9PmJ/upzm6devGgQMHADjvvPN45JFHOH78ODt37uTVV1/ljDPOaFQGYN++fZSUlNCuXTseeuihJifTrYVFQBlQXV3Nyy+/3CTfBMgwjNamT58+nHPOOYwePZqpU6cyduxYxo0bh4hw9913M2DAAPr06UP79u0ZN24c11xzDd/+9rf50pe+xD/+8Q/OP//80F6aZwKUAc888wz19fVN8k2ADMMIg/gp1r/61a8aLRcWFja5UX758uWx9F133QXA0KFDWblyZZasbIoNwWXA008/HUuPGDEiREsMwzDyFxOgNDl06BDPP/98bHnatGmxtEVAhmEY/jEBSpOXXnqJgwcPAk70M3LkyNg6EyDDaHu0hf99tnwMVYBEZIqIrBWR9SIyM8H6U0TkTRE5LCI/8uQPFpGXRWSNiKwSke951t0mIltFpML9XBykzd7ht2nTpjW6+awtHIiGYXxCUVERNTU1kf7vN7wPKBtTtUObhCAiBcD9wEVAJbBYRBao6mpPsd3Ad4FL46ofA25S1XdFpBuwVEQWeereq6r3BG2zqvLf//3fseVp06axfv36RusNw2g7DBo0iMrKSnbu3Jlw/aFDh0K7xyZIGt6IGjRhzoI7A1ivqh8CiMg8YDoQEyBVrQaqReTz3oqqWgVUuekDIrIGKPXWzQbr1q1jy5YtgDP3/pxzzmHDhg1eu7K5ecMwcozCwsJm3xJaXl5uz4hshjAFqBTY4lmuBM5MtxERGQpMAN72ZN8oIl8DluBESnsS1JsBzAAoLi6mvLw85bbmz58fS48ZM4Z//vOfrF27NpZXVVXlq53Wpra2NiftCgrzL7+Jsn9R9i0IwhSgRE/uSyuEEJGuwOPA91W14b0IvwfucNu6A/g1cF2TDanOAeYAjBw5UsvKylJu73e/+10s/S//8i+UlZWxcePGWF7//v3x005rU15enpN2BYX5l99E2b8o+xYEYU5CqAQGe5YHAdv8VhaRQhzx+Zuqxt6LoKo7VPW4qtYDf8AZ6msxx48f56WXXootX3jhhQ12xPJsCM4wDMM/YQrQYmCEiAwTkQ7AFcACPxXF6fX/CKxR1X+PW1fiWbwMCOS23oqKCvbscUbyBgwYwKhRoxq2FytjAmQYhuGf0IbgVPWYiNwIPAcUAH9S1VUicr27fraIDMC5jtMdqBeR7wOjgLHAV4EVIlLhNnmrqi4E7haR8ThDcJuAbwVhr/cxFhdccEFMeEyADMMwMiPUZ8G5grEwLm+2J70dZ2guntdJfA0JVf1qkDY28Nprr8XS559/fixtAmQYhpEZ9iQEH6gqb775Zmz53HPPjaWz+RZEwzCMKGMC5IN169ZRU1MDQO/evTn55JMTlrMIyDAMwz8mQD544403Yumzzz67UdRjQ3CGYRiZYQLkA68AffrTn260zgTIMAwjM0yAfOC9/mMCZBiGEQwmQCnYu3cvq1atAqCgoIDTTz+90XoTIMMwjMwwAUrBkiVLYsIyduzYJu9ONwEyDMPIDBOgFFRUVMTS8dEPmAAZhmFkiglQCrwCNH78+CbrTYAMwzAywwQoBa0hQPPnz2fUqFF86Utf4vjx4xm1YRiGkW+YADVDXV0d77//PuAIzZgxY5qUaakAPfDAA3zxi19kzZo1PPHEE40e+WMYhhFlQn0WXK6zatWqWEQyYsQIunbt2qRMSwTo7rvv5pZbbmmUt2/fvgwsNRJx8OBBli9fTkVFBe+99x7Lli2jW7duPPDAAwwZMiRs8wyjzWMC1AzvvfdeLJ1o+A0yexacqnLrrbcya9asJuvq6+vTbs+AnTt38t5778XEpqKigg8++CDh/pw9ezZ33XVXCFYahuHFBKgZUl3/icdPBKSq/OQnP0koPn7bMD5h7ty53HzzzezcudN3nb1792bPIMMwfGMC1Ax+BCjdIbjbb7+90dn3tGnTOH78OAsXOm+lsAjIP/X19Xz3u99l//79Cde3a9eOkSNHMn78eA4cOMAzzzwTq2cYRviYACVBVVmxYkVsedy4cQnLpSNAv/zlL7n99ttjy5dccgmPPfYYX/nKV3y3YXzCsWPHGonPmWeeyfjx45kwYQLjx49nzJgxdO7cGYA5c+bEBMj2sWHkBiZASdi2bRsHDhwAoGfPnpSUlCQs51eAHnroIWbOnBlbnjJlCo899hgdOnRo1IadnfvHu6+Kiop46623kpa1fWwYuYdNw05Cw/RrgFNOOSXpZAM/AvTCCy9w3XXXxZbLysp44okn6NixI+AMFTVgnaN/vPvKuw8T4V1vEZBh5AYmQEmIF6BkpBKgFStW8MUvfpFjx44BMGbMGJ566ik6derkuw0jMekIkEVAhpF7hCpAIjJFRNaKyHoRmZlg/Ski8qaIHBaRH/mpKyK9RWSRiKxzv3tlYtuaNWti6U996lPN+RBLx4vHnj17uPTSS2NDeaWlpSxcuJAePXo0KmcRUGZYBGQY+U1oAiQiBcD9wFRgFHCliIyKK7Yb+C5wTxp1ZwIvquoI4EV3OW1aGgEdP36cr3zlK3z44YcAdO3alWeffZZBgwb5bsNonkwFyETeMHKDMCOgM4D1qvqhqh4B5gHTvQVUtVpVFwNH06g7HZjrpucCl2ZiXEsF6Pbbb+fZZ5+NLc+dOzfho3zAOsdMyXQIzkTeMHKDMGfBlQJbPMuVwJkB1O2vqlUAqlolIv0SNSAiM4AZAMXFxZSXl8fW1dXVsXXrVsB5Cd3mzZvZtm1bQkO8U7VramooLy/n7bff5o477ojlf/nLX6Z3796NtuGluro6ll6zZk3ScplSW1sbeJu5gPexRcePH2/WR+8JRVVVVV7tj6j+fg1E2b8o+xYEYQpQomllfk9NW1LXKaw6B5gDMHLkSC0rK4ut84rK8OHDufDCC5O2c/jw4Vi6V69enHrqqVxxxRWxvAsvvJAHH3yQgoKCpG3MnTs3lj755JPx2hIE5eXlgbeZC3iFu2PHjs366D2B6NevX17tj6j+fg1E2b8o+xYEYQ7BVQKDPcuDgMRhRnp1d4hICYD7XU2aNFy3ARg2bFizZeOHdr75zW+yY8cOAPr378/DDz/crPiAXSDPFLsGZBj5TZgR0GJghIgMA7YCVwBfDqDuAuBqYJb7PT9dwzZu3BhLDx8+3He9F154odHyn//8Z4qLi1PWsynCmWHXgAwjvwlNgFT1mIjcCDwHFAB/UtVVInK9u362iAwAlgDdgXoR+T4wSlX3J6rrNj0LeFREvg5sBi5P1zZvBJRKgJLdoHrDDTcwdepUX9uzCCgzLAIyjPwm1EfxqOpCYGFc3mxPejvO8Jqvum5+DXBBS+wKQoAmT57se3sWAWWGRUCGkd/YkxASkM4QnJ9H9KTCIqDMsAjIMPIbE6A4VLWRAKUzCcFLqg4xWRvWOfrHHsVjGPmNCVAcO3bsoK6uDnCmVffs2bPZ8kFHQNY5+scexWMY+Y0JUBybN2+OpYcMGZKyfBACZNcnMsMiIMPIb0yA4qisrIylEz23LZ4ghuAsAsoMi4AMI78xAYqj4RE80DIBskkI2ccmIRhGfmMCFIdXgEpLS1OWD3oIzjpH/3j3Var9bcOchpF7mADF4R2Ca4kAZToEZ52jfywCMoz8xgQojjCG4CwCygy7EdUw8hsToDjSHYJLhl0Dyj4WARlGfmMC5EFVQxmCswgoMywCMoz8xgTIw969e2M3oXbt2pXu3bunrGOz4MLDIiDDyG9MgDzED7/5ERG7BhQeFgEZRn5jAuTB+9bMgQMH+qpjN6KGh0VAhpHfmAB58L7iuX///r7q2BBceFgEZBj5jQmQh507d8bSft5kCjYEFyYWARlGfmMC5MErQH379vVVx25EDQ97GKlh5DfNvhFVRCb6aOOoqq4IyJ5Q2bVrVyxtEVDuYw8jbcrx48fZvXs3O3fupLq6mp07d3LkyBGmTp1K7969wzbPMBqR6pXcrwCLgeZ602HA0KAMCpOwhuDaSucYNG0tAqqqquL1119n586djQTGm66pqUno3+TJk3n11VdDsNowkpNKgBar6mebKyAiL2W6cRGZAvwGKAAeUNVZcevFXX8xcBC4RlXfFZGRwCOeosOBn6nqfSJyG/BNoEFNblXVhX7sCVKAcvlG1IMHD7JixQpGjBiR12fFbSkCWrt2LePGjePw4cMZ1X/nnXcCtsgwWk6zApRKfPyWSYSIFAD3AxcBlcBiEVmgqqs9xaYCI9zPmcDvgTNVdS0w3tPOVuBJT717VfWedG2KcgRUXV3N008/zfz581m0aBGHDh2itLSU9evXU1RUlJVtZpu2NAlh0aJFvsWnZ8+e9OvXjz59+vDmm28C+Sm6RvRp0TUgVX23Bds+A1ivqh+625oHTAe8AjQdeFCdf89bItJTREpUtcpT5gJgg6p+1AJbgMwEKBm5cA1o7dq1zJ8/nwcffJDVq1c36YS2bt3KypUrOe200wLbZmvSlqZhHz9+PJYeO3Ys06ZNo7i4mH79+lFcXBxL9+3bl8LCQgCOHj1Khw4dgPwUXSP6pBqC+7X7XQScBizDuR40FngbOLcF2y4FtniWK3GinFRlSgGvAF0B/D2u3o0i8jVgCXCTqu6J37iIzABmgCM2L7zwAnv27GlYx/LlyykoKEjpxPr16xPmV1RUcOzYsZT1ATZu3BhLb968mfLycl/1muO+++5j/vz5Kcu988471NbWtnh7YbBs2bJYes+ePc3uN+/vtH///kD2cWtRW1vLunXrYssjRozgwgsvbFRmz5497Nmzh7Vr18byvKJVX1+fsz7X1tbmrG0tJcq+BUGqIbjzIRadzGiY7SYio4EftXDbiUKE+FPTZsuISAfgC8CPPet/D9zhlrsDR0Sva9KI6hxgDsDIkSN19OjRsXW9e/fmggsu8OVEsmsokyZN4uyzz/bVxpIlS2Lp0tJSysrKfNVLRl1dHU8//XSjvHbt2nHOOecwffp07r777thNt+PGjeOcc85p0fbCYv/+/bF0cXFxs/utV69esXSXLl1avI9bk/LycoYPHx5bHjx4sC/7vVGPquasz+Xl5TlrW0uJsm9BkCoCauAU71RrVV0pIuNbuO1KYLBneRCwLc0yU4F3VXWHx7ZYWkT+ADzjx5hMh9+CvgYUxFDJkSNHYu0UFhbywx/+kJtuuinm11NPPRUTIO9Zcr7Rlq4BeYcN/R5b+T7saEQfv1O11ojIAyJSJiKfcTv2NS3c9mJghIgMcyOZK4AFcWUWAF8Th7OAfXHXf64kbvhNREo8i5cBK/0YE7QAhXkjqreD7dy5M1OmTGnkk3doMR874wba0jWglgpQfBuGkQv4jYCuBW4Avucuv4oz1JUxqnpMRG4EnsOZhv0nVV0lIte762cDC3GmYK/HmYZ9bUN9EemMM4PuW3FN3+1GZwpsSrA+IWFGQEFPQkjVWXkFyCKg/CATAQLH7wZ/6+vrfV3XNIzWwpcAqeoh4F73Exju/TkL4/Jme9IKfCdJ3YNAnwT5X83Elt27d8fS6dwbk4vTsE2AmhKlCCjTe8zy0W8j2jR7JIvInFQN+CmTD3gvaPfs2dN3vVyMgFJ1zG1RgNpyBNRAPvptRJtUEdClInKomfUCnB+gPaHhFSA/b0JtIBevAVkE1JR8jwRMgIwokkqAbvbRxmtBGBI2+/bti6V79Ojhu16uR0AmQE3X52NHnKkA5bvwGtEm1X1Ac1vLkLAJOgLKlWtANgTnkO8dsUVARhSx9wG5BB0Bhfkw0lSdlde2tiJA+d4Rp4pqk5HvfhvRxgTIJUoRkE1CaEq+v47BhuCMKGIC5OKNgNIRoGTYfUDZpy29jsGG4Iwo4us+IBE5GWdCwhBvnUxfxZCLeCOg1h6CC7qTMAFqSr53xCZARhTx+ySEfwCzgT8A+dtjNYMNweUfbXUSgt2IakQFvwJ0TFVb9OidXCdK07DTiYDy+azYIqDU5LvfRrTxeyr1tIh8W0RKRKR3wyerlrUiqsqRI0cAaN++fVpvCM3FG1HtPqCm5HskYJMQjCjiNwK62v323piqwPAEZfMOb0fWo0ePjP/gfvJTlQ06ArIhuKbr8zESsAjIiCJ+H0Y6LNuGhIn3j5nuDLhcvAZk9wE1Jd8jARMgI4r4nQX3Gs4rGF4D/qmqB7JqVSvj7YTTuf4DuXkjqk1CaEq+d8SZ3oia78JrRBu/veTVwFrgS8AbIrJERAJ9NUOYtLUIqC0KUL53xBYBGVHE7xDchyJSBxxxP+cDn8qmYa1J/DWgdMj1a0AmQE3X52NHbAJkRBFfEZCIbACeAvoDfwRGq+qULNrVqmQjAsqVWXA2BOfQViOgfPfbiDZ+e8nfApuBK4HvAleLyIlZs6qVCTsCsichZEZbjYDCfMqGYQSJryNZVX+jqpcDFwJLgduAD7JoV6vSkggoGWGepVoE1JR8jwRsCM6IIn6H4H4tIm8DbwPjgZ8BI1q6cRGZIiJrRWS9iMxMsF5E5Lfu+uUiMtGzbpOIrBCRChFZ4snvLSKLRGSd+90rlR3eP2bXrl3T9SFhvj0LLvu01QjIhuCMqOD3RtS3gLtVdUdQGxaRAuB+4CKgElgsIgtUdbWn2FQcoRsBnAn83v1u4HxV3RXX9EzgRVWd5YraTOCW5mzxdkidO3dO14+08lOVbe33AeVjZ9yARUCpicpvbUQTv0Nw/wDOFJF73M+0ALZ9BrBeVT9U1SPAPGB6XJnpwIPq8BbQU0RKUrQ7HWh4k+tc4NJUhnj/3GEIkE1CyAyLgFKT734b0cbvjah34QjG39ys74rIp1X1xy3YdimwxbNcSePoJlmZUqAK51FAz4uIAv+lqnPcMv1VtQpAVatEpF8Sn2YAMwA6duwYy//oo48oLy/37cTevXsT5r/55pv07NnTVxsrVqyIpXft2pXW9hOxdu3aWPrjjz+mtra2UZsfffRRLL1x48YWby8sNm7cGEun+t2OHTsWSx8/fjyvfK6trWXz5s2x5Q0bNvi2v66uLpZ+++232bEjsEGMwIg/PqNElH0LAr9DcJ8HxqtqPYCIzAXeA1oiQIlO4+JP/5src46qbnMFZpGIvK+qr/rduCtYcwC6deumhw8fBmDSpEmUlZX5bYZdu+JHAB3OPfdc+vbt66uNQ4cOxdK9evVKa/uJ6NKlSyzdrVs3unbt2qjNxYsXx9KlpaUt3l5YPP/887H0iSee2KwfXgEC8srn8vJySktLY8sjRozwbb/3muZpp53G6NGjgzavxZSXl+fV75EOUfYtCNJ5I2pPTzq9ucqJqQQGe5YHAdv8llHVhu9q4EmcCA1gR8MwnftdncqQKF8DsiG4puvzcSjKJiEYUcSvAN0FvCcif3Gjn6XA/2vhthcDI0RkmIh0AK4AFsSVWQB8zZ0Ndxawzx1W6yIi3QBEpAvwv4CVnjoNT+++GpifyhDvH7NTp05pOZGLN6LaLLim5HtHbNeAjCji91E8fxeRcuB0nGGxW1R1e0s2rKrHRORG4DmgAPiTqq4Skevd9bOBhcDFwHrgIHCtW70/8KT7R2wPPKyq/+OumwU8KiJfx7l59vJUtkQtArJJCE2J/z1UNa3fKGxMgIwo0qwAee+7cal0vweKyEBVfbclG1fVhTgi482b7Ukr8J0E9T4ExiVpswa4IE07YukozIKzCCgxIhLbN/ksQPZKbiMqpIqAfu1+FwGnActwIqCxODelnps901qPbERAuXwjahTeB1RfX0/DxBHwt7/btWsX87e+vj6t3yhsLAIyokizAqSq5wOIyDxghqqucJdHAz/KvnmtQzauAdmjeNJDVTlw4ADV1dXs3LmT6urqZtO7du1qNLPNbwTk3V4+YQJkRBG/07BPaRAfAFVdKSLjs2NS69OSIbhk2MNIHerr66mpqaGqqqrZz/bt2zl48GDG2+nfv3/KMvncGdssOCOK+BWgNSLyAPBXnPtwrgLWZM2qVibsIbhsTkJIJUBBd8QVFRX88pe/ZMOGDVRVVbFjxw6OHj0a6DYa6N69O8XFxYwaNYovfOELKcvnc2ec6RtR81l0jejjV4CuBW4Avucuv4rzXLZI0b59ewoLC9Oqk+uTEFp7CO7aa6+loqIio7pFRUX079+ffv360a9fP4qLi5Omi4uLKSoqApyb/fwMneZzZ2xDcEYU8StAn8Z53E1kXsOdiHSv/0BuTsMOawhOVVm1alWT/B49elBSUpLy071796zOTMvnCMiG4Iwo4leArgFmi0gN8Jr7eV1V92TLsDDwPhPOL7l4I2pYkxBqa2tjw21FRUWsWrWKAQMGBHZdraXkczRgEZARRfzeiPo1ABEZCPxvnNcoDPRbP18IUoDaYgRUU1MTS/ft25fhw4cH1nYQ5HM0EIQA5ZvPRvTx+zTsq4DJwBhgF/A7nCgoUnTo0CHtOrl+Dag17wPavXt3LN2nT5/A2g2KfI4GgrgRNd98NqKP3wjmPmADMBt4WVU3ZcugMAkrAgq6YwxrCM4bAfXu3TuwdoMin6MBG4IzoojfF9L1Ba7DeSLCnSLyjog8lFXLQiCsCCjooaFcGILLxQgon6MBm4RgRBG/Q3DdgROAIcBQnNcx5Nc/2Ae5MASXixHQ0aNHqaqqYtu2bWzfvr3ZTwO5KED5HA1YBGREEb9DcK97Pr9T1coU5fOSoIbg0p1KHGYEdPjwYTZs2EBlZWXss3Xr1kbL27dvT9uuE044IXMHskQ+RwN2I6oRRfzOghubbUNygaAioHQfchnmo3hee+01TjrppBZvs4GuXbty9tlnc9111wXWZlDkc2dsQ3BGFPE7BFcM/BtwKs51IABU9bNZsisUMomAEtGSCCiTjvHAgQOsW7eO6upqzj777JRDcA1PEEjHvgEDBsRuGB0wYECTT0lJCf3792/0CuhcI587YxuCM6KI3yG4vwGPAJcA1+O8aXRntowKi6AioHQFyO/srH379rFmzRpWrVrF6tWrWb16NatWrWLLli2xMpMmTeKnP/1ps7ZMmDCBsWPHsnz5cgoKChg4cCCDBg1K+ikpKUn7EUW5SD53xiZARhTxK0B9VPWPIvI9VX0FeEVEXsmmYWEQ1DWgdIfg4iMgVWXz5s0sXbqUd999l6VLl7JixQq2bt2asq2lS5dSW1vbrH1FRUW899577N69m169ejUakosybTECymefjejjV4AaHmdcJSKfB7YBg7JjUnjkQgS0adMmiouLG01pTkX79u0bvRvH+/TpZGLYrl07+vbtm5ad+U4+RwOZ3oiazz4b0cevAP1CRHoANwH/AXQHfpA1q0IiEwFKREuuAR05ciSp+BQWFjJy5EhOPfVURo0axahRozj11FM56aST6NevH3v37gVoJEb59NrpbJML0UB9fT2HDh2irq6OQ4cOpUwfPXqUDh062KN4jEiSUoBEpAAYoarPAPuA84PauIhMAX4DFAAPqOqsuPXirr8YOAhco6rvishg4EFgAM79SHNU9TdunduAb/LJNapbVXWhH3synYQgIhmfoQKUlpbSpUsXPv7441hez549mThxYuwzYcIETjrpJNq3T/yTeYfR0n1TaFshWTRQX19PXV0dBw8e5ODBg0nT6S7X1dXFPl5BSZdOnTpx9tlnx5bDfM6gYQRJSgFS1eMi8gUg0FcxuMJ2P3ARUAksFpEFqrraU2wqMML9nInzDqIzgWPATa4YdQOWisgiT917VfWedG3KNAKKF6B0o47OnTvzyiuvsGjRIk488UQmTZrEsGHDMj7TtQgoMd59cdZZZ3H06FEOHjzI4cOHQ7QqNXV1daxZ88n7H20SghEV/A7BvSEiv8OZCRc7TVfVd1uw7TOA9ar6IYCIzAOmA14Bmg48qE7v/paI9BSRElWtAqpcGw6IyBqgNK5u2rQkAmpu2Q+TJk1i0qRJGW0fkkdAJkCf4J3JV11dHZodRUVFFBUV0alTp2bTr7zyCrt27QIaP7HCJiEYUSGdF9IB/NyTp0BL7gMqBbZ4litxoptUZUpxxQdARIYCE4C3PeVuFJGvAUtwIqUm7y0SkRnADG9eVVUV5eXl6frRhPr6+kDaSQev6KxduzaW3rVrF7W1ta1uT2vi17/TTz+dDz74IOG6oqIiOnToQFFRER07dox9mlv2W6djx4506NCBDh06UFhY6FtA1q1bFxOgQ4cOxfJXrVrl+1FHDfUBVq5cSb9+/XzVa02ifHxG2bcg8PskhMCu+3hI9C+MP0VrtoyIdAUeB76vqvvd7N8Dd7jl7gB+jfMg1caNqM4B5rjtKMDJJ59MWVlZWk6AM8zhPUMtLCzMqJ2W4H2b65AhQ2LphptDW9ue1qS8vNyXf5/5zGe48847OXz4MJ07d6Zz586xaCMXI8VevXrF0l77xowZ4/v3nDNnTix9yimn5ORx4Pf3y0ei7FsQNCtAIvLD5tar6r+3YNuVwGDP8iCc6d2+yohIIY74/E1Vn/DYtKMhLSJ/AJ7xa1BLrgE1t9wa2BBcakSkkTjnOsne22RDcEZUSDVFqpv7OQ24AWf4qxTnaQijWrjtxcAIERkmIh2AK4AFcWUWAF8Th7OAfapa5c6O+yOwJl4ERaTEs3gZsNKvQUFdAwpj5pnNgoseQQiQTUIwcplmIyBVvR1ARJ4HJqrqAXf5NuAfLdmwqh4TkRuB53CmYf9JVVeJyPXu+tnAQpwp2OtxpmFf61Y/B/gqsEJEKty8hunWd4vIeJwhuE3At/zaFNZ9QEFgs+CihwmQEXX8TkI4ATjiWT6C816gFuEKxsK4vNmetALfSVDvdRJfH0JVv5qpPVEcgrMIKH9JJkCZvpLbhuCMXMOvAD0EvCMiT+JEFpcBc7NmVUjk8xCcRUDRwyIgI+r4nQV3p4g8C0x2s65V1feyZ1Y4RDECMgHKX5KdyNijeIyo4DcCarjptCU3nuY8Yd6I2lJsCC56BCFA9igeI5ex3slDUBGQDcEZQRB0BGQCZOQaJkAeLAIycgkbgjOijvVOHuwakJFL2BCcEXVMgDzYEJyRS9gQnBF1TIA82BCckUsEHQHZEJyRa1jv5CGfh+AsAooeyX47eyW3ERVMgDzk842odg0oetgQnBF1TIA85POz4GwILnrYEJwRdax38pDP14C8ndXRo0dDtcUIBouAjKhjAuQhn2fBWQQUPUyAjKhjvZOHqExCyPTBlUZuYUNwRtQxAfKQz0Nw3gjIhuCigUVARtQxAXJp165do048HWwIzsgG9igeI+pY7+SSafQDuREB2RBc9LBH8RhRxwTIpSVTsHNBgGwILnokEyC7EdWICiZALvkeAdkQXPSwITgj6oTaO4nIFBFZKyLrRWRmgvUiIr911y8XkYmp6opIbxFZJCLr3O9efmwJMgKyh5EaQWBDcEbU8f1G1KARkQLgfuAioBJYLCILVHW1p9hUYIT7ORP4PXBmirozgRdVdZYrTDOBW1LZYxGQkWsEHQG99tprzJo1K2N7TjzxRC6//PKM6xtGPKEJEHAGsF5VPwQQkXnAdMArQNOBB9UZO3hLRHqKSAkwtJm604Eyt/5coBwfApTv14BsEkL0CDoCevnll3n55Zcztufzn/+8CZARKGEKUCmwxbNciRPlpCpTmqJuf1WtAlDVKhHpl2jjIjIDmNGwfOTIEcrLy9P3AhgyZAhbtnxiTklJScZtZUp1dXUs/fHHH8fSH330EbW1ta1uT2sSVf+2bt2aMH/p0qXU1tb6aqN9++D+4jU1NVnZz1H9/SDavgVBmAKU6DQu/ippsjJ+6jaLqs4B5gD07t1bb7jhBsrKytJpIsaTTz7JQw89xK5duxg4cCBXXXUVPXr0yKitTJk3b14s7T1zHj58OF27ds3Yt3ygvLw8kv49++yzCfNPP/10Jk6cmHBdPOeddx4jR47k7bffbrE9p5xySlb2c1R/P4i2b0EQpgBVAoM9y4OAbT7LdGim7g4RKXGjnxKgmhQUFxdz0003pWn+J/Tt25cf/OAHGdcPAhuCix5BXQO67LLLuOyyy4IyyzACI8wr1IuBESIyTEQ6AFcAC+LKLAC+5s6GOwvY5w6vNVd3AXC1m74amJ9tR3KBZPcB2SSE/CUIATKMXCa0CEhVj4nIjcBzQAHwJ1VdJSLXu+tnAwuBi4H1wEHg2ubquk3PAh4Vka8Dm4E2cdXUXkgXPYK4EdUwcpkwh+BQ1YU4IuPNm+1JK/Adv3Xd/BrggmAtzX3sPqDoYRGQEXXsVCoi2H1A0cMEyIg61jtFhGSPXLHOKn9J9tvZb2pEBROgiJDsVRLWWeUvFgEZUccEKCIkEyAbgstfTICMqGO9U0Swzip62G9qRB0ToIhgEVD0MAEyoo71ThHBOqvoYb+pEXVMgCKCTUKIHiZARtQxAYoINgQXPexJCEbUsSM5ItjZcvSw39SIOiZAEcGG4KKHCZARdUyAIoINwUUPEyAj6ljvFBGss4oe9psaUccEKCJYBBQ9TICMqGO9U0To27dvwvzCwsJWtsQIChMgI+qYAEWESy65hH/913+lV69esbySkhIuuuiiEK0yWoIJkBF1Qn0hnREcnTp1Yt68eQAcOXKEmpoa+vbtS2FhIevXrw/ZOiMTTICMqGMCFEE6dOhASUlJ2GYYLcRuRDWijh3JhpGjWARkRJ1QBEhEeovIIhFZ5373SlJuioisFZH1IjLTk/8rEXlfRJaLyJMi0tPNHyoidSJS4X5mt5JLhhE4JkBG1AkrApoJvKiqI4AX3eVGiEgBcD8wFRgFXCkio9zVi4DRqjoW+AD4safqBlUd736uz6YThpFNTICMqBOWAE0H5rrpucClCcqcAaxX1Q9V9Qgwz62Hqj6vqsfccm8Bg7JrrmG0PiZARtQJaxJCf1WtAlDVKhHpl6BMKbDFs1wJnJmg3HXAI57lYSLyHrAf+KmqvpbIABGZAcwAKC4upry8PG0n8oXa2lrzLw9ZvXp1wvw33niDHj16tLI12SOqvx9E27cgyJoAicgLwIAEq37it4kEeRq3jZ8Ax4C/uVlVwAmqWiMik4CnRORUVd3fpCHVOcAcgJEjR2pZWZlPs/KP8vJyzL/8o6amJmH+ueeeS58+fVrZmuwR1d8Pou1bEGRNgFT1wmTrRGSHiJS40U8JUJ2gWCUw2LM8CNjmaeNq4BLgAlVVd5uHgcNueqmIbABOBpa01B/DaG1sCM6IOmFdA1oAXO2mrwbmJyizGBghIsNEpANwhVsPEZkC3AJ8QVUPNlQQkWJ38gIiMhwYAXyYNS8MI4uYABlRJywBmgVcJCLrgIvcZURkoIgsBHAnGdwIPAesAR5V1VVu/d8B3YBFcdOtzwOWi8gy4DHgelXd3VpOGUaQ2I2oRtQJZRKCqtYAFyTI3wZc7FleCCxMUO6kJO0+DjwenKWGER4WARlRx06lDCNHMQEyoo4JkGHkKCZARtQxATKMHMUEyIg6JkCGkaOYABlRxwTIMHIUEyAj6pgAGUaOYgJkRB0TIMPIUUyAjKhjAmQYOYoJkBF1TIAMI0cpKChImG8CZEQFEyDDyFEmTJjAgAGNHyg/duxYCgsLQ7LIMIIlrPcBGYaRgo4dO7J27Vruu+8+hgwZgqoyderUsM0yjMAwATKMHKZ79+6cd9559k4ZI5LYEJxhGIYRCiZAhmEYRiiYABmGYRihYAJkGIZhhIIJkGEYhhEKJkCGYRhGKJgAGYZhGKEgqhq2DaEjIgeAtWHbkUX6ArvCNiKLmH/5TZT9i7JvACNVtVumle1GVIe1qnpa2EZkCxFZYv7lL+Zf/hJl38DxryX1bQjOMAzDCAUTIMMwDCMUTIAc5oRtQJYx//Ib8y9/ibJv0EL/bBKCYRiGEQoWARmGYRihYAJkGIZhhEKbFyARmSIia0VkvYjMDNueTBCRP4lItYis9OT1FpFFIrLO/e7lWfdj19+1IvK5cKz2h4gMFpGXRWSNiKwSke+5+VHxr0hE3hGRZa5/t7v5kfAPQEQKROQ9EXnGXY6MbwAisklEVohIRcO05Kj4KCI9ReQxEXnf/Q+eHahvqtpmP0ABsAEYDnQAlgGjwrYrAz/OAyYCKz15dwMz3fRM4JduepTrZ0dgmOt/Qdg+NONbCTDRTXcDPnB9iIp/AnR104XA28BZUfHPtfmHwMPAM1E6Nj3+bQL6xuVFwkdgLvANN90B6Bmkb209AjoDWK+qH6rqEWAeMD1km9JGVV8FdsdlT8c5eHC/L/Xkz1PVw6q6EViPsx9yElWtUtV33fQBYA1QSnT8U1WtdRcL3Y8SEf9EZBDweeABT3YkfEtB3vsoIt1xTm7/CKCqR1R1LwH61tYFqBTY4lmudPOiQH9VrQKnEwf6ufl567OIDAUm4EQJkfHPHaKqAKqBRaoaJf/uA/4NqPfkRcW3BhR4XkSWisgMNy8KPg4HdgJ/dodQHxCRLgToW1sXIEmQF/V56Xnps4h0BR4Hvq+q+5srmiAvp/1T1eOqOh4YBJwhIqObKZ43/onIJUC1qi71WyVBXk76Fsc5qjoRmAp8R0TOa6ZsPvnYHmdo//eqOgH4GGfILRlp+9bWBagSGOxZHgRsC8mWoNkhIiUA7ne1m593PotIIY74/E1Vn3CzI+NfA+7wRjkwhWj4dw7wBRHZhDO8/VkR+SvR8C2Gqm5zv6uBJ3GGnaLgYyVQ6UbkAI/hCFJgvrV1AVoMjBCRYSLSAbgCWBCyTUGxALjaTV8NzPfkXyEiHUVkGDACeCcE+3whIoIzBr1GVf/dsyoq/hWLSE833Qm4EHifCPinqj9W1UGqOhTnv/WSql5FBHxrQES6iEi3hjTwv4CVRMBHVd0ObBGRkW7WBcBqgvQt7FkWYX+Ai3FmVm0AfhK2PRn68HegCjiKcxbydaAP8CKwzv3u7Sn/E9fftcDUsO1P4du5OGH8cqDC/VwcIf/GAu+5/q0EfubmR8I/j81lfDILLjK+4VwnWeZ+VjX0IVHxERgPLHGPz6eAXkH6Zo/iMQzDMEKhrQ/BGYZhGCFhAmQYhmGEggmQYRiGEQomQIZhGEYomAAZhmEYoWACZBiGYYSCCZBhBIz7CPtve5YHishjWdjObSKyVUR+nmT9JhHpKyKd3FcFHBGRvkHbYRiZYgJkGMHTE4gJkKpuU9X/naVt3auqP2uugKrWqfOsuVx95IvRRmkftgGGEUFmASe6T7heBNyP8xSA0SJyDc7j6wuA0cCvcd6z8lXgMHCxqu4WkRPdesXAQeCbqvp+cxsVkT44T8UoxnkESqKHQxpGzmARkGEEz0xgg6qOV9WbE6wfDXwZ56GVdwIH1Xna8JvA19wyc4D/o6qTgB8B/+lju/8f8Lrb1gLghJa5YRjZxSIgw2h9Xlbn5XoHRGQf8LSbvwIY67564tPAP5xnsQLOWyZTcR7wRQBV/W8R2ROs2YYRLCZAhtH6HPak6z3L9Tj/yXbAXve6TbrYwx2NvMGG4AwjeA4A3TKtrM4L9zaKyOXgvJJCRMb5qPoq8BW3zlScJxcbRs5iAmQYAaOqNcA/RWSliPwqw2a+AnxdRBoe8z/dR53bgfNE5F2c99JsznDbhtEq2OsYDCNPEZHbgFpVvcdn+U3Aaaq6K5t2GYZfLAIyjPylFpiR7EbUBhpuRAUKca4zGUZOYBGQYRiGEQoWARmGYRihYAJkGIZhhIIJkGEYhhEKJkCGYRhGKPz/qSLCEWtFWW4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xw, yw, x0, y0 = 100, 50, 0.5, -3 # assigning all at once, well and obs point coordinates\n", "r = np.sqrt((xw - x0)**2 + (yw - y0)**2)\n", "\n", "# set up the plot\n", "plt.title('Superposition in time, one well, many switches')\n", "plt.xlabel('time [d]')\n", "plt.ylabel('drawdown [d]')\n", "plt.xlim((0, 600))\n", "plt.grid()\n", "\n", "# Initialize the head changes to an array of all zeroes\n", "s = np.zeros_like(times)\n", "\n", "# Do the superposition\n", "for st, Q in zip(swt, dQ):\n", " I = times > st # logical array I telling which times > st\n", " # print(I) # shows the booleans\n", " \n", " # effect of this well only, note times[I] are used onlyu\n", " ds = Q/(4 * np.pi * kD) * exp1(r**2 * S /(4 * kD * times[I] - st))\n", " #plt.plot(times[I], ds, label='st = {:.0f} d'.format(st))\n", " \n", " s[I] += ds # add them to s pertaining to the right times[I] !!\n", "plt.plot(times, s, 'k', lw=3, label='total')\n", "plt.legend()\n", "plt.show()\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's extremely useful to get used to logical indexing of arrays. A logical index is an array of booleans (i.e. True/False values) of the same shape as the target arrays that can be used as indices. Like so" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some elegant logical indexing, which is extremely useful" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:02.567871Z", "start_time": "2019-01-17T22:54:02.561515Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[False True False True True True True False True]\n" ] } ], "source": [ "a = np.array([-3, 2, -1, 1.2, 3.2, 0.7, 3.2, -3, 5.1])\n", "I = a > 0.15 # this yields a boolean array I\n", "print(I)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then use the boolean array to pick values" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:02.960709Z", "start_time": "2019-01-17T22:54:02.953916Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2. 1.2 3.2 0.7 3.2 5.1]\n", "[2. 1.2 3.2 0.7 3.2 5.1]\n" ] } ], "source": [ "b = a[I]\n", "print(b)\n", "print(a[a > 0.15]) # read \"a where a > 0.15\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We applied logical indexing above to select the times that were greater than ts, and did something with those times only.\n", "\n", "Note that while the boolean array is as along as the target array. So boolean array `I` is as long as floating point array `a`, `a[I]` is smaller, i.e. it's as long as the number of True values in `I`" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:03.316839Z", "start_time": "2019-01-17T22:54:03.310892Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a.shape = (9,)\n", "I.shape = (9,)\n", "Number of True values in I: 6\n", "Number of values in array a[I]: (6,)\n", "Number of values in array a[I]: 6\n" ] } ], "source": [ "print('a.shape = ', a.shape) # yields a tuple with the array dimensions here, just 9\n", "print('I.shape =', I.shape)\n", "print('Number of True values in I:', np.sum(I)) # True becomes 1 in computations\n", "print('Number of values in array a[I]: ', a[I].shape) # shape yields a tuple\n", "print('Number of values in array a[I]: ', len(a[I])) # len yiels an integer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's an example of elegant usage of logical indexing to select subareas." ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:03.979704Z", "start_time": "2019-01-17T22:54:03.669152Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAkUlEQVR4nO2de3xcZZ3wv79M0gu0SbiWFjANSwu2haSlFCIWUgsvrCAgWIVFoVuRXQQRQYSusuKLUFbRRUXZt8rVLZaLuBTWK7UDKINAbQKUcimml0BboJJpAzRpMs/7xzmTTqZzn3PmXOb3/Xzmk5lzfZ6Z9vmd312MMSiKoihKLmq8HoCiKIrif1RYKIqiKHlRYaEoiqLkRYWFoiiKkhcVFoqiKEpeVFgoiqIoeVFhoSgVQkRmi8grOfZPFBEjIrUlXn++iPyp9BEqSnZUWCihQkSiIvKuiIz0eizpGGOeNMYclvwsIutE5EQvxmJ/Txd6cW8lmKiwUEKDiEwEZgMGOD3PsZFKjElRwoIKCyVMnA88DdwFXJC6Q0TuEpHbROTXIvIeMEdEJojIL0XkbRHpEpHLUo6fJSIxEekRkU0icquIjMh0UxG5W0SutN8faJuSvmh/PlRE/i4W7SLSbW//OfAh4BER6RWRr6Vc8jwR2SAi74jI17NNVkT2EZFlIrJNRJ4B/iFt/0dE5FkRidt/P2JvvwFLqN5q3/tWe3z/KSJv2cc/LyLTCvrWlerAGKMvfYXiBawFvggcBewExqXsuwuIA8dhPSTtAawE/h0YARwC/A042T7+KOBYoBaYCKwBLs9y3wXAI/b7fwJeB+5L2few/b4d6E45bx1wYsrniVha0U+B0UAL0Ad8OMt9lwL3A3sC04A3gD/Z+/YG3gU+Z8/hXPvzPvb+KHBhyrVOtr+PRkCADwPjvf5N9eWfl2oWSigQkY8CTcD9xpiVWAv2P6Ud9rAx5s/GmARwBLCfMeb/GmP6jTF/w1qkzwEwxqw0xjxtjBkwxqwD/h9wQpbbPw7MFpEa4HjgO1hCCfucx4uczreMMR8YYzqBTiyhkT7fCHA28O/GmPeMMS8Cd6cccirwmjHm5/YcfgG8DHwiyz13AmOBwwExxqwxxmwqctxKiFFhoYSFC4DfG2PesT/fS5opCtiY8r4JmGCbmXpEpAf4N2AcgIhMFpFHRWSziGwDbgT2zXRjY8zrQC/QimXeeRR4U0QOozRhsTnl/fvAmAzH7IelMaTOaX3K+wlpn5P7D8x0Q2PMH4FbgR8DW0RksYjUFzluJcSosFACj4iMBj4NnGAv7puBrwAtIpL6VJ5aYnkj0GWMaUx5jTXGfNzefxvWk/gkY0w9liCRHMN4HPgUMMIY84b9+XxgL6AjyznllHx+GxgADk7Z9qGU929iCUTS9r+R7d7GmB8aY44CpgKTgavKGJ8SMlRYKGHgTGAQmIL1dN+KZXN/EmvBzsQzwDYRuVpERotIRESmicjR9v6xwDagV0QOBy7OM4bHgUuBJ+zPUeBLWD6EwSznbMHylRSNfc2HgOtEZA8RmcJwTerXwGQR+ScRqRWRz2B9P49mureIHC0ix4hIHfAesAPrO1UUQIWFEg4uAO40xmwwxmxOvrDMKudlSnKzF9tPYAmWLuAd4GdAg33IV7F8HtuxfBn35RnD41gCJiks/oTlRH8i6xmwCPiGbQb7at5Z7s6lWCaqzVgO/DuTO4wxW4HTgCuBrcDXgNNSzHQ/AD5l56T8EKjHmue7WOaqrcDNJYxJCSlijDY/UhRFUXKjmoWiKIqSFxUWiqIoSl5UWCiKoih5UWGhKIqi5KWkUshBYN999zUTJ04s6dz33nuPPffc09kBeURY5hKWeYDOxa+EZS7lzmPlypXvGGP2S98eWmExceJEnnvuuZLOjUajtLe3OzsgjwjLXMIyD9C5+JWwzKXceYhIeuY/oGYoRVEUpQBUWCiKoih5UWGhKIqi5EWFhaIoipIX14SFiNxhd916MWXb3iLyBxF5zf67V8q+hSKyVkReEZGTU7YfJSIv2Pt+KCK5Kn8qiqIoLuCmZnEXcEratmuA5caYScBy+zN2xcxzsEojnwL8JKVH8m3ARcAk+5V+TUVRFMVlXBMWxpgngL+nbT6DXd287sYqLZ3cvtQY02eM6cJqjzlLRMYD9caYmLEqHt6Tco6SgVgMFi2y/iqKojhFpfMsxiVbNRpjNonI/vb2A4GnU47rtrfttN+nb1cyEIvB3LnQ3w8jRsDy5V6PSHGS2MYY0XVR2ie203Zwm9fDUaoMvyTlZfJDmBzbM19E5CIskxXjxo0jGo2WNJje3t6Sz/WSJUs+RF9fM4mE0NeX4I471nHGGdnmshqrgVsrlvXP3wT1N8lEKXNZHV/Nlc9fyc7ETupq6vjekd9jaoP3v1u1/y5+xK15VFpYbBGR8bZWMR54y97ezfD2kAdhtYXstt+nb8+IMWYxsBhg5syZptQsxqBmco4cCUuWJDWLGhYsOIS+vg27zSUej9HZeRWJRD81NSNoaVlOQ4O/n1SD+ptkopS5xJ6MMWAGSJBgwAywbe9ttM/e/RqV1j6q/XfxI27No9Khs8vY1frxAuDhlO3niMhIEWnGcmQ/Y5ustovIsXYU1Pkp5yhptLVZpqfrr7f+tmVZK3p6oiQS/cAgiUQ/PT3RSg5TKYH2ie2MiIwgIhFGREbQPrF9t2NiG2PMvWcu1664lrn3zCW2UR1XinO4plmIyC+AdmBfEekGvgncBNwvIp8HNgDzAIwxq0XkfuAlrCb0l6T0Lb4YK7JqNPAb+6Vkoa0tu5BI0tjYTk3NiCHNorGxvSJjU0qn7eA2lp+/PKfWEF0XpX+wn0EzSP9gP9F10Yr6NuKxOD3RHhrbG2loa8h/ghIoXBMWxphzs+yam+X4G4AbMmx/Dpjm4NCqnoaGNlpaltPTE6Wxsd1xE1Q8HnPt2tVGulkp1+Kf1D76B/uzah9uEY/F6ZzbSaI/Qc2IGlqWt6jACBl+cXArFaahoc2Vhdzyh8wNlD/EryTNSsnFf/n5y3MKi0K0D7foifaQ6E/AICT6E/REe1RYhAwVFj4jFoNoFNrb85uT/Egmf4gKi9IoxayUT/solXwmpsb2RmpG1AxpFo3tjY6PoVjULOYsKiwcpNSFPnnePvvA5ZcPz5MImsBQf4hzeGlWSqUQE1NDWwMty1t8szirWcx5VFg4RKaEuEIW+tTzRCCRsF79/ZYA8UpYlOp3cNsfUk14aVZKpVATU0Nbg28WZDWLOY8KizRiMfj+9ydx331w/vmFL9bRqLXADw4Wt9CnnldTA5GIJTRGjLA0lGLH7oQJq1y/g1v+kGrELbNSMfjFxFSMWckvYw4TKixSiMWshba/fwIAd94JK1YUtvC2t1sLfFKzKHShTz/vlltg69bSTFmlaDaZ8IPfQSOq/EO5JiYnfAfFmpX8ZhYLAyosUohGYedOSFYZKUZDSCbEFftkX+p56ZSq2WTCa79DJs1G8ZZSTUxO+Q5KMSv5ySwWBlRYpNDeDnV10N9vlaUq1hRUSEKck+elUqpmkwmv/Q6ZM8xVuwgiTvkOvDAraTTVcFRYpNDWZj2RL1r0JgceeGBRPguvQ16d0lCSJP0O8XiM9esXOS40cpmZMms2fY7dW6kcTi3ylTYraTTV7qiwSKOtDa644jXa2wuvhO6kv6AcnNBQUnErwS7fdTNrNtGy76tUHicX+UqYlZLaxI4NOzSaKg0VFg7gpL/AT7jl6C7kuk5GVKmz3FsC4ztYDZ1XWdqERASpFQxGo6lsVFg4gJP+Aj/hlqO7kg50LT+iFEwHQ9qEwTD+C+MZ9aFR6rOwUWFRILl8Ek77C/yCW47uSjrQ/RAGHHZC4whuZZh/5YDzDwj2fBxGhUUBFOKTcNpf4CTlmGHcSrCrVOKe12HAoSfFdBN4R/BUNDcjByosCiDYPonVgeuK5yRehwGHng5C5QgOjH/FA1RYFEC5Pglvw2o7qsoMk0mL0vIjLtKK6/kPoTFzBRwVFgVQjk/C+7Da1qoxw6gz2wNcNt1ovoN/UGFRIKX6JLw3YU2tGjOMOrO9wU3TTdirxwZJa1Jh4TJ+CKsNkhmmHGe8OrPDR5irxwZNa1Jh4TJhDat1AydKo1eLFlUthLl6bCatKbndj3NVYVEB/BxW6yecMCNVUouKbYx53pioGghrhFK61lS3T52vNQ0VFkrJFGsyynd8kMxIsY0x5t4zd6jl6fLzlzsiMFQAVQ/pWpPf/TMqLJSSKNZkVMjxQTIjRddF6R/sZ9AM0j/YT3RdtOzF3S0BlIsgOVizEeQ5pGtNfvbPqLBQSqJYk1GhxwfFGd8+sZ0RkRFDC3v7xPayr+mGAMpF0BysmQjDHJL43T+jwkIpiWJNRkEyMRVC28FtLD9/uaMmIzcEUC5KMXv47Sne76abYvGzf0aFhVJSuGqxJqMgmZgKpe3gNkef/N0QQLkoNiw101O814Q5tNZvqLCocrL1uy5kUS/WZOSlicntnhZOOaadFkC5KNbskTHU02OZ73fTTZhQYVHlpPsSNm++hy1b7s7qiI7HY2zefA8ABxxwfiA0BLfLgHjhmHaKYsweGZ/ifdDttpg5+M2MVgpezUGFRZWT7ksAsjqi4/EYHR3tGNMPwObNd9LausL3AsPtMiCVdkx7Rcan+KjXoyqMeCzO5ns2s/nOzZgBE1hnuJcOfRUWVU66LwEYplmkOqJ7eqIYs3PoszH+rr+UND3V1e3jqnO90o5pL/GzAzYbQwvsjgQYa1tQneFeOvRVWCi7+RKyOaIbG9sRqRvSLET8G9WUbno69NBb2Llzqys+i0o7ppXiGFpgbUGBEFhnuJcOfU+EhYh8BbgQ6+d7AfhnYA/gPmAisA74tDHmXfv4hcDngUHgMmPM7yo1Vm97UXhDNkd0Q0Mbra3RQPgs0k1PO3dupalpoWv3q6RjWimO1AVWIsIBCw4IbMtULx36FRcWInIgcBkwxRjzgYjcD5wDTAGWG2NuEpFrgGuAq0Vkir1/KjABeExEJhtjBt0cZywG99wDd94JAwNe9aIonEoJtaAkzYUtr0PZRbEO3oIX2NWwPra+pEW4kk5nr0yBXpmhaoHRIrITS6N4E1gItNv778ZynV0NnAEsNcb0AV0ishaYBcTcGtzq1fVcdRXs2AHGVl393E7V+wZL/iOMeR1BwO1Fs1QHb74FNh6Lw5XQNdBVtOM4TFnkuaip9A2NMW8ANwMbgE1A3Bjze2CcMWaTfcwmYH/7lAOBjSmX6La3uUZHRyP9/bsEhYh3vSgKIVODJcUSGE1NC1VQVIjkotl1bRedczutBdhhspX1duK67KSk67o1Jr/hhRlqLyxtoRnoAR4Qkc/mOiXDNpNhGyJyEXARwLhx44iWuGoedlgttbVNGCNEIoZTTtnMySdvoa9vW1EL8erV9XR0NNLa2sPUqdtKGksh1NfXU1vbgjFCba2hvr6TaNS6X29vb8nfg58IyzwgxHNZgpV3kYBEX4JVd6xyPg+jHmvVMtbfrvouuqJdpV9vNdBhXdfUGmRAir9utjElr92KZUSvEG79+/LCDHUi0GWMeRtARB4CPgJsEZHxxphNIjIeeMs+vhs4OOX8g7DMVrthjFkMLAaYOXOmaS9ZFYiyYkUkxQdwIMUqM7EYXHVVZUxD7e0wY0aqz2LG0L5oNErp34N/CMs8ILxziY+M07kkxRyzwAVzTDvEZ+Q3dRViDovH4nRetWu8iUsTNO/bXLwJLcOY0q9dSdOUW/++vBAWG4BjRWQP4ANgLvAc8B5wAXCT/fdh+/hlwL0i8n0sB/ck4Bm3B1luw6JK997WBkuK11QqUqcQ/0MhPoR08xHboOl7TY6MKWwFDsEDYWGM+YuIPAj8FRgAVmFpA2OA+0Xk81gCZZ59/Go7Yuol+/hL3I6EcgI/9N5WlErjh6S9Qhfq9JyFRGvCsTGEscChJ9FQxphvAt9M29yHpWVkOv4G4Aa3x+Uk2ntbUbyh0IU6XRNa1bfKsTGEscChZnC7iJqGFKXyFLNQD9OEos6PIwxCIokKC0XJg/bFLp5y8y3KPT9sC7UfUGGhKDkIcvlxrygnSS0s1WHDSMWT8hQlE/F4jPXrFxGPu5aYXxKZyo8ruSk1SS0pZDb9v02YPhP6JLegoZqF4jn5mhNZAmQJ8fjIimdjV1P5cacoNRIoTNVhw4gKC5epxqq1xZKrOVFSkEAfnZ1LHO9yl49iyo974dso5Z5u128qNRIoTNVhy8WPHf1UWLiIFvgrjFwVYncJkoQrXe4KoZDy4174Nkq5Z6WK3pXiYPZzuGklF2+/FiZUn4WLuF3gLxaDRYusv0EmWSG2ufn63TSHpCCBGl+XGvfCt1HKPf1e9K6hrYGmhU2+WByTVKJAYip+/Y1Us3ARN7O4w6a15Gq41NKynFWr7qClZYFvK8h64dso5Z5hzCx2m0qX7vDrb6TCwkXczOKudO0pL7EERJ9vBQV401q1lHv62dTjVyq9ePv1N1Jh4TJuZXFr7Sn/4UVr1VLuGfaENaf9C14s3n78jVRYBBStPeUc8XhMO+q5RKWjetxwDvsxMskLVFgEGK09VT75cjwyUWi4arWXCfEiqsdp/4JfI5O8QIWFEkic0gZy5XhkotBwVTdDaYMihLzo6eC0fyGMfSlKRYWFEjhK0QaykSvHIxOZwlUzLdiFHpekGG0lKLWqvIjqcdq/4NfIJC9QYaEEjmK1gVwkQ3ML1VIKDVctJqy1GAFQrBDyEq+iepx0Dvs1MskLVFgogSIej7FjxwZEajEGRxL1suV4ZKLQcNViwlqLEQBBq1Xlx6ieYgnDHJxAhYUSGFLNTyIRxo//AgcccL4vy38Uc1wxAsCLfA5FARUWvkSLD2Ym1fxkDIwa9aFQhLoWKwC8yOdwEg1FDSYqLHyGl2U8/J5vUKwzOkgEQQA4schrKGpwUWHhM7wq4+FkhJFbFOuMVpzDqUVeQ1GDi1ad9RnJMh6RSGXLeGSKMPIjDQ1tNDUtLFtQ+LUzn19xqhJqMhSViDY2ChqqWfgMr8p4hNnEk04QtKhicTtRz6l8Aw1FDS4qLHyIF2U8gm7iKcbf4mSehh+oRKKek4u8hqIGExUWVUauRbWYfAM/UaymEDYtqtA8jXId1LrIVzcqLFzCj+GvYTS/QPGaQtC1qHQKydPQKKTC0dDezKiwcAG/drELm/klSSmaQlC1qEwUkqfhRjXWnmgP1APtJV/Gd6hQzY4KCxfwaxe7sJlfkoRNUyiFfHka6Q7qun3qWL9ofUlPz6kLKrUQnxEPzYKqob3ZUWHhAn7tYhfmRTVMmoIbpDqo6/apY+3la0t+ek5dUDGEakHVKrPZUWHhAn7uYqeLavWSdFCvX7S+rKfn1AWVWkK1oGpob3ZUWLiEdrFT/Eq5T8+pC2pXfVfoFlQ3o76C7Dz3RFiISCPwM2AaYIAFwCvAfcBEYB3waWPMu/bxC4HPYym+lxljflfxQStKSHDi6Tm5oHZFu1wYYTgJuvPcq3IfPwB+a4w5HGgB1gDXAMuNMZOA5fZnRGQKcA4wFTgF+ImIRDwZtU0sBosWWX8VJYg0tDXQtLApUItV0HGqZIpXVFyzEJF64HhgPoAxph/oF5Ez2BWEdzcQBa4GzgCWGmP6gC4RWQvMAjxZqv0aFhtU/F7pVnGWIJthyiXoznMvzFCHAG8Dd4pIC7AS+DIwzhizCcAYs0lE9rePPxB4OuX8bnubJ/g1LDaIeJEkqMLJO4JuhimGTEIx6M5zL4RFLTAD+JIx5i8i8gNsk1MWJMM2k/FAkYuAiwDGjRtHNBotaYC9vb1Zz62vr6e2tgVjhNpaQ319J9HotpLuUwl2zWU10AG0Yln0/MASoA9IkEj0sWrVHfbn3cn1mxTOauBKYCdQB3wPL74LZ+biD4qay66fm0RfglV3rMr2c3uCY79Lvn9mbVjzduBWmXDr35cXwqIb6DbG/MX+/CCWsNgiIuNtrWI88FbK8QennH8Q8GamCxtjFgOLAWbOnGnaS0xwiEajZDu3vR1mzEgNi51R0j0qRTQaZfr0kXR2XuW7Mh/x+Eg6O5ekjGtB1nHl+k0KZf36GF1dA0ACGKC5eRtNTeVdsxScmItfKGYu8ZFxOpekaBYL/KVZOPW7rI+tp2ugK/nPjOZtzTS1N5V93UJx699XxYWFMWaziGwUkcOMMa8Ac4GX7NcFwE3234ftU5YB94rI94EJwCTgmUqPO5WghcX6tcxHpZMEw5rBHhSCboYplKD7JrLhVZ7Fl4AlIjIC+Bvwz1iRWfeLyOeBDcA8AGPMahG5H0uYDACXGGMGvRl2MPHzIlnJJMEwZ7AHhWqoXBtWoeiJsDDGdAAzM+yam+X4G4Ab3BxTmNFFchfFCCe3GwpVE9UWBZUUivFYvOQaXH5DM7irBC3zURyVaChULVRTFFQqYZu39uBWlAxkaiiklIZfk9GST/3xWNyV6/t13qWimoWiZKCQhkJKYfjR4Zvpqd9p/DjvcsgqLERkWQHn/90YM9+54QQTP3bF8xPxeIzNm+8B4IADzg+EOayQhkJhxWn/gh8dvhmf+h3+if0473LIpVl8GLgwx34BfuzscIKHlv/ITTweo6OjHauqC2zefCetrSsCIzCqSUiAe3Z2v0VBZXzqLzBBsBhh6rd5l0MuYfF1Y8zjuU4WkW85PJ7AEcTyH06XvMh1vZ6eKMbsHPpsjH/yPJTdqZZOcRmf+qP5zwub07oYsgoLY8z9+U4u5Jiw49eueNlwuh5Tvus1NrYjUjekWYj4K88jFa0bFT47ey5KeeqvFmGaibwObhGZCXwdaLKPF8AYY450eWyBwM9d8TLhdDZ3vus1NLTR2hr1vc/Ci6KGfiRsdnanqSZhmk4h0VBLgKuAF7CqnShpBKn8h9PZ3IVcLwg5Hn4tieIFYbKzO001C9NChMXbxphCIqOUIvEiisrpbO6wZIf7uSSK4i/yCdOwZqsXIiy+KSI/w+peNxQvYIx5yLVRVQHFRFE5LVScftIPguaQD6eFnpYKqU7C7AAvRFj8M3A4VmX2pBnKACosyqDQKCoNza0cTgk9LRVSPGF5Gg+zA7wQYdFijDnC9ZFUGYVGUQUxNLfayVQqRIVFdsL0NB5mB3ghwuJpEZlijHnJ9dFUEfmiqJKmp3322V2o+DljXMNPw1EqpJJP+mF6Gg+zA7wQYfFR4AIR6cLyWWjorENki6JKNz3dcgts3bpL+/CrWUrDTy2ylQpJ+jHq4/W00+7tIFNIFwyVftIP29N4WKPJChEWp7g+CmUY6aanrVth4UJr36JF/jVLafjpLtJLhaT6MWqllhkzZuwmRLxwhmcSDJV+0k8+jW++Z7Nr91DKJ6+wMMasr8RAqpVMJqVc/gw/Z4xr+Gl2Uv0YxpghP4bXzvBMgsGrJ/0td28h0Z9gy91bXNdm3DCzhcVJn41cVWf/aoyZkevkQo5RspMt0imXP8PPGeNhyblwg1Q/Rq3UDvkxvHaGZxIMXtjdK6nNuGFmC5OTPhs5q86KyPM59gsQrm+jwuSKdMqVFe7njPEw5Fy4Qaofo/7v9UMCwWtneDbBUGm7eyW1GTcEU5ic9NnIJSwOL+D8QacGUo342aSkOE/SjxGNRodt87pvhh8csuVqM8WYgNwQTGFz0mciV9VZ9VW4jJ9NSkrlqMa+GZkoVWgVawJyw8wW5pDZJNpW1WP8bFJSlCBQignIDW3KDxqam9R4PQBFUZRySJqAiBBaE5AfKKSfxaXAEmPMuxUYj6IoSlF4YQIKe5hsJgoxQx0APCsifwXuAH5njDHuDktRFKVwKmkCcjpMNiiCJ68ZyhjzDWAScDswH3hNRG4UkX9weWxKGrGYlcEdi3k9EsXPxGNx1i9aTzwW98V1wkYmH0khZPo+k4Kn69ouOud2+vq7LsjBbYwxIrIZ2AwMAHsBD4rIH4wxX3NzgIqF26XKtQBgOHDqqbcaksxKpZQw2WzfZ5DyMwrxWVwGXAC8A/wMuMoYs1NEaoDXABUWFcDNUuVaADA8OLX4BGkRqzSl+EiyfZ+pgkciwo4NO3hz8Zvs3LpzWGFHP5ipCtEs9gXOSs+7MMYkROQ0d4alpONmAp8WAAwPmZ56S1lsqiHJrByK9ZFk+z5TiyhuvnMzmxZvYlNiE9RAzcgaDr3lUNZevtYXGl4hhQT/Pce+Nc4OR0klvcigWwl8WgCwePzaNjX9qRcoyZxUDUlmTlGIMM71fSbNUWbA7OpFmrA0kLd/+bZvNDxNyvMpuYoMOo0WACyOXJVi/SBEUp961y9aX/JiE/YkMycoxreT6/sc0jz6EpbAqLFyRvY7ez/iT8Z9oeF5JixEJAI8B7xhjDlNRPYG7gMmAuuATydzO0RkIfB5rFpUlxljfufJoCtIpdupZioAqE7vzGSrFOt1ufF04rE4OzbsQGoFg/F8sSkFv9jrs+GUbydV86jbp26Yz2LPI/b0xXfgpWbxZWANUG9/vgZYboy5SUSusT9fLSJTgHOAqcAE4DERmWyMCWQRw0JbonpdZFCd3tnJVinWyXLj5WooqU+8EhHGf2E8B5x/gC8X3GwEISLLSd9ONs3DLxqeJ8JCRA4CTgVuAK6wN58BQ70m7waiwNX29qXGmD6gS0TWArOAwGUbFBP+6nWRQXV6ZydbpVinyo0XqqHkEiipT7wGw6gPjcq74PjtKT4IEVnV5NvxSrO4BSvkdmzKtnHGmE0AxphNIrK/vf1A4OmU47rtbYGjWNOSl0UG1emdm0yVYp0qN16IhpJPoBT7xOvHp/igRGT55cnfbSouLOxw27eMMStFpL2QUzJsy1huREQuAi4CGDdu3LC+AcXQ29tb8rm5qK+vp7a2BWOE2lpDfX0n0eg2x++TSnlz+S7QQSLRyqpVfVjKnje49ZtkYnV8NR3xDlobWpnaMLXo89too+/1PqKvRzPuzzeX+ng9tVKLMYZaqaX+7/W7Hb9kwxL6BvpIkKBvoI87VtxB34f6hl/I+vlItCZY1bcq98+3BOjDisLpS7DqjlXW5zy4/rsUM4cyqeS/MTdxbR7GmIq+gEVY2sE6rIzw94H/Bl4BxtvHjAdesd8vBBamnP87oC3ffY466ihTKitWrCj53Hw89ZQxN95o/a0Ebs6lkhQzj56ep8y6dTeanp7iv+SnNjxlRn97tIl8K2JGf3u0eWqD8z9UIXN5asNT5sYnbsx6f6fH2fNUj3l89ONmRWSFeXz046bnqZ6CzgvLvy9jwjGXnqd6zIoLVxT8+2UCeM5kWFMrrlkYYxbaAgBbs/iqMeazIvJdrEzxm+y/D9unLAPuFZHvYzm4JwHPVHjYjlFt/SsqHVFVrmPe657YSfI1RHK6w1412d7DStKUSB90Lul03JTopzyLm4D7ReTzwAZgHoAxZrWI3A+8hFWX6hIT0EioasOLiKpyHfNe98QuhkI67BXjtK4W23tYGQoISLgTEOCpsDDGRLGtkMaYrcDcLMfdgBU5pQQILyKqynXM+6EntlP40WmtuEdqYp8bAQF+0iyUkOFFRJUT2ehh6YkdhNBTxTmSpsRVd6yiZYHzDwYqLBTX8KqMSKZs9GokKKGninM0tDVAH648FKiwUFzFTwt3tZUvUae14iQqLJSqoFrLl6jTWnGKvG1VFSUMZHK2K0pQ8aLlrWoWSlWg5UuUsOBVlJsKC6Uq0J4dSljwKspNhYVSNRTibPdD8yI38VtlWaV4vIpyU2GhKDZ+a17kNJqkFw68inJTB7ei2GSqC1UIsY0xFj25iNhGf7dYyWS+UIJJQ1sDTQubKirsVbNQFJtS6kIFSRvRJL3gEo/F2XzPZgDPOh6qsFAUm1LqQvmlSm0haJJeduKxOCyB+Mi4776XeCxOR3sHpt9q47P5zs20rmit+DhVWChKCsXWhQpSlVrQJL1MuF3au1x6oj2Ynbv6vZl+40mdLxUWAScW865PtxKuKrXVitulvculsb0RqZMhzUJGCI3tjRWPbFNhUQZeL9SxGMyda/XzHjECli93dxzVVlupUMJSpbZacbu0d7k0tDXQGm0d5rMAKh7ZVlXCYufOnXR3d7Njx46cxzU0NLBmzZqcx/T1wbZtMH269bejA0aOdHCwBSACDz00/HP6sAuZSyEkEn30928DpvP3v29jxIgOampGMmrUKA466CDq6urKvkc1Eva8jiDgdmlvJ0g3H65ftL7iiXlVJSy6u7sZO3YsEydORESyHrd9+3bGjh2b81qbNsHAwK7P48bB+PFOjbQwenvh1VchkYCaGpg8GcaMGX5MIXMphL6+TfT375rwiBHjGDHiALZu3Up3dzfNzc1l36PaCFIkVdhxs7S3G3gR2VZVwmLHjh15BUWhjB1rLdDJhdqB9bhoxoyxBMT27db90wWFk0QiY7HSchJADZHIWESEffbZh7ffftu9G4eYIEVSKf7Ci8i2qhIWgCOCAiq7UOcbRyXuXVs7htGjJzM4uJ1IZCy1tdZNnfo+q5GgRVIp/qLSkW1VJyycpFILtV+orR0zJCSU8tFIquBRzbW1VFgoiof4JZKqmhfBXKR+L1D5CCQ/ocIiD16HxyqK22iBwcykfy/jLhjnSWnwfGOslJDXQoI5SOYxXHut9TdWZp249957j1NPPZWWlhamTZvGfffdx8qVKznhhBM46qijOPnkk9m0aRMAP/3pTzn66KNpaWnh7LPP5v3333dgRkq5BKVoYDGUVWBwNRXv2FYp0r8XgJoRNRDBF/kYSWHWdW0XnXM7Xf8NVFjkIBq1Et4GB62/0Wh51/vtb3/LhAkT6Ozs5MUXX+SUU07hS1/6Eg8++CArV65kwYIFfP3rXwfgrLPO4tlnn6Wzs5MPf/jD3H777WXPRymPZKjrtSuuZe49c0MjMJJhmMUugvFYHK6kYotVpUn/Xg44/wBalrfQfH2zL7SvSlcRVjNUDtrbrczo/n6oq4NZs8q73hFHHMFXv/pVrr76ak477TT22msvXnzxRU466SQABgcHGW8na7z44ot84xvfoKenh97eXk4++eQyZ6OUS1hDXUsNw+yJ9sBOfFsmo1yyfS9+mWOlcy1UWOSgrQ0eeQR+9SuYMQP22stKhCs1Amry5MmsXLmSX//61yxcuJCTTjqJqVOnEstg35o/fz7/8z//Q0tLC3fddRfRctUapWzCHOpaShhmY3sj1AED/jDLuIGfCy9WOtdChUUepkyBBvs3SCSsvIpShcWbb77J3nvvzWc/+1nGjBnD4sWLefvtt4nFYrS1tbFz505effVVpk6dyvbt2xk/fjw7d+5kyZIlHHjggcOu1dvrfY5HtaGhrsNpaGuA70HztmaNovKISgozFRZ5cDJT+4UXXuCqq66ipqaGuro6brvtNmpra7nsssuIx+MMDAxw+eWXM3XqVK6//nqOOeYYmpqaOOKII9i+ffvQdQop86G4g19CXX3DVGhqbxq2yakInbCH8wZtfios8uBkpvbJJ5+c0ffwxBNP7Lbt4osv5uKLL854ne3bLUEB5Ws7iuIkToXhFnKdoC22qQQxXFmjoQpgzBirSKBfFuSktgPe1aVSlEyUEqETj8V3C7/Nd51Kh406TRD7oatmEUD8UpdK8Rd+KHdebIROtifsfNfJtNjmezL3kyYSxH7oFRcWInIwcA9wAFYJ08XGmB+IyN7AfcBEYB3waWPMu/Y5C4HPA4PAZcaY31V63H7DL3WpBgZ6GRiIE4/HtCGSh/il3HmxETrZFv1813FKKHlFEPuhe6FZDABXGmP+KiJjgZUi8gdgPrDcGHOTiFwDXANcLSJTgHOAqcAE4DERmWyMGfRg7EoKAwO9fPDBqwwM9NDZeRYtLctVYHiEn3JAionQybXo57qOU0LJS/wclpuJigsLY8wmYJP9fruIrAEOBM4A2u3D7gaiwNX29qXGmD6gS0TWArOAcKTPBpjBwe1YyiEkEv309ERVWHhEUHNAynnCdkooKYXhqc9CRCYC04G/AONsQYIxZpOI7G8fdiDwdMpp3fY2xWN2NUSCmpoRNDa2ezoeN/GDPyAXQc4BqcQTdhDNPn5DjDHe3FhkDPA4cIMx5iER6THGNKbsf9cYs5eI/BiIGWP+295+O/BrY8wvM1zzIuAigHHjxh21dOnSYfsbGho49NBD845tcHCQSCRS+uQK5MYbb2TMmDFcdtllrt3D/bl8wNq1rxKPd2JZCt2ht7eXMR45aVbHV3Pl81eyM7GTupo6vnfk95jaUNxcV8dX0xHvoLWhlaZIE2PGjBm2rdjr+QUvfxenCctcyp3HnDlzVhpjZqZv90SzEJE64JfAEmPMQ/bmLSIy3tYqxgNv2du7gYNTTj8IeDPTdY0xi4HFADNnzjTt7e3D9q9Zs6agftSpfavj8Rg9PVEaG9uHmVicyKAeOXIkI0eOdKRHdjYK6cE9MDBAbW2p/xTGMmrUVqZPvyTrd+UE0WiU9N+zUsSejDFgBkiQYMAMsG3vbbTPLnwssY0xrrrnqiET0XenfZcZ/zBj2LZUh7TftZhUvPxdnCYsc3FrHhXPsxCrD+ftwBpjzPdTdi0DLrDfXwA8nLL9HBEZKSLNwCTgmUqMNR6P0dk5l66ua+nsnEs8brlJkhnUb7xh/e3tLfyaN9xwA4cddhgnnngir7zyytD29vZ2nnvuOQDeeecdJk6cCMBdd93FWWedxSmnnMKkSZP42te+NnTO73//e9ra2pgxYwbz5s2jN8NA7rrrroylzufPn88VV1zBnDlzuPrqq3n99dc55ZRTOOqoo5g9ezYvv/wyAI888gjHHHMM06dP58QTT2TLli1FfVdhIOkPiEikJH9AuvO5I96R0SEN4a1s6yWZ8jiqcQzl4kVS3nHA54CPiUiH/fo4cBNwkoi8Bpxkf8YYsxq4H3gJ+C1wSaUioXp6oiQS/cDgkAMXMmdQF8LKlStZunQpq1at4qGHHuLZZ58t6LyOjg7uu+8+XnjhBe677z42btzIO++8w7e//W0ee+wx/vrXvzJz5ky+//3v73buJz7xiaylzl999VUee+wxvve973HRRRfxox/9iJUrV3LzzTfzxS9+EYCPfvSjPP3006xatYpzzjmH73znO0V9V2Eg6Q+4fs71JYWkpgub1obWrAIomxBRSsMPyXt+GIMTeBEN9SdAsuyem+WcG4AbXBtUFhob26mpGUEi0T/MgVtqvagnn3yST37yk+yxxx4AnH766QWdN3fuXBrsaoZTpkxh/fr19PT08NJLL3HccccB0N/fT1uGVn5r1qzhc5/7XMZS5/PmzSMSidDb28tTTz3FvHnzhvb19fUB0N3dzWc+8xk2bdpEf38/zc3NGceY7bsql6RpC+rZFSxXecqpCZXufO57vS+rQzqoUU1+xQ8hs34YgxNoBncOGhraaGlZvpsdvpwMassKtzu1tbUkbHVlx44dw/aNHDly6H0kEmFgYABjDCeddBK/+MUvct7v4osv5uGHH85Y6nzPPfcEIJFI0NjYSEdHx27nf+lLX+KKK67g9NNPJxqNct1112W8T7bvqhySpi1LY6klHp9BQ0Obq74Rt0gVNtHXo7ttSz0uXYgEyYfhN/wQMuuHMTiBCos8NDS0ZVyQSsmgPv7445k/fz7XXHMNAwMDPPLII/zLv/wLABMnTmTlypXMmjWLBx98MO+1jj32WC655BLWrl3LoYceyvvvv093dzeTJ08edly+UucA9fX1NDc388ADDzBv3jyMMTz//PO0tLQQj8eHzrn77rtzjinbd1UqqaYtMEOmraQAqakZEcpEwFQh4pfM7KDih5DZSowhtZSJW2ghwQoyefIMPv7xz3Dkka2cffbZzJ49e2jfV7/6VW677TY+8pGP8M477+S91n777cddd93Fueeey5FHHsmxxx475JRO5Rvf+AbHHHMMJ510EocffnjW6y1ZsoTbb7+dlpYWpk6dysMPW/EF1113HfPmzWP27Nnsu+++Jcy6dJKmLYgAdTQ2tofaN5IJ9WGUT0NbA00Lmzwv7+HWGNJ9Iqx2/BaAh3kWbjNz5kyTjC5KsmbNGj784Q/nPbeQcNNi8aoHhRtzSafQ77UUkianrq562tsvGWaaCqpmUUxoo981i7CEm0Jw57J+0Xq6ru2yFPAI8M/Q/tP2kq8nIv7Js6hGtAdFaSRNW11d0aHPTvpG/O7/CHJmtlIZ0n0iidaEK/dRYVEhnOy4V+045RsJipai3fmql0LKqqf7RFb1rXJlLCosKoT2oNgdr5/qM/k//CgslOqkmLLqw+prRd0ZjwqLClKpHhSppUj8SqaneqCiwsOt3BBFcQK/5WeosAgZ6Y70gw6q8aXQSH+q37z5HrZsubuiJiE3ckMUxSn8lp+hwiIgFFq4MN2R/v77/vyJ05/qgbwmITfMVk7nhihKPgpt7+qHHJFU/LmSKMMoJux27Fj41rfm89GPnsZJJ32KPfYYAEZmPthD0p/qgWGaxe4modV0dl7le2e0ouSi2Paufuqmp8IiD34otVBM2O2YMVBfD3vvbQkVY9wJo3OC9Kf63CahDnVGK4HHb36IYtAM7hy4US563bp1HH744Vx44YVMmzaN8847j8cee4zjjjuOSZMm8cwzz/Dee++xYMECjj76aKZPn87jjz9MTQ28+eY6vvCF2fzjP85gxowZPPXUUwAYY7j00kuZMmUKp556Ku+++xZ77WUJjmj0d7S2TmPatKksWLCAvr4+nnnmGc466ywAHn74YUaPHk1/fz87duzgkEMOAeCHP/whU6ZM4cgjj+Scc84pe96F0NDQRlPTwixCoHUom1ud0dVL0Et9J/0QRPCFH6IYVLPIQaZSC05oF2vXruWBBx5g8eLFHH300dx777386U9/YtmyZdx4441MmTKFj33sY9xxxx309PQwa9YsnnzyRBoa9ucPf/gD++47itdee41zzz2X5557jl/96le88sorvPDCC2zZsoUpU6awYMECenvf4eKLL2XZsh8zaVIzF198M7fddhuXXnopq1ZZsdhPPvkk06ZN49lnn2VgYIBjjjkGgJtuuomuri5GjhxJT09P2XMun6mhcUavjq8m9mRMk+yKpFgTjh/xmx+iGFRY5MCtctHNzc0cccQRAEydOpW5c+ciIhxxxBGsW7eO7u5uli1bxs033wxYVWj//vcNTJgwgUsvvZSOjg4ikQivvvoqAE888QTnnnsukUiECRMm8LGPfQyANWs6aGqawKRJTUCC8847k5/+dCmXX345hx56KGvWrOGZZ57hiiuu4IknnmBwcHCoXtWRRx7Jeeedx5lnnsmZZ57pyLzLJQzO6MUrF/Plzi9jMIyMjPRd+Q4/E2QTTip+8kMUgwqLHLhVaiG15HhNTc3Q55qaGgYGBohEIvzyl7/ksMMOG3beddddx7hx4+js7CSRSDBq1KihfZlKn4vskfKphkhk1+fZs2fzm9/8hrq6Ok488UTmz5/P4ODgkID63//9X5544gmWLVvG9ddfz+rVq8tovTocr5PxvCK2McYlv76EQbt3V99gn2PaajXgt1BSv5GMsnKr9Yv6LPLQdnAbC2cvrOh/6JNPPpkf/ehHJIs8Jk1G8Xic8ePHU1NTw89//nMGB61F5/jjj2fp0qUMDg6yadMmVqxYAcC0aTPYsOEtNm7sZ/Toydx774OccMIJQ+fccssttLW1sd9++7F161Zefvllpk6dSiKRYOPGjcyZM4fvfOc7Q42TnCDM7VfzEV0XHepZAhCRSKCbG8VjcVhCxfwHSRNO8/XNgTRBuUlq5VmudOc3Uc3Ch1x77bVcfvnlHHnkkRhjmDhxIo8++ihf/OIXOfvss3nggQeYM2fOUPOiT37yk/zxj3/kiCOOYPLkyUMCYdSoUfzkJ7dx3nmXMDAwwNFHH82//uu/AnDMMcewZcsWjj/+eMAyO+2///6ICAMDA3z2s58lHo9jjOErX/kKjY2NjsytmktstE9sZ2TtSHYM7CBSE+HWj98aWK0iuTjRB51LOotevAvNNUgnqCYct0k10WFwxUSnJcozkKusd6HJcX7BbyXKSy3eF9Ty0enENsa4Y8UdLJizILCCAnYvi918fTNNC5sKOtevjuog/xtL/U6phekrppf8nWqJcgfwqidFmKj2EhttB7fR96G+QAsKSPEf9BXvPwiLo9pPpEZZddV3ufJ9qrAoAu1J4QxhiGoqFD8kdbpBcnFadccqWhYUpxmoo9odkia6rmiXK9dXYVEE2pPCf2SLrPJDxJWbXe78IIQa2hqgj6KfYoOca5BOqb6XIKLCogjc7Enhd1/IwEAvg4PbiUTGUlvrjwFm83/4pamRW0mdfm+1WghhcFT71ffiFho6WyRjxsD48c4LildfhTfesP46FKXqGAMDvXzwwav097/BBx+8ysCAPwaYKbIq1/ZKk0zqjEjE0aTOTEIoiAS9dEcm30uYUc3CB/jdFzI4uB1I5gckGBzc7gvtIlvzIr80NXIrqdOtygKVJAxP5dXme1Fh4QMy+ULmz5/Paaedxqc+9Smvh0ckMhZLCU1gZYJXxlmT6nfIRDKyavPmezJuz+ezqIRfw43+2W4JoUriVkRUJX0IYfK9FIIKizxU4h9fui9k1KgBV+5TKrW1Yxg9enJFfRbpfgf4LtlqGCT7YGzZcveQfyJfxJVf/Bql4oYQqiRuPJV7oa2EwfdSKOqzyEFqCn3n3E5HbKvr1q1jypRpbNpk+SZuvvlmrrvuOk47rZ0f/ejfOPXUE/jBD34AwGOPPcbs2bOZPHkyjz766ND5s2fPZsaM4WXKkwlFn/rUpzj88MM577zzcDLhsrZ2DCNHjq+Y+Snd7wAdBR1XqH/CL36NasWN0h3V5kOoNKpZ5MANVfm996C/33Jm19RAX1/K/Xp6ePzxxwHLDLVu3Toef/xxXn/9debMmcPatWvZf3+rTPmoUcPLlINVQ2r16tVMmDCB4447jj//+c+0tLSUNd5SKdfEk+53SCRaCzquUP+EX/wa1YwTT+Wpmn+1+RAqjQqLHLjxj++99yD5wJ9IWMKixtbvPvOZzww79tOf/jQ1NTVMmjSJQw45hJdffpnm5uaMZcoBZs2axUEHHQRAa2sr69at80RYOGHiSfc7rFrVV9Bxhd6n2jPJC8XPeQSZzE7V5EOoNCoscuCGA6uhoXao1aklJHYM7UsWBkySXnZcRPjP//zPrGXKU0ufRyIRBgby+z7cyJ9wqljgcL9DtMDjSr1+ZvyQ/OYVfo9YyqT5Ny1s8tUYw0RgfBYicoqIvCIia0Xkmkrdt6Gtoah/gL29DPkjMnHIIeOIx99i9OitNDX18fvfP5r1Wg888ACJRILXX3+dv/3tbxx22GFZy5SXglv5E0kTT9BboLrRVrecsSx6clFFx+B3H0CQW5QGkUBoFiISAX4MnAR0A8+KyDJjzEvejmw4hRQarKur45vf/HfOOOMYmpubOfzww7Ne77DDDuOEE05gy5Yt/Nd//RejRo3KWqa8FNzKnwiLicetDOxi8Spj2+8+gGoLXfWaQAgLYBaw1hjzNwARWQqcAfhKWBSaXHfZZZdx2WWX5bzWXXfdlXH7pEmTeP7554c+L1q0CID29vZh5ZVvvfVWe0zbs97DzfyJMBQL9Evym1dCKwiLcTWFrnpNUITFgcDGlM/dwDEejSUrQSs06EX+RJDwS/Kbl0JLF2MlSSCaH4nIPOBkY8yF9ufPAbOMMV9KO+4i4CKAcePGHbV06dJh12loaODQQw/Ne7/BwUEikUhJY/3ggxref7+WPfYYYPToRP4TXKacuRTK2rVricfdre/T29vLGD/VQCmDUuayOr6ajngHrQ2tTG2Y6tLIiqfafxc/Uu485syZE+jmR93AwSmfDwLeTD/IGLMYWAxWp7z0rldr1qxhzJgxu0UZpVNOd7ldp43MdVjFcLtTnjGGUaNGMX36dNfuAcHuYpZOKXNpz5K97jXV/rv4EbfmEZRoqGeBSSLSLCIjgHOAZcVeZNSoUWzdutXRzOZqxhjD1q1bh4XvKooSTgKhWRhjBkTkUuB3QAS4wxizutjrHHTQQXR3d/P222/nPG7Hjh2hWQDdnsuoUaOGEgEVRQkvgRAWAMaYXwO/LucadXV1NDc35z0uGo26blapFGGai6Io3hEUM5SiKIriISosFEVRlLyosFAURVHyEog8i1IQkbeB9SWevi/wjoPD8ZKwzCUs8wCdi18Jy1zKnUeTMWa/9I2hFRblICLPZUpKCSJhmUtY5gE6F78Slrm4NQ81QymKoih5UWGhKIqi5EWFRWYWez0ABwnLXMIyD9C5+JWwzMWVeajPQlEURcmLahaKoihKXlRYKIqiKHmpamEhIt8VkZdF5HkR+ZWINKbsW2j3+35FRE5O2X6UiLxg7/uh5Kt37hFe9SwvFRE5WERWiMgaEVktIl+2t+8tIn8Qkdfsv3ulnJPxN/IDIhIRkVUi8qj9OajzaBSRB+3/J2tEpC3Ac/mK/W/rRRH5hYiMCspcROQOEXlLRF5M2Vb02Mtav4wxVfsC/g9Qa7//D+A/7PdTgE6sphTNwOtAxN73DNAGCPAb4B+9nkeGeUXsMR8CjLDnMsXrceUZ83hghv1+LPCq/Tt8B7jG3n5NIb+RH17AFcC9wKP256DO427gQvv9CKAxiHPB6rbZBYy2P98PzA/KXIDjgRnAiynbih57OetXVWsWxpjfG2MG7I9PYzVVAqu/91JjTJ8xpgtYC8wSkfFAvTEmZqxv/h7gzEqPuwCGepYbY/qBZM9y32KM2WSM+av9fjuwBus/+BlYCxb23zPt9xl/o4oOOgsichBwKvCzlM1BnEc91iJ1O4Axpt8Y00MA52JTC4wWkVpgD6wGaoGYizHmCeDvaZuLGnu561dVC4s0FmBJWsjc8/tA+9WdYbvfyDb+QCAiE4HpwF+AccaYTWAJFGB/+zA/z/EW4GtAal/dIM7jEOBt4E7bpPYzEdmTAM7FGPMGcDOwAdgExI0xvyeAc0mh2LGXtX6FXliIyGO2jTL9dUbKMV8HBoAlyU0ZLmVybPcbQRnnbojIGOCXwOXGmG25Ds2wzfM5ishpwFvGmJWFnpJhm+fzsKnFMn3cZoyZDryHZe7Ihm/nYtvzz8Ayy0wA9hSRz+Y6JcM2X8ylAFxZvwLT/KhUjDEn5tovIhcApwFzbdUMsvf87maXqSp1u98oqGe53xCROixBscQY85C9eYuIjDfGbLLV6Lfs7X6d43HA6SLycWAUUC8i/03w5gHW2LqNMX+xPz+IJSyCOJcTgS5jzNsAIvIQ8BGCOZckxY69rPUr9JpFLkTkFOBq4HRjzPspu5YB54jISBFpBiYBz9iq3nYROdaOIjgfeLjiA8+PIz3LK4n9fd4OrDHGfD9l1zLgAvv9Bez6vjP+RpUabzaMMQuNMQcZYyZife9/NMZ8loDNA8AYsxnYKCKH2ZvmAi8RwLlgmZ+OFZE97H9rc7H8YkGcS5Kixl72+uWVd98PLyzHz0agw379V8q+r2NFEbxCSsQAMBN40d53K3YWvN9ewMexIopeB77u9XgKGO9HsVTi51N+j48D+wDLgdfsv3vn+4388gLa2RUNFch5AK3Ac/bv8j/AXgGey7eAl+3/vz/HihYKxFyAX2D5WnZiaQifL2Xs5axfWu5DURRFyUtVm6EURVGUwlBhoSiKouRFhYWiKIqSFxUWiqIoSl5UWCiKoih5UWGhKIqi5EWFhaK4jIhMFJEPRKSjyPM+Y5eSftSloSlKwaiwUJTK8LoxprWYE4wx9wEXujMcRSkOFRaKUgYicrRYzbNGiciednOdaXnOmWg3E/qZXdRyiYicKCJ/thvZ+Kmst6IAVVBIUFHcxBjzrIgsA74NjAb+2xjzYp7TAA4F5gEXYdXy+ieskienA/+GP/ukKFWMCgtFKZ//i7Xg7wAuK/CcLmPMCwAishpYbowxIvICMNGVUSpKGagZSlHKZ29gDFY72FEFntOX8j6R8jmBPsQpPkSFhaKUz2LgWqzmWf/h8VgUxRX0CUZRykBEzgcGjDH3ikgEeEpEPmaM+aPXY1MUJ9ES5YriMnZP8UeNMTmjpLKc2w581RhzmsPDUpSiUDOUorjPINBQSlIe8BPgXTcGpSjFoJqFoiiKkhfVLBRFUZS8qLBQFEVR8qLCQlEURcmLCgtFURQlL/8fcbQERKLj0gMAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(-200, 1000, 241) #; print(x)\n", "y = np.random.rand(len(x)) * 1000 # randon number as many as there are x-values\n", "\n", "# logical indexing, specifying subarea, names are intuitive\n", "sea = x < 0\n", "dunes = np.logical_and(x > 0, x <=300)\n", "grass = np.logical_and(x > 300, x < 500)\n", "urban = x >= 500\n", "\n", "plt.title('Area with dots')\n", "plt.xlabel('x [m]')\n", "plt.ylabel('y [m]')\n", "plt.grid()\n", "\n", "plt.plot(x[sea], y[sea], 'b.', label='sea') # b. means blue as dots\n", "plt.plot(x[dunes], y[dunes], 'y.', label='dune area') # y. means yellow as dots\n", "plt.plot(x[grass], y[grass], 'g.', label = 'meadows') # g. means green as dots\n", "plt.plot(x[urban], y[urban], 'm.' , label= 'urban') # m. means magenta as dots\n", "\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## List comprehensions to generate and filter lists" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with an `iterable`, here an array, but it could by any `iterable` like a `tuple` or an `list`. An iterable is any object that can be iterated over (list, tuple, array, dictionary, even a tring)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:04.312566Z", "start_time": "2019-01-17T22:54:04.308864Z" } }, "outputs": [], "source": [ "a = np.array([-3, 2, -1, 1.2, 3.2, 0.7, 3.2, -3, 5.1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By way of trivial start, generate a `list` from the `items` in this `iterable`" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:04.664246Z", "start_time": "2019-01-17T22:54:04.659271Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b = [-3.0, 2.0, -1.0, 1.2, 3.2, 0.7, 3.2, -3.0, 5.1]\n" ] } ], "source": [ "b = [y for y in a] # the most basic list comprehension, b will be the same as a\n", "print('b = ', b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that `b` is not an array, it is now a `list`. Use `np.array(b)` to make it an `array` if you need it.\n", "\n", "Next, select only the items from the array `a` that are larger than some value. This is called `filtering` the `iterable`:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:05.014818Z", "start_time": "2019-01-17T22:54:05.009522Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b = [2.0, 1.2, 3.2, 3.2, 5.1]\n" ] } ], "source": [ "b = [x for x in a if x > 1.1] # same but filtering out some values according to condition\n", "print('b = ', b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And of course, you now got a list which is a bit shorter.\n", "\n", "While get the items one by one in the comprehension, you may be tempted to do something with them before putting them in the new list. Here we take the logarithm of each item and we only take the items that fulfill the if-criterion of the comprehension." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:05.384364Z", "start_time": "2019-01-17T22:54:05.378797Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.3862943611198906, 0.36464311358790924, 2.3263016196113617, -0.7133498878774649, 2.3263016196113617, 3.25848107946056]\n" ] } ], "source": [ "c = [np.log(x**2) for x in a if x > 0]\n", "print(c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this example it becomes clear that comprehensions can be really advanced and are very flexible.\n", "\n", "As another example, we take only the items of type `str` and then convert these strings into their uppercase variants before putting them in the list:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:05.720606Z", "start_time": "2019-01-17T22:54:05.713112Z" } }, "outputs": [ { "data": { "text/plain": [ "['AHA', 'JOHN']" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = ['aha', 1, np.array([3, 4.2]), 'John'] # list with objects of different type\n", "[x.upper() for x in a if isinstance(x, str)] # get the strings from the list of items in uppercase" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clearly, comprehensions are extremely flexible for generating new lists from iterables (an iterable can be a list, tuple, array, string, i.e. is somethign you can iterate over like in `for this in that`. Then `that` is the iterable and `this` is an item from the iterable.\n", "\n", "You can always turn the obtained list into an array if it contains only numbers:" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:06.051327Z", "start_time": "2019-01-17T22:54:06.046627Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.34002931 0.60527907 0.19226882 0.34173995 0.05403698 0.90028315\n", " 0.78349603 0.86749492 0.54073142 0.90522101]\n" ] } ], "source": [ "a = np.random.rand(10) # get 10 random numbers from the uniform distribuion\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "ExecuteTime": { "end_time": "2019-01-17T22:54:06.121350Z", "start_time": "2019-01-17T22:54:06.113483Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b = [0.340029308382628, 0.6052790719936957, 0.3417399455487897, 0.9002831456929633, 0.783496028777145, 0.8674949162157715, 0.5407314164689259, 0.9052210060737772]\n", "\n", "c = [0.34002931 0.60527907 0.34173995 0.90028315 0.78349603 0.86749492\n", " 0.54073142 0.90522101]\n", "\n", "d = ['y = 0.34', 'y = 0.19', 'y = 0.34', 'y = 0.05']\n", "joined into a single string: y = 0.34; y = 0.19; y = 0.34; y = 0.05\n" ] } ], "source": [ "b = [y for y in a if y > 0.2] # filter the values > 0.2 from the list into a new list\n", "print('b = ', b) # list\n", "print()\n", "c = np.array([y for y in a if y > 0.2]) # turn the list inte an array\n", "print('c = ', c) # an array\n", "print()\n", "\n", "# generate a list of strings. Using the floating point values obtained from the \n", "d = ['y = {:.2f}'.format(y) for y in a if y < 0.5] # an array of strings\n", "print('d = ', d)\n", "\n", "# If you have a list of strings, you can join them into a single string\n", "print('joined into a single string: ', '; '.join(d)) # this is now a string" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* We showed strings, tuples, lists and numerical arrays (numpy ndarrays), which are all iterables (meaning: they can be iterated over).\n", "* We showed some basic logical indexing.\n", "* And used them in loops.\n", "* We looped using several iterables in parallel to get values of them that belong together like the x, y and Q of successive wells.\n", "* We have only used one-dimensional arrays here. That is, we only scratched the surface of the possibilities multi-dimensional arrays. \n", "* We did not use list of lists or list of tuples or tuples of lists either to keep things simple. \n", "* We only used basic examples of list-comprehensions. List comprehensions are extremely powerful to generate lists on the fly.\n", "* We did not look at dicts and dict comprehensions (generating dictionaries, just as powerful as list comprehensions).\n", "* We only showed the simplist of logical indexing, which is a extremely powerful way of indexing in arrays.\n", "* There is still a world of plotting possibilities to be explored (see matplotlib gallery for examples ready to use).\n", "* But it's a good start.\n", "\n", "You mway want to google for \"Exploratory computing with Python by Mark Bakker\" to get more excercises focused on engineering students like you.\n", " " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.7 64-bit ('miniconda3')", "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.9.7" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "nbTranslate": { "displayLangs": [ "*" ], "hotkey": "alt-t", "langInMainMenu": true, "sourceLang": "en", "targetLang": "fr", "useGoogleTranslate": true }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false }, "vscode": { "interpreter": { "hash": "1ae27910742444e13e8c72a4e4606f2a7bceb52340d5c783cbc8f639475be36a" } } }, "nbformat": 4, "nbformat_minor": 2 }