{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Foundations of Computational Economics #14\n", "\n", "by Fedor Iskhakov, ANU\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Vectors and matrixes (Numpy)\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "[https://youtu.be/QLp3PEziRJE](https://youtu.be/QLp3PEziRJE)\n", "\n", "Description: NumPy arrays data types and differences to native Python, operations on the arrays, solving linear systems of equations." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Scientific stack in Python\n", "\n", "Collection of modules (libraries) used in scientific numerical computations:\n", "\n", "- **``NumPy``** is widely-used scientific computing package for implements fast array processing — vectorization \n", "- **``SciPy``** is a collection of functions that perform common scientific operations (optimization, root finding, interpolation, numerical integration, etc.) \n", "- **``Pandas``** is data manipulation package with special data types and methods \n", "- **``Numba``** is just in time (JIT) compiler for a subset of Python and NumPy functions \n", "- **``Matplotlib``** serves for making figures and plots " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### NumPy library\n", "\n", "\n", "\n", "- **Vectorization in Python** \n", "- **NumPy** is a widely-used scientific computing package for brings fast array processing to Python \n", "- Reference guide: [https://docs.scipy.org/doc/numpy/reference/](https://docs.scipy.org/doc/numpy/reference/) \n", "- Runs fast compiled code written in C & Fortran under the hood " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Importing modules in Python\n", "\n", "- `import LIBRARY as ref`, then call library functions as `ref.function` \n", "- `from LIBRARIY import function` or `from LIBRARIY import *`, then call library functions directly \n", "- keeping conventional reference is a good idea for making your code understood by others! " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# import libraries\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Power of NumPy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.28 ms ± 63.6 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)\n", "575 µs ± 26 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)\n", "NumPy is on avarage 9.193 faster\n" ] } ], "source": [ "N = 1000000\n", "data_list = list(range(N)) # Native Python list\n", "t1 = %timeit -n10 -r10 -o mean1 = sum(data_list)/N\n", "\n", "import numpy as np\n", "data_array = np.array(range(N)) #NumPy array\n", "t2 = %timeit -n10 -r10 -o mean2 = data_array.mean()\n", "\n", "print('NumPy is on avarage %2.3f faster' % (t1.average/t2.average))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "SI orders of magnitude https://en.wikipedia.org/wiki/Micro-" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Scientific computing: more advanced treatment of numbers\n", "\n", "Inherited from low lever C implementation\n", "\n", "- int8, uint8 (signed and unsigned) \n", "- int16, uint16 \n", "- int32, uint32 \n", "- float16 \n", "- float32 \n", "- float64 \n", "\n", "\n", "Full list of types:\n", "[https://numpy.org/doc/stable/user/basics.types.html](https://numpy.org/doc/stable/user/basics.types.html)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Array initialization with type" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x , takes 99 bytes\n", "y , takes 102 bytes\n", "z , takes 120 bytes\n" ] } ], "source": [ "import sys\n", "x = np.array([-1,0,1.4],dtype='bool')\n", "y = np.array([-1,0,1.4],dtype='int16')\n", "z = np.array([-1,0,1.4],dtype='float64')\n", "print('x %s, takes %d bytes' % (type(x[0]),sys.getsizeof(x)))\n", "print('y %s, takes %d bytes' % (type(y[0]),sys.getsizeof(y)))\n", "print('z %s, takes %d bytes' % (type(z[0]),sys.getsizeof(z)))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "NumPy array hold data of **the same type**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Integer overflow" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[255 0 1]\n" ] } ], "source": [ "x = np.array([0,0,0],dtype='uint8')\n", "x[0] = 255\n", "x[1] = x[0] + 1\n", "x[2] = x[1] + 1\n", "print(x)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Infinity and not-a-number" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-inf nan inf inf]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in true_divide\n", " after removing the cwd from sys.path.\n", "/usr/local/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in true_divide\n", " after removing the cwd from sys.path.\n" ] } ], "source": [ "# Inf and NaN\n", "# np.seterr(all=None, divide='ignore', over=None, under=None, invalid='ignore')\n", "x = np.array([-1,0,1,10],dtype='float64')\n", "print( x / 0 )\n", "# y = 10 / 0 # core Python" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True False True\n" ] } ], "source": [ "# Comparing nans and infs\n", "a = (np.inf == np.inf)\n", "b = (np.nan == np.nan) # Can not compare nan to nan!\n", "c = np.isnan(np.nan)\n", "print (a, b, c)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### What is inside array?\n", "\n", "First, import module `obj_explore.py`: this is not trivial because jupyter notebooks require imported modules to be on the known `PATH`" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# add path to the modules directory\n", "import sys\n", "sys.path.insert(1, './_static/include/')\n", "# import the obj_explore() function\n", "from obj_explore import *" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[, , , , ]\n" ] }, { "data": { "text/plain": [ "array([1., 2., 3., 4., 5.])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.array([1,2,3,4,5],dtype='float64')\n", "# a = np.arange(5,dtype='uint8') + 1\n", "print([type(a_element) for a_element in a])\n", "a" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 2. 3. 4. 5.]\n", "------------------------------------------------------------\n", "Object report on object = array([1., 2., 3., 4., 5.])\n", "Objec class : \n", "Parent classes : \n", "Occupied memory : 136 bytes\n", "PUBLIC PROPERTIES\n", "T = array([1., 2., 3., 4., 5.]) \n", "base = None \n", "ctypes = \n", "data = \n", "dtype = dtype('float64') \n", "flags = C_CONTIGUOUS : True\n", " F_CONTIGUOUS : True\n", " OWNDATA : True\n", " WRITEABLE : True\n", " ALIGNED : True\n", " WRITEBACKIFCOPY : False\n", " UPDATEIFCOPY : False\n", " \n", "flat = \n", "imag = array([0., 0., 0., 0., 0.]) \n", "itemsize = 8 \n", "nbytes = 40 \n", "ndim = 1 \n", "real = array([1., 2., 3., 4., 5.]) \n", "shape = (5,) \n", "size = 5 \n", "strides = (8,) \n", "PUBLIC METHODS\n", "all \n", "any \n", "argmax \n", "argmin \n", "argpartition \n", "argsort \n", "astype \n", "byteswap \n", "choose \n", "clip \n", "compress \n", "conj \n", "conjugate \n", "copy \n", "cumprod \n", "cumsum \n", "diagonal \n", "dot \n", "dump \n", "dumps \n", "fill \n", "flatten \n", "getfield \n", "item \n", "itemset \n", "max \n", "mean \n", "min \n", "newbyteorder \n", "nonzero \n", "partition \n", "prod \n", "ptp \n", "put \n", "ravel \n", "repeat \n", "reshape \n", "resize \n", "round \n", "searchsorted \n", "setfield \n", "setflags \n", "sort \n", "squeeze \n", "std \n", "sum \n", "swapaxes \n", "take \n", "tobytes \n", "tofile \n", "tolist \n", "tostring \n", "trace \n", "transpose \n", "var \n", "view \n" ] } ], "source": [ "obj_explore(a,'public')\n", "# help(a)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Memory footprint" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "import sys\n", "def memory_usage(var,grow,steps=10):\n", " \"\"\"Returns data on memory usage when var is grown using supplied function for given number of steps\"\"\"\n", " d={\"x\":[],\"y\":[],\"v\":[]} # dictionary for x, y data, and values\n", " for i in range(steps):\n", " var=grow(var) # next value\n", " d[\"v\"].append(var)\n", " d[\"x\"].append(i+1)\n", " d[\"y\"].append(sys.getsizeof(var))\n", " return d" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3gV1bn48e+bKyQkXAKEcBOQcL+JEVC8gCgCVVGPoD2e46Va2vOzVautYmur1lpta731nKOHHi162qJovQAVlaqobRUEIdxRQISQEBJACAm57vf3x6wdNjEhOyE7s5O8n+fZz55Zs2bmTdB5s9bMrCWqijHGGNMYMX4HYIwxpuWyJGKMMabRLIkYY4xpNEsixhhjGs2SiDHGmEaL8zuASOjatav269fP7zCMMabF2FlYzP4vtxSqareG7Ncqk0i/fv1YtWqV32EYY0yLcKC4nHEP/o39D33jy4bua91ZxhjTxi3dkEdloHHvDFoSMcaYNm5xdi4DuiU3at+IJhER+YGIbBSRDSKyQETaiUh/EVkhIp+LyIsikuDqJrr1bW57v5Dj3O3Kt4rIRZGM2Rhj2pL8w6Ws+OIAl47u2aj9I3ZPRER6AbcAw1T1qIgsBK4GZgCPqeoLIvI0cCPwlPs+qKoDReRq4FfAVSIyzO03HOgJ/E1EBqlqVUPiqaioICcnh9LS0ib7GVuLdu3a0bt3b+Lj4/0OxRjTzJasy0MVLh7Vkx80Yv9I31iPA9qLSAWQBOQB5wP/6rY/B9yHl0RmumWAl4H/FBFx5S+oahnwhYhsA8YBHzUkkJycHFJSUujXrx/eYQ2AqrJ//35ycnLo37+/3+EYY5rZ4uxchmWkMrB7h0btH7HuLFXdAzwC7MJLHoeA1cBXqlrpquUAvdxyL2C327fS1U8LLa9ln2oiMkdEVonIqoKCgq/FU1paSlpamiWQGkSEtLQ0a6EZ0wbtPlDC2t1fcemYxnVlQQSTiIh0xmtF9MfrhkoGptdSNfhIQG1Xdz1B+fEFqvNUNUtVs7p1q/0xZ0sgtbPfizFt06LsXAAuHpXR6GNE8sb6BcAXqlqgqhXAK8BZQCcRCXaj9QZy3XIO0AfAbe8IHAgtr2UfY4wxjbQ4O5fTT+lM785JjT5GJJPILmCCiCS5extTgE3Ae8CVrs51wOtueZFbx21/V73JThYBV7unt/oDmcDKCMYdMbGxsYwZM4YRI0Ywa9YsSkpK6qy7c+dO/vznP1evz58/n+9973vNEaYxpg34PL+ILXuLuOQkWiEQ2XsiK/BukH8KrHfnmgfcBdzubpCnAc+4XZ4B0lz57cBcd5yNwEK8BPQmcHNDn8yKFu3bt2ft2rVs2LCBhIQEnn766Trr1kwixhjTlBZn5xIjMCNakwiAqt6rqkNUdYSq/ruqlqnqDlUdp6oDVXWWe+oKVS116wPd9h0hx3lQVU9V1cGqujSSMTeXc845h23btvHTn/6UJ554orr8Jz/5CU8++SRz587lww8/ZMyYMTz22GMA5ObmMm3aNDIzM7nzzjur91mwYAEjR45kxIgR3HXXXdXlHTp04Cc/+QmjR49mwoQJ5OfnN98PaIyJWqrK4nV5nHlqGt1T2p3UsVrl2Fn1uX/xRjblHm7SYw7rmcq9lwwPq25lZSVLly5l2rRpTJ8+nSuuuIJbb72VQCDACy+8wMqVKxk1ahSPPPIIS5YsAbzurLVr17JmzRoSExMZPHgw3//+94mNjeWuu+5i9erVdO7cmalTp/Laa69x2WWXUVxczIQJE3jwwQe58847+f3vf88999zTpD+3Mabl2bDnMF8UFvOdcwec9LHaZBLxy9GjRxkzZgzgtURuvPFGEhISSEtLY82aNeTn53PaaaeRlpZW6/5TpkyhY8eOAAwbNowvv/yS/fv3M2nSJIJPpF1zzTV88MEHXHbZZSQkJHDxxRcDcPrpp7Ns2bJm+CmNMdFu8bpc4mOFaSN6nPSx2mQSCbfF0NSC90Rquummm5g/fz579+7lW9/6Vp37JyYmVi/HxsZSWVmJ9+xB7eLj46sf3w3WN8a0bYGAsiQ7l3Mzu9EpKeGkj2cDMEaByy+/nDfffJNPPvmEiy7yhgZLSUmhqKio3n3Hjx/P+++/T2FhIVVVVSxYsIDzzjsv0iEbY1qo1bsOknuolEsaOVZWTW2yJRJtEhISmDx5Mp06dSI2NhaAUaNGERcXx+jRo7n++uvp3LlzrftmZGTw0EMPMXnyZFSVGTNmMHPmzOYM3xjTgizOzqVdfAwXDktvkuPJibpDWqqsrCytOSnV5s2bGTp0qE8RnVggEGDs2LG89NJLZGZm+hJDNP9+jDFNo7IqwISH3mF8/zT+65qxX9suIqtVNashx7TuLJ9t2rSJgQMHMmXKFN8SiDGmbfhox34Kj5RzyeiTezcklHVn+WzYsGHs2LGj/orGGHOSFmfn0iExjkmDuzfZMa0lYowxbUBZZRVLN+xl6vB02sXHNtlxLYkYY0wb8MFnhRSVVjZ6BsO6WBIxxpg2YHF2Lp2T4pk4sGuTHteSiDHGtHIl5ZUs25TP9JEZxMc27WXfkkgzEhHuuOOO6vVHHnmE++67z7+AjDFtwjub93G0oqrJu7LAkkizSkxM5JVXXqGwsNDvUIwxbcii7FzSUxM5o1+XJj+2JZFmFBcXx5w5c6qHdg91/fXX8/LLL1evd+jQAYDly5dz3nnnMXv2bAYNGsTcuXP505/+xLhx4xg5ciTbt2+v3v+73/0u55xzDoMGDaoe/fecc845bryuiRMnsm7dukj+mMaYKHLoaAXvby3gGyN7EhvT9FNht833RJbOhb3rm/aYPUbC9IfrrXbzzTczatSo4+YDqU92djabN2+mS5cuDBgwgJtuuomVK1fyxBNP8Lvf/Y7HH38c8Cayev/999m+fTuTJ09m27Zt1YM7Pv7443z22WeUlZUxatSoRv+YxpiW5e2NeymvCnDpmKbvygJriTS71NRUrr32Wp588smw9znjjDPIyMggMTGRU089lalTpwIwcuRIdu7cWV1v9uzZxMTEkJmZyYABA9iyZQuzZs1iyZIlVFRU8Oyzz3L99dc38U9kjIlmi7Jz6dslidG9O0bk+BFriYjIYODFkKIBwM+A5115P2AnMFtVD7p52J8AZgAlwPWq+qk71nVAcDalX6jqcycVXBgthki67bbbGDt2LDfccEN1WVxcHIFAAPBmHSsvL6/eFjoEfExMTPV6TEzMccO7B4d9D11PSkriwgsv5PXXX2fhwoXUHFPMGNN6FR4p45/b9/Pd8wZ87frQVCI5x/pWVR2jqmOA0/ESw6t4c6e/o6qZwDtuHWA6kOk+c4CnAESkC3AvMB4YB9wrIrUPadtCdOnShdmzZ/PMM89Ul/Xr14/Vq1cD8Prrr1NRUdHg47700ksEAgG2b9/Ojh07GDx4MODNV3LLLbdwxhln0KVL099YM8ZEp6Ub9lIV0CYb9r02zdWdNQXYrqpfAjOBYEviOeAytzwTeF49HwOdRCQDuAhYpqoHVPUgsAyY1kxxR8wdd9xx3FNa3/72t3n//fcZN24cK1asIDk5ucHHHDx4MOeddx7Tp0/n6aefpl07b+7k008/ndTU1ONaPsaY1m/x2lwyu3dgcHpKxM7RXDfWrwYWuOV0Vc0DUNU8EQmOBNYL2B2yT44rq6u8xTly5Ej1cnp6OiUlJcetf/zxx9XrDz30EACTJk1i0qRJ1eXLly+vXq65beLEibU++ZWbm0sgEKi+l2KMaf3yDh1l5c4D3HHhoIh1ZUEztEREJAG4FHipvqq1lOkJymueZ46IrBKRVQUFBQ0PtJV6/vnnGT9+PA8++CAxMfYchTFtxZLsPAAujmBXFjRPS2Q68Kmq5rv1fBHJcK2QDGCfK88B+oTs1xvIdeWTapQvr3kSVZ0HzANvUqqm/AFagvnz59dafu2113Lttdc2bzDGGN8tXpfLyF4d6d+14V3jDdEcf5p+k2NdWQCLgOvc8nXA6yHl14pnAnDIdXu9BUwVkc7uhvpUV9ZgrXEWx6ZgvxdjWpedhcWsyzkUkWFOaopoS0REkoALge+EFD8MLBSRG4FdwCxX/gbe473b8J7kugFAVQ+IyAPAJ67ez1X1QENjadeuHfv37yctLS2i/YMtjaqyf//+6pvwxpiWb3F2LgDfGNV0MxjWJaJJRFVLgLQaZfvxntaqWVeBm+s4zrPAsycTS+/evcnJycHul3xdu3bt6N27t99hGGOayOJ1uYzr14WendpH/FxtZtiT+Ph4+vfv73cYxhgTUVv3FvFZ/hEemDm8Wc5nj+sYY0wrsih7D7ExwvSRke/KAksixhjTaqgqi7PzOOvUNLp2SKx/hyZgScQYY1qJ7JxD7DpQEtFhTmqyJGKMMa3E4uxcEmJjuGh4j2Y7pyURY4xpBQIBZcm6XM4b3I2O7eOb7byWRIwxphVYufMA+YfLmrUrCyyJGGNMq7A4O5f28bFcMLR7/ZWbkCURY4xp4SqqAizdsJcLhqWTlNC8r/9ZEjHGmBbuH9sKOVBcziXNMMxJTZZEjDGmhVucnUdKuzjOG9yt2c9tScQYY1qw0ooq3t64l2nDe5AYF9vs57ckYowxLdjyrQUUlVU2+1NZQZZEjDGmBVu8Lpe05ATOOjWt/soRYEnEGGNaqOKySt7ZnM+MkRnExfpzObckYowxLdTfNudTWhHg0jH+dGWBJRFjjGmxFmfnktGxHaf37exbDPUmERFJFpEYtzxIRC4VkeYbmMUYY8zXfFVSzvufFXDxqAxiYvyb8juclsgHQDsR6QW8gzf3+fxwDi4inUTkZRHZIiKbReRMEekiIstE5HP33dnVFRF5UkS2icg6ERkbcpzrXP3PReS6hv+YxhjTury1cS8VVcqlo3v5Gkc4SUTcXOlXAL9T1cuBYWEe/wngTVUdAowGNgNzgXdUNRMvKc11dacDme4zB3gKQES6APcC44FxwL3BxGOMMW3V4uw8+qUlMaJXqq9xhJVERORM4Brgr66s3sFZRCQVOBd4BkBVy1X1K2Am8Jyr9hxwmVueCTyvno+BTiKSAVwELFPVA6p6EFgGTAvrpzPGmFZoX1Ep/9xeyCWjeyLiX1cWhJdEbgPuBl5V1Y0iMgB4L4z9BgAFwB9EZI2I/K+IJAPpqpoH4L6DQ072AnaH7J/jyuoqP46IzBGRVSKyqqCgIIzwjDGmZVq6fi8BhUt9esEwVL1JRFXfV9VLgf906ztU9ZYwjh0HjAWeUtXTgGKOdV3VprZ0qicorxnnPFXNUtWsbt2af/wYY4xpLouycxnSI4XM9BS/Qwnr6awzRWQT3v0MRGS0iPx3GMfOAXJUdYVbfxkvqeS7birc976Q+n1C9u8N5J6g3Bhj2pycgyWs/vKgb8Oc1BROd9bjePcl9gOoajbevY4TUtW9wG4RGeyKpgCbgEVA8Amr64DX3fIi4Fr3lNYE4JDr7noLmCoind0N9amuzBhj2py/rssD4JJR0ZFEwpq9RFV317h5UxXm8b8P/ElEEoAdeI8HxwALReRGYBcwy9V9A5gBbANKXF1U9YCIPAB84ur9XFUPhHl+Y4xpVRZl5zK6Tyf6piX5HQoQXhLZLSJnAeqSwS24rq36qOpaIKuWTVNqqavAzXUc51ng2XDOaYwxrdX2giNszD3MTy8O9y2LyAunO+u7eBf3Xnj3J8YA/y+SQRljjPm6Jdl5iMDFPsxgWJdwWiKDVfWa0AIRmQj8IzIhGWOMqUlVWZS9h/H9u5Ce2s7vcKqF0xL5XZhlxhhjImRzXhHbC4qj5qmsoDpbIu4t9bOAbiJye8imVKD552A0xpg2bFF2LnExwvQR0dOVBSfuzkoAOrg6oW+0HAaujGRQxhhjjlFVFmfncnZmV7okJ/gdznHqTCKq+j7wvoi8oqrrmzEmY4wxIdbs/oo9Xx3l9gsH+R3K14RzT+QpEVkpIv9PRDpFPCJjjDHHWbQ2l4S4GKYOT/c7lK8JZ+yss4F/wxt6ZJWI/FlELox4ZMYYY6gKKH9dn8f5g7uT0i765gMMa3pcVf0MuAe4CzgPeNJNNHVFJIMzxpi2bsUX+ykoKou6p7KCwhmAcZSIPIb3lvr5wCWqOtQtPxbh+Iwxpk1bnJ1LckIs5w/pXn9lH4TzsuF/Ar8HfqyqR4OFqporIvdELDJjjGnjyisDLN2wlwuHpdM+ITrfrKg3iajquW7MrCEiosBWVS132/4v0gEaY0xb9fdtBXxVUhG1XVkQ3jS3M4D/AbbjTRDVX0S+o6pLIx2cMca0ZYuz8+jYPp5zMqN3or1wurMeBSar6jYAETkVb651SyLGGBMhpRVVvL1xL5eM7klCXFjPQPkinMj2BROIs4NjsxEaY4yJgHe37KO4vCqqu7LgxGNnBR/f3SgibwAL8eY2n8WxCaKMMcZEwOLsXLp2SGTCgDS/QzmhE3VnXRKynI/3fghAAdA5YhEZY0wbV1Rawbtb9vHNcX2JjZH6d/DRicbOuuFkDy4iO4EivOl0K1U1S0S6AC8C/YCdwGxVPSje/LtP4E2RWwJcr6qfuuNch/eyI8AvVPW5k43NGGOi1bJN+ZRVBqK+KwvCfGP9JE1W1TGqGpwmdy7wjqpmAu+4dYDpQKb7zAGeAnBJ515gPDAOuFdErCVkjGm1Fmfn0qtTe8b2jf7hCv245T8TCLYkngMuCyl/Xj0fA51EJAO4CFimqgdU9SCwDJjW3EEbY0xzOFhczoefF3Lx6Ay8DproFukkosDbIrJaROa4snRVzQNw38F3+XsBu0P2zXFldZUfR0TmiMgqEVlVUFDQxD+GMcY0j6Ub9lIZUC5tAV1ZEN57IscRkZnAXlVdEUb1iW54lO7AMhHZcqJD11KmJyg/vkB1HjAPICsr62vbjTGmJVicncuAbskMy0j1O5SwNKYlMh64R0TqfdlQVXPd9z7gVbx7Gvmumwr3HXznJAdvuPmg3kDuCcqNMaZVyT9cysdf7OeSUT1bRFcWNCKJqOqPVfUSVZ1+onoikiwiKcFlYCqwAVgEXOeqXQe87pYXAdeKZwJwyHV3vQVMFZHO7ob6VFdmjDGtyl/X5aFKi3gqKyicsbNWAX8A/uxubIcrHXjVZdM4t/+bIvIJsFBEbgR24b28CPAG3uO92/Ae8b0BQFUPiMgDHHvB8eeqeqABcRhjTIuweF0uwzJSGdi9g9+hhC2ceyJX413QPwlJKG+r6gnvO6jqDmB0LeX7gSm1lCtwcx3HehZ4NoxYjTGmRdp9oIQ1u75i7vQhfofSIOFMj7tNVX8CDAL+jHcx3yUi97t3OIwxxpykxeu8W73fGJnhcyQNE9Y9EREZBfwW+A3wF+BK4DDwbuRCM8aYtmPR2lzG9u1Eny5JfofSIOHcE1kNfAU8A8xV1TK3aYWITIxkcMYY0xZ8nl/Elr1F3HfJML9DabATJhERiQH+oqq/rG27ql5RW7kxxpjwLV6XR4zAjFEtqysL6unOUtUANsSIMcZEjKqyODuXCQPS6J7Szu9wGiyceyLLROSHItJHRLoEPxGPzBhj2oCNuYf5orC4xQxzUlM4j/h+y32HPn6rwICmD8cYY9qWxdm5xMcK00b08DuURqk3iahq/+YIxBhj2ppAwOvKOjezG52SEvwOp1Hq7c4SkXgRuUVEXnaf74lIfHMEZ4wxrdmnuw6Se6i0RQ1zUlM43VlPAfHAf7v1f3dlN0UqKGOMaQsWZeeSGBfDBcPS/Q6l0cJJImeoaujwJe+KSHakAjLGmLagsirAG+vzuGBoOh0SGzwrR9QI5+msKhE5NbgiIgPw5kw3xhjTSB/vOEDhkXIuGd3y3g0JFU76+xHwnojswJsg6hTcCLvGGGMaZ1H2HjokxjFpcPf6K0excJ7OekdEMoHBeElkS8jQJ8YYYxqorLKKNzfsZerwdNrFx/odzkkJZ+ysWOAioJ+rP0VEUNVHIxybMca0Sh9+Vsjh0soW/VRWUDjdWYuBUmA9EIhsOMYY0/q9np1L56R4zh7Y1e9QTlo4SaS3qo6KeCTGGNMGfLxjP0vW5XLjxP7ExzZ4hvKoE85PsFREpjb2BCISKyJrRGSJW+8vIitE5HMReVFEElx5olvf5rb3CznG3a58q4hc1NhYjDHGT0WlFdyxMJu+XZL4wYWD/A6nSYSTRD7Gmyv9qIgcFpEiETncgHPcCmwOWf8V8JiqZgIHgRtd+Y3AQVUdCDzm6iEiw/Cm6B2ON6Lwf7v7NMYY06I8sGQTeYeO8ujsMSS34HdDQoWTRH4LnAkkqWqqqqaoamo4BxeR3sA3gP916wKcD7zsqjwHXOaWZ7p13PYprv5M4AVVLVPVL4BtwLhwzm+MMdHi7Y17Wbgqh/+YdCqnn9LZ73CaTDhJ5HNgg6pqI47/OHAnx27IpwFfqWqlW88BernlXsBuALf9kKtfXV7LPsYYE/UKj5Rx9yvrGZaRyq1TWkc3VlA47ak8YLmILAWq3w+p7xFfEbkY2Keqq0VkUrC4lqpaz7YT7RN6vjnAHIC+ffueKDRjjGk2qsrdr6ynqLSSP397DAlxLf9meqhwksgX7pPgPuGaCFwqIjOAdkAqXsukk4jEudZGbyDX1c8B+gA5IhIHdAQOhJQHhe5TTVXnAfMAsrKyGtNqMsaYJvfy6hyWbcrnJzOGMrhHit/hNLlw3li/H0BEklW1ONwDq+rdwN1u30nAD1X1GhF5CbgSeAG4Dnjd7bLIrX/ktr+rqioii4A/i8ijQE8gE1gZbhzGGOOX3QdKuH/xJsb378KNZ7fOqZnCmU/kTBHZhHvCSkRGi8h/17PbidwF3C4i2/DueTzjyp8B0lz57cBcAFXdCCwENgFvAjerqg0AaYyJaoGA8sOXvAHPH5k1mpiY2nrmW75wurMexxv2ZBGAqmaLyLkNOYmqLgeWu+Ud1PJ0laqWArPq2P9B4MGGnNMYY/z0zN+/YMUXB/j1laPo0yXJ73AiJqw7PKq6u0aRtQSMMaYOW/cW8Zu3tnLhsHRmnd7b73AiKpyWyG4ROQtQ93b5LRz/8qAxxhinvDLAD15cS2r7OB66YiTe626tVzgtke8CN+O9m5EDjHHrxhhjanjinc/YlHeYX14+kq4dEv0OJ+LCeTqrELimGWIxxpgWbfWXB3hq+XZmnd6bqcN7+B1OswhnPpH+wPc5Np8IAKp6aeTCMsaYlqW4rJLbF2bTs1N7fnbJML/DaTbh3BN5De/x28XYfCLGGFOrX76xmV0HSljw7QmktIv3O5xmE04SKVXVJyMeiTHGtFDvbd3Hn1bs4tvn9GfCgDS/w2lW4SSRJ0TkXuBtjh8769OIRWWMMS3EweJy7nx5HYPTU7hj6mC/w2l24SSRkcC/4w3hHuzOUrdujDFtlqpyz2sb+KqknPk3nEG7+LY31VE4SeRyYICqlkc6GGOMaUkWZefy1/V5/OiiwQzv2dHvcHwRznsi2UCnSAdijDEtSd6ho/z0tQ2M7duJ75w7wO9wfBNOSyQd2CIin3D8PRF7xNcY0yYFAsqPXlpHZUB5dPYY4mJb1xwhDRFOErk34lEYY0wL8vxHO/n7tkIevHwE/bom+x2Or8J5Y/395gjEGGNagm37jvDQ0i1MGtyNfx1ns6i23TaYMcY0UEVVgDsWrqV9Qiy//pdRrX5wxXCE051ljDEG+K/3tpGdc4j/+texdE9t53c4UcFaIsYYE4bs3V/xu3e3cdmYnnxjVIbf4USNcAZgnAjcB5zi6gugqtp2n2kzxrQppRVV/GDhWrp1SOT+mSP8DieqhNMSeQZ4FDgbOAPIct8nJCLtRGSliGSLyEYRud+V9xeRFSLyuYi86Ca6QkQS3fo2t71fyLHuduVbReSihv+YxhjTeA8v3cKOgmIemTWaju3bzuCK4QgniRxS1aWquk9V9wc/YexXBpyvqqPxJrKaJiITgF8Bj6lqJnAQuNHVvxE4qKoDgcdcPURkGHA1MByYBvy3iLS9sQWMMb74++eFzP/nTq4/qx9nZ3b1O5yoE04SeU9EfiMiZ4rI2OCnvp3Uc8StxrtPcMytl135c8BlbnmmW8dtnyLeow8zgRdUtUxVvwC2AePC+eGMMeZkHDpawY9ezmZAt2TumjbE73CiUjhPZ41331khZWENwOhaDKuBgcB/AduBr1S10lXJwZt2F/e9G0BVK0XkEJDmyj8OOWzoPqHnmgPMAejb157dNsacvPsWbWRfURmv/MdZtE+wDpDahPOy4eTGHlxVq4AxItIJeBUYWls1913bA9d6gvKa55oHzAPIysr62nZjjGmIN9bn8eqaPdw6JZPRfWz4wLrUmURE5N9U9Y8icntt21X10XBPoqpfichyYALQSUTiXGukN5DrquUAfYAcEYkDOgIHQsqDQvcxxpgmt+9wKT9+dT2jenfke+cP9DucqHaieyLBAWFS6vickIh0cy0QRKQ9cAGwGXgPuNJVuw543S0vcuu47e+qqrryq93TW/2BTGBlWD+dMcY0kKpy51/WcbS8ikdnjyG+DQ+uGI46WyKq+j/u+/5GHjsDeM7dF4kBFqrqEhHZBLwgIr8A1uA9Qoz7/j8R2YbXArnanX+jiCwENgGVwM2um8wYY5rcn1fuYvnWAu69ZBgDu3fwO5yoJ94f+61LVlaWrlq1yu8wjDEtzM7CYqY/8SFjT+nE/31rPDExbWtsLBFZrapZ9dc8xtppxhgDVAWUO17KJi5W+M2Vo9tcAmksG4DRGGOAp9/fzuovD/L4VWPo2am93+G0GPW2REQkXUSeEZGlbn2YiNxY337GGNNSbMw9xON/+4wZI3swc0xPv8NpUcLpzpoPvAUEf7OfAbdFKiBjjGlOpRVV3P5iNp2SEnjwspE2R0gDhZNEuqrqQiAA3tvkgD0dZYxpFR5d9hlb84v49ZWj6Jyc4Hc4LU44SaRYRNJwb4m7QRQPRTQqY4xpBh/v2M/vP9zBv47vy+TB3f0Op0UK58b6HXgv/J0qIv8AugGzIhqVMcZEWFFpBXcszKZvlyR+MqO2EZlMOMIZO2u1iJwHDMYbx2qrqlZEPDJjjImgB5ZsIu/QUV767pkkJ5qrH+kAABtLSURBVNqDqo0VztNZ24GbVHWjqm5Q1QoRWdIMsRljTES8vXEvC1fl8N3zTuX0U7r4HU6LFs49kQpgsoj8ITgLIbUMxW6MMS1B4ZEy7n5lPUMzUrntgkF+h9PihZNESlT1KrzBEz8UkVOoZSh2Y4yJdqrK3a+sp6i0ksevGkNCnA3acbLC6QgUAFX9tYisxntnxNp/xpgW5+XVOSzblM+PZwxhcI96ByM3YQgnifwsuKCq74jIRRwbst0YY1qE3QdKuH/xJsb178KNZw/wO5xW40STUg1R1S3AnlrmVLcb68aYFiMQUH74Ujaqym9njSbWBldsMidqidyON2f5b2vZFtYc68YYEw2e+fsXrPjiAL/+l1H06ZLkdzityokmpZrjvhs9x7oxxvht694ifvPWVi4Yms6srN5+h9PqhPOeyCwRSXHL94jIKyJyWuRDM8aYk5N/uJRbX1hDSrs4Hv4XG1wxEsJ5vu2nqlokImcDFwHPAU/Xt5OI9BGR90Rks4hsFJFbXXkXEVkmIp+7786uXETkSRHZJiLrQu/DiMh1rv7nImI39Y0x9Xpncz7THv+AL/eX8Mjs0XTtkOh3SK1SOEkkOGLvN4CnVPV1IJyhLiuBO1R1KDABuFlEhgFzgXdUNRN4x60DTAcy3WcO8BR4SQe4FxgPjAPuDSYeY4ypqbSiivsWbeTG51aR0bE9i79/tg2uGEHhPOK7R0T+B7gA+JWIJBJG8lHVPCDPLReJyGa8N91nApNcteeA5cBdrvx59SZ9/1hEOolIhqu7TFUPAIjIMmAasCDMn9EY00Z8nl/E9xesYcveIr41sT93TR9MYlys32G1auEkkdl4F+1HVPUrd2H/UUNOIiL9gNOAFUC6SzCoap6IBP9E6AXsDtktx5XVVV7zHHPwWjD07du3IeEZY1o4VWXByt38fMlGkhPi+MP1ZzB5iLU+mkM4o/iWAK+ErFe3MMIhIh2AvwC3qerhE9zYqm2DnqC8ZpzzgHkAWVlZNiyLMW3EVyXlzP3Let7cuJezB3bl0dmj6Z7azu+w2oyIjn8sIvF4CeRPqhpMRPkikuFaIRnAPleeA/QJ2b03kOvKJ9UoXx7JuI0xLcPKLw5w2wtr2FdUxt3Th/DtcwYQYy8SNquIjT4mXpPjGWCzqj4asmkRx4ZNuQ54PaT8WveU1gTgkGv1vAVMFZHO7ob6VFdmjGmjKqsCPLbsM66e9xHxcTH85T/O4jvnnWoJxAeRbIlMBP4dWC8ia13Zj4GHgYUiciOwi2OzJL4BzAC2ASXADQCqekBEHgA+cfV+HrzJboxpe3IOlnDbC2tZ9eVBrhjbi5/PHEEHm1TKN+I9DNW6ZGVl6apVq/wOwxjTxP66Lo+5r6xDFX5x2QguO82mNmpKIrJaVbMaso+lb2NM1Cspr+SBJZtYsHI3o/t04smrx3BKWrLfYRksiRhjotym3MN8f8Gn7Cgs5j8mncrtFw4iPtYmk4oWlkSMMVFJVZn/z5089MYWOiXF88cbxzNxYFe/wzI1WBIxxkSd/UfK+NHL63h3yz6mDOnOr68cRZqNfRWVLIkYY6LK3z8v5AcL13LoaAX3Xzqca888xUbfjWKWRIwxUaG8MsBvl21l3gc7OLVbB57/1jiGZqT6HZaphyURY4zvvtxfzC0L1pCdc4hvjuvLzy4eRvsEGzixJbAkYozx1Suf5vDT1zYQGyM8dc1Ypo/M8Dsk0wCWRIwxvigqreBnr2/k1TV7GNevC49dPYZendr7HZZpIEsixphml737K255YQ27D5TwgwsGcfPkU4mzdz9aJEsixphmEwgo//PBDn779lbSU9vx4nfO5Ix+XfwOy5wESyLGmGax73Apty/M5u/bCpkxsgcPXT6KjknxfodlTpIlEWNMxL2zOZ8fvbyOkvJKHr5iJFed0cfe/WglLIkYYyKmtKKKh5duYf4/dzI0I5XffXMMA7un+B2WaUKWRIwxEbFtXxHfX7CWzXmHuWFiP+6aNoR28fbuR2tjScQY06RUlRc+2c39izeSlBDHs9dncf6QdL/DMhFiScQY02QOlVRw96vreGP9Xs4e2JVHZ4+me2o7v8My9QkEYPfHjdo1YklERJ4FLgb2qeoIV9YFeBHoB+wEZqvqQTcf+xN40+OWANer6qdun+uAe9xhf6Gqz0UqZmNM4+QdOspra3J5/qOdFBSVMXf6EOacM8DmPI92BZ/Buhdh3UI4tKtRh4hkS2Q+8J/A8yFlc4F3VPVhEZnr1u8CpgOZ7jMeeAoY75LOvUAWoMBqEVmkqgcjGLcxJgxFpRUs3bCX19bs4aMd+1GF00/pzNP/djqj+3TyOzxTlyMFsOEvsO4FyF0DEgMDJsOUn8L9VzX4cBFLIqr6gYj0q1E8E5jklp8DluMlkZnA8+pN+P6xiHQSkQxXd5mqHgAQkWXANGBBpOI2xtStoirAh58X8OqaXN7euJeyygCnpCVx65RMLj+tl01ZG63KS2DrG16rY9s7oFXQYxRc9EsYcSWkBO9ZRVESqUO6quYBqGqeiHR35b2A3SH1clxZXeVfIyJzgDkAffv2beKwjWm7VJX1ew7xyqd7WJydy/7icjolxTM7qw+Xj+3FaX062Tsf0SgQgJ0feolj0yIoL4LU3jDxFhh1FXQf2iSniZYb67X9F6gnKP96oeo8YB5AVlZWrXWMMeHbfaCE19fu4dU1e9heUExCXAwXDO3O5af15rxB3UiIs7GuolL+Jq+rav3LcHgPJKTA8Jle4jjlbIhp2n+35k4i+SKS4VohGcA+V54D9Amp1xvIdeWTapQvb4Y4jWmTDh2t4I31eby6Zg8rvzgAwLj+XbjpnAHMGJlBx/Y2TElUKtoL61/yWh1710NMHAy8AKY+AINnQHzkRkdu7iSyCLgOeNh9vx5S/j0ReQHvxvohl2jeAn4pIp1dvanA3c0cszGtWnllgOVb9/Ha2j38bfM+yisDDOiWzA+nDmLmmF706ZLkd4imNmVHYMtfvVbHjuWgAeg5Fqb/Gkb8CyR3bZYwIvmI7wK8VkRXEcnBe8rqYWChiNwI7AJmuepv4D3euw3vEd8bAFT1gIg8AHzi6v08eJPdGNN4qsqnu77itTV7WLIul4MlFaQlJ/Cv4/pyxdhejOzV0e5zRKNAlZcw1r0Im5dARTF06gvn3OF1V3XNbPaQxHsgqnXJysrSVatW+R2GMVHny/3FvLpmD6+t2cPO/SUkxsUwdXgPrjitF2dndiXe5vSIPqpeF9W6F70uqyP50K4jDL8cRl0NfcY32X0OEVmtqlkN2SdabqwbYyLkYHE5S9bn8eqnOXy66ytE4MwBadw8eSDTRvQgpZ3d54hKh/bA+oXei4D7NkFMPAy6yGtxZE6F+OgYCcCSiDGtUFllFe9u3scra/awfOs+KqqUwekpzJ0+hJljepLR0aahjUqlh2HzYu8+xxcfAuq1NL7xqNfySIq+CbwsiRjTSgQCyqovD/Lqmhz+ui6Pw6WVdE9J5Pqz+nHZab0YlpFq9zmiUVUFbH/PSxxb3oDKo9BlAEyaC6Nme8tRzJKIMS3c9oIjvPrpHl5bu4ecg0dJSohl2vAeXHZaLyYO7EqsjV8VfVS9IUfWvei9z1FSCO27wGnXePc5emdBC0n4lkSMaWHKKwNsLzjCih37eXXNHrJzDhEjMHFgV+6YOoipw3qQnGj/a0eNqko4sB3yN0D+Rti7wbtRXpQLsYkweJqXOAZeAHEJfkfbYPZfmjFRbF9RKVvyitiy9zCb84rYnHeY7QVHqKjynqoclpHKPd8YyqWje9qQ69Gg5MDxySJ/AxRsgcpSb3tMHHQdBP0mQr9zYNhMaN+yB6u0JGJMFCirrGLbviNscYliy14vcRQeKa+uk9GxHUN6pHD+kO4MyUhlZK+O9O9qAx76oqoS9n/uJYv8DS5hbPRaF0FJXaHHCDjjJkgf4S13HQRxif7FHQGWRIxpRqrKvqIyNud5LYstew+zJa+I7QVHqAx4rYvEuBgGu2QxNCOVIT1SGZqRQqekltfV0SoUF9bSutgKVWXe9ph46DYY+p8L6cO9T4+R0KH7iY/bSlgSMSZCSiu81sWmvMMhXVKHOVhSUV2nV6f2DOmRwoXD0hmSkcKQHqn075psN8P9UFUBhZ+5ZLHetTI2wpG9x+p0SPdaFQMmHWtdpGW2yHsZTcWSiDEnSVXJO1Rafd9iy16vS+qLwmKqXOuiXXwMg3ukctHwHq514SWMjkn2op8vjuyr0brY6N27CLgEH5sA3YbAqee7lsUI6D4cOnTzN+4oZEnEmAY4Wl7FZ/nH3+jesreIQ0ePtS56d27PkB6pzBjRgyEuYZySZq2LZldVAcUF3jAhBVuPb10U7ztWLyXDa1UMnOJ1Q6UPh7SBEGsJPhyWRIxxAgHlq6MVFB4po7CojIIjZRQUlVF4pJzdB0rYvPcwOwuLcY0LkhJiGdwjhRkjMxiWkcKQjFQG90gh1YYRiZzKci8xFO/zpnk9kn9suXif18IIJo6jNWbRjk30JmLKnHp86yI5zZ+fpZWwJGJatUBAOVhSTuGRcgqrk4KXIAqLvLJg+YHi8uqb26HiY4WMjt69i4tH9fQSRo9U+nZJIsZaFyevssxd+EMSQPVyyPeRfCj9qvZjJKR4XU3J3d0jtGd7y9VlmdDlVIi1S15Ts9+oaXGqqhODlwgKjpRWJ4SCI17LIdiSOFBcXn1fIlR8rNC1QyLdUhJJT23H8J6pdO2Q6H1SEunWIZFuKQl07ZBIx/bxNlxIQ1WWuQRQV4shpKz0UO3HSEyF5G7eU07dh3hPP3Xofqwsubv33aF7RCddMidmScT4qrIqQHF5FUfLqygur+RoeRUHS8qrWwyhCaHwSLlrMZRRS14gIS6Gbh0S6dohgYyO7RjZqyNdXSLolpJYnSS6dUgktX2cJYYgVe+iX1EC5UegvBjK3XJFiVt3n4qQ5drqlBV5j8SW1ZUYOh5rHaQPgw6Tj28xBJNCcjdLDC2EJRETloqqACVlVZRUVFJcduyiX1JeSUl5lbetvPK4hODVr6KkzNWpub28ivLKwAnPmxgXU9066NWpHaN7dzwuIXTtkEBXt57arhUmhkAAqsrdp8J7eii4XFlW48Je86J+xJXVUqe8+PiLv1aFH1NMPCQkQ0IHSEjyluOTvcdfuwxwLYWQpJAcmhjsrfrWxpJIlAsElIpAgIoqpaIyQEVVgIrAseXyKm9bZchyRWWAykCA8tB9gtuqgvt5y0fdxd27yB9bLi6rdBd776JfXnXii32oGIHkhDiSEmNJToijfYL33SkpgV6dY0lKiCMpIfQ7luREb7l9fCydkhLo2iGBbimJdEisJTGoelOBasCb6U2rIHAUSqu8i65WufIay9V13bcGatQP/Q4cv3+g8usX8dCLe1Ut5YHKGnVq1q0lMdSs05CLe02xie4i3wHi3cU+IRlSe7oLf9LxiaBmveAnPric5C234XcizNe1yiRSfjCHXS/+CFUloN5z/Ip638F1b4EAx9ZVqbWeN/mjElDFW/SOiwZQOLashGz3jq/VxwfVgHeMgLc9uBwsD1SfL+DqeOcVgn03ikD1utSxHqwbixIHtK+jboxAbAzEixAXo8TFQJwrixPvHmRcvBLb4VhZjECsKLHuO0YgFiXGLce4eEUVCF7sFSoCUKFwxK1Xb/N+d8fV1dCLfh0X/2gUm+i9XxAb574TvMdEYxO8v96Dy7HxEN/RrccfXzcm/vj9TnS8uIQTX/jtJrJpBi3mvzIRmQY8AcQC/6uqD9dVN/5oIemb/gBQffn1luVryw3ZTo3tIlRfkmtSBAHU/RUtwWNW/1XtLUvIsvcJWXfn8L4l+IvwlkW8MwTXXd1jyyH1XN3jy0CI8c4nMe6cMS6GGuvHbau5HhMSf83jNKaugMRCTKzbFrIcExuyHiyLOb7suLon2lbzmHXUDdaJiavl4h5y0Y+JbTFDdxvTlFpEEhGRWOC/gAuBHOATEVmkqptqq1/SdQQbvvUesTFCrAgxMYQse9+xMd6y95f1sfKYGG9bXHA/V7fV9bUbY0wTaBFJBBgHbFPVHQAi8gIwE6g1iSQnxHH6KZ2bMTxjjGmbWkoS6QXsDlnPAcaHVhCROcAct1omIhuaKbaG6AoU+h1EDRZTeCym8EVjXBZTeAY3dIeWkkRq60s67naFqs4D5gGIyCpVzWqOwBoiGuOymMJjMYUvGuOymMIjIqsauk9MJAKJgBygT8h6byC3jrrGGGOaSUtJIp8AmSLSX0QSgKuBRT7HZIwxbV6L6M5S1UoR+R7wFt4jvs+q6sYT7DKveSJrsGiMy2IKj8UUvmiMy2IKT4NjEtVaBiEyxhhjwtBSurOMMcZEIUsixhhjGq1VJREReVZE9kXTOyIi0kdE3hORzSKyUURujYKY2onIShHJdjHd73dMQSISKyJrRGSJ37EEichOEVkvImsb8whkJIhIJxF5WUS2uP+2zvQ5nsHu9xP8HBaR2/yMycX1A/ff+AYRWSAivg8jLCK3ung2+vk7qu16KSJdRGSZiHzuvut9a7tVJRFgPjDN7yBqqATuUNWhwATgZhEZ5nNMZcD5qjoaGANME5EJPscUdCuw2e8gajFZVcdE0XP9TwBvquoQYDQ+/85Udav7/YwBTgdKgFf9jElEegG3AFmqOgLvoZyrfY5pBPBtvFE4RgMXi0imT+HM5+vXy7nAO6qaCbzj1k+oVSURVf0AOOB3HKFUNU9VP3XLRXj/s/fyOSZV1SNuNd59fH/CQkR6A98A/tfvWKKZiKQC5wLPAKhquarWMW+sL6YA21X1S78DwXsCtb2IxAFJ+P9+2VDgY1UtUdVK4H3gcj8CqeN6ORN4zi0/B1xW33FaVRKJdiLSDzgNWOFvJNXdRmuBfcAyVfU9JuBx4E4g2sZ6V+BtEVnthtfx2wCgAPiD6/r7XxFJ9juoEFcDC/wOQlX3AI8Au4A84JCqvu1vVGwAzhWRNBFJAmZw/IvUfktX1Tzw/gAGute3gyWRZiIiHYC/ALep6mG/41HVKtf10BsY55rZvhGRi4F9qrrazzjqMFFVxwLT8bojz/U5njhgLPCUqp4GFBNGt0NzcC8DXwq8FAWxdMb7y7o/0BNIFpF/8zMmVd0M/ApYBrwJZON1ebdYlkSagYjE4yWQP6nqK37HE8p1gyzH/3tJE4FLRWQn8AJwvoj80d+QPKqa67734fXzj/M3InKAnJDW48t4SSUaTAc+VdV8vwMBLgC+UNUCVa0AXgHO8jkmVPUZVR2rqufidSd97ndMIfJFJAPAfe+rbwdLIhEm3kQkzwCbVfVRv+MBEJFuItLJLbfH+59ti58xqerdqtpbVfvhdYe8q6q+/tUIICLJIpISXAam4nVJ+EZV9wK7RSQ44uoU6pgWwQffJAq6spxdwAQRSXL/H04hCh7aEJHu7rsvcAXR8/sCbzip69zydcDr9e3QIoY9CZeILAAmAV1FJAe4V1Wf8TcqJgL/Dqx39yAAfqyqb/gYUwbwnJvsKwZYqKpR80htlEkHXnWTksUBf1bVN/0NCYDvA39y3Uc7gBt8jgfXx38h8B2/YwFQ1RUi8jLwKV6X0RqiY6iRv4hIGlAB3KyqB/0IorbrJfAwsFBEbsRLwrPqPY4Ne2KMMaaxrDvLGGNMo1kSMcYY02iWRIwxxjSaJRFjjDGNZknEGGNMo1kSMVFJRG5zj4z6ce7fuBFWf9PEx80SkSfrqXO9iPxnHduO1FYeKSKyXEROOOCkn/9OJjq0qvdETKtyG/BHvNFgm9t3gG6qWtZUBxSROFVdBUTFUPJNyM9/JxMFrCVifOXeCP+rm9tkg4hcJSK34I119J6IvOfqTRWRj0TkUxF5yY1FFpzr41dufpSVIjLQlc9yx8sWkQ9qOa+4FscG8eYKucqVLwKSgRXBMlce487VKaRsm4iki8glIrLCDYb4NxFJd9vvE5F5IvI28LyITBI3T4qIjBORf7p9/hny9jlAHxF5U0S2isi9dfzefiQin4jIOqljPpjQlouIXCki893yfBF5WkQ+FJHP3LhliEh7EXnBHfNFoH3I/k+JyCoJmX+mIf9OphVTVfvYx7cP8C/A70PWO7rvnUBXt9wV+ABIdut3AT8LqfcTt3wtsMQtrwd6ueVOdZx3Gd4cE+l4b+dmuG1H6oj1CeAGtzwe+Jtb7syxF3dvAn7rlu8DVgPt3fqkkPhSgTi3fAHwF7d8Pd6Is2l4F/ENePNhVMeFN/TKPEDw/hBcApxbS7xHQpavBOa75fl4g//FAJl4Y3G1A24HnnV1RuG95R08dxf3HYs31tqohvw72af1fqwlYvy2HrjAtSbOUdVDtdSZAAwD/uGGjrkOOCVk+4KQ7+AMf/8A5ovIt/EufDWdDSxQbzTjfLx5Hc6oJ9YXgWDr5Gq3Dt5IyG+JyHrgR8DwkH0WqerRWo7VEXhJvFnlHquxzzJV3e/2e8XFGmqq+6zBG9JjCF4yaIiFqhpQ1c/xhk0ZgjdHyR8BVHUdsC6k/mwR+dSdczjev0dN9f07mVbI7okYX6nqZyJyOt68Cg+JyNuq+vMa1QTvwvrNug5Tc1lVvysi4/EmuVorImNUdX+NYzbUR8BAEemGN1nPL1z574BHVXWRiEzCa4EEFddxrAeA91T1cvHmmVlex89T27oAD6nq/9QTb+h+NaeFrescXxsHSUT6Az8EzlDVg65brLZpZuv7dzKtkLVEjK9EpCdQoqp/xJtAKDikeRGQ4pY/BiaG3O9IEpFBIYe5KuT7I1fnVFVdoao/Awr5+sQ/HwBXiTc5Vze8v8JXnihWVVW8oeAfxRuVOZiUOgJ73PJ1te1bi9B9rq+x7ULx5rpuj5es/lFj+1vAt0LuC/USNzJsDfkiMlREYvj67Hmz3H2eU/EmudqK9zu5xh1zBF6XFnhdb8XAIXe/Z3rIcRry72RaIWuJGL+NBH4jIgG8UU3/w5XPA5aKSJ6qThaR64EFIpLott8DfOaWE0VkBd4fRcG/gn8j3tzVgjdXdHaN876K1/WVjffX953qDbFenxeBTzj+wn8fXtfUHrwLaf8wjvNrvJGUbwferbHt78D/AQPxRg0+7okuVX1bRIYCH4k3uvAR4N/4+twPc/Hul+zGu7cSepN7K14XXjrwXVUtFZGn8GZLXAesxSVVVc0WkTXARryur9Ck1pB/J9MK2Si+pkUTbxKrLFUt9DuWlsJ1Ry1R1Zf9jsW0fNadZYwxptGsJWKMMabRrCVijDGm0SyJGGOMaTRLIsYYYxrNkogxxphGsyRijDGm0f4/IeIZXAWC+tMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = [0,] # Python list\n", "grow = lambda x: [0,]*len(x)*2\n", "d1 = memory_usage(x,grow,steps=10)\n", "x = np.array([0])\n", "grow = lambda x: np.array([0,]*x.size*2,dtype='int8')\n", "d2 = memory_usage(x,grow,steps=10)\n", "\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.plot(d1[\"x\"],d1[\"y\"],label='Python')\n", "plt.plot(d2[\"x\"],d2[\"y\"],label='Numpy')\n", "plt.axis([min(d1[\"x\"]),max(d1[\"x\"]),0,max(d1[\"y\"])+1])\n", "plt.ylabel('size in memory, bytes')\n", "plt.xlabel('steps of variable update')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Creating arrays\n", "\n", "- From lists \n", "- Using functions for standard cases " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 3. 5. 7.]\n" ] } ], "source": [ "a = np.array([1,3,5.0,7])\n", "print(a)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0 0 0 0 0 0 0 -96 0 0 0 0 0 0 0 -96 41 0\n", " -18 10 1 0 0 0 9]\n", "\n", "[0. 0. 0. 0. 0.]\n", "\n", "[1. 1. 1. 1. 1.]\n", "\n", "[0 1 2 3 4 5 6 7 8 9]\n", "\n", "[2. 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. ]\n" ] } ], "source": [ "a = np.empty(25,'int8') # not initialized !\n", "b = np.zeros(5) # initialized with zeros\n", "c = np.ones(5)\n", "d = np.arange(10)\n", "e = np.linspace(2, 3, 11) # fill between 2 and 3 with 10 points\n", "print(a,b,c,d,e,sep='\\n\\n')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Note that uninitialized array may have garbage (random state of the memory)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Matrices and multidimensional arrays" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1. 0. 0. 0. 0.]\n", " [0. 1. 0. 0. 0.]\n", " [0. 0. 1. 0. 0.]\n", " [0. 0. 0. 1. 0.]\n", " [0. 0. 0. 0. 1.]]\n", "[[1. 1. 1.]\n", " [1. 1. 1.]]\n" ] } ], "source": [ "a = np.eye(5) # identity matrix\n", "b = np.ones((2,3))\n", "print(a)\n", "print(b)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3. 3. 3.]\n", " [3. 3. 3.]]\n", "\n", "[[3. 3. 3.]\n", " [3. 3. 3.]]\n", "\n" ] } ], "source": [ "b=b+2\n", "c = np.asmatrix(b) # matrix !\n", "print(b)\n", "print(type(b))\n", "print(c)\n", "print(type(c))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3. 3. 3.]\n", " [3. 3. 3.]]\n", "[[9. 9. 9.]\n", " [9. 9. 9.]]\n" ] }, { "ename": "ValueError", "evalue": "shapes (2,3) and (2,3) not aligned: 3 (dim 1) != 2 (dim 0)", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# element by element\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# matrix multiplication\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/usr/local/anaconda3/lib/python3.7/site-packages/numpy/matrixlib/defmatrix.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 218\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 219\u001b[0m \u001b[0;31m# This promotes 1-D vectors to row vectors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 220\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0masmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 221\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misscalar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'__rmul__'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 222\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mdot\u001b[0;34m(*args, **kwargs)\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: shapes (2,3) and (2,3) not aligned: 3 (dim 1) != 2 (dim 0)" ] } ], "source": [ "print(b)\n", "print(b*b) # element by element\n", "print(c*c) # matrix multiplication" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[27. 27.]\n", " [27. 27.]]\n" ] } ], "source": [ "# c*c\n", "print(c*c.transpose())" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Indexing into arrays\n", "\n", "Several types of indexes:\n", "\n", "- scalar index x[0] (getitem) \n", "- slicing like strings x[1::2] \n", "- numerical indexing \n", "- masks " ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0.14285714 0.28571429]\n", " [0.42857143 0.57142857 0.71428571]\n", " [0.85714286 1. 1.14285714]\n", " [1.28571429 1.42857143 1.57142857]\n", " [1.71428571 1.85714286 2. ]]\n" ] } ], "source": [ "z = np.linspace(0, 2, 15)\n", "z = np.reshape(z,[5,3])\n", "print(z)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0.14285714 0.28571429]\n", " [0.42857143 0.57142857 0.71428571]\n", " [0.85714286 1. 1.14285714]\n", " [1.28571429 1.42857143 1.57142857]\n", " [1.71428571 1.85714286 2. ]]\n", "\n", "[0.42857143 0.57142857 0.71428571]\n", "0.42857142857142855\n", "[[1.71428571 1.85714286 2. ]]\n", "[0.57142857 1. ]\n", "[[0.57142857 0.71428571]\n", " [1. 1.14285714]]\n", "[0.14285714 0.57142857 1. 1.42857143 1.85714286]\n" ] } ], "source": [ "print(z,end='\\n\\n')\n", "print( z[1] ) # scalar index: returns row array\n", "print( z[1,0] ) # scalar index: returns number\n", "print( z[-1:] ) # slicing: returns ?\n", "print( z[1:3,1] ) # slicing + scalar index\n", "print( z[1:3,1:] ) # slicing\n", "print( z[:,1] ) # slicing to get the column" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 4.2 9.3 0.28571429]\n", " [ 5.2 9.3 0.71428571]\n", " [ 6.2 5. 7. ]\n", " [ 7.2 -2. 1.57142857]\n", " [ 8.2 1.85714286 2. ]]\n" ] } ], "source": [ "# Assigning elements of an array\n", "z[1,0] = -1\n", "z[2] = [4,5,7] # assign whole row from a list\n", "z[:,0] = np.array([4.2,5.2,6.2,7.2,8.2]) # assign column from nparray\n", "z[:2,1]=9.3\n", "z[3][1]=-2 # note double bracket indexing\n", "print(z)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0.18181818 0.36363636]\n", " [0.54545455 0.72727273 0.90909091]\n", " [1.09090909 1.27272727 1.45454545]\n", " [1.63636364 1.81818182 2. ]]\n", "\n", "[0. 1.27272727 1.09090909]\n", "[1.09090909 1.27272727 1.45454545 1.63636364 1.81818182 2. ]\n", "[[False False False]\n", " [False False False]\n", " [ True True True]\n", " [ True False False]]\n", "\n", "[1.09090909 1.27272727 1.45454545 1.63636364]\n" ] } ], "source": [ "z = np.linspace(0, 2, 12)\n", "z = np.reshape(z,[4,3])\n", "print(z,end='\\n\\n')\n", "\n", "print( z[[0,2,2],[0,1,0]] ) # numerical (element by element) indexing\n", "print( z[z>1.0] ) # boolean indexing (masking)\n", "mask = np.logical_and(z>1.0,z<1.75)\n", "print(mask,end='\\n\\n')\n", "print( z[mask] ) # boolean indexing (masking)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Operation broadcasting\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Operation broadcasting" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[11 12 13]\n", "[5 7 9]\n" ] } ], "source": [ "a = np.array([1,2,3])\n", "b = np.array([4,5,6])\n", "print(a + 10)\n", "print(a + b)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5 6 7]\n", "\n", "[[1. 1. 1.]\n", " [1. 1. 1.]\n", " [1. 1. 1.]]\n", "\n", "[[6. 7. 8.]\n", " [6. 7. 8.]\n", " [6. 7. 8.]]\n" ] } ], "source": [ "x = np.arange(3) + 5\n", "y = np.ones((3,3))\n", "print(x,y,x+y,sep='\\n\\n')" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0]\n", " [1]\n", " [2]\n", " [3]\n", " [4]]\n", "\n", "[[0 1 2 3 4]]\n", "\n", "[[0 1 2 3 4]\n", " [1 2 3 4 5]\n", " [2 3 4 5 6]\n", " [3 4 5 6 7]\n", " [4 5 6 7 8]]\n" ] } ], "source": [ "x = np.arange(5).reshape((5,1)) # or\n", "x = np.arange(5)[:,np.newaxis]\n", "print(x,end='\\n\\n')\n", "print(x.transpose(),end='\\n\\n')\n", "print(x + x.transpose())" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0 1 2 3]\n", " [ 4 5 6 7]\n", " [ 8 9 10 11]]\n", "\n", "[0 1 2 3]\n", "\n", "[[ 0 2 4 6]\n", " [ 4 6 8 10]\n", " [ 8 10 12 14]]\n" ] } ], "source": [ "x = np.arange(12).reshape((3,4))\n", "y = np.arange(4)\n", "print(x,y,x+y,sep='\\n\\n')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### Broadcasting works with:\n", "\n", "- addition ($ + $) \n", "- subtraction ($ - $) \n", "- multiplication ($ * $) \n", "- division ($ / $) \n", "- integer division ($ // $) \n", "- power ($ ** $) \n", "- all *universal functions* " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### ufunc\n", "\n", "Functions provided by NumPy which support vectorization and broadcasting\n", "\n", "- Act on array element-wise \n", "- Efficient implementation with low level code \n", "\n", "\n", "[https://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs](https://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.02 ms ± 15.4 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)\n", "101 µs ± 11.6 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)\n", "NumPy is on avarage 10.158 faster\n" ] } ], "source": [ "import math\n", "N = 10000\n", "data_list = list(range(N)) # Native Python list\n", "t1 = %timeit -n10 -r10 -o sin1 = [math.sin(x) for x in data_list]\n", "\n", "import numpy as np\n", "data_array = np.array(range(N)) #NumPy array\n", "t2 = %timeit -n10 -r10 -o sin2 = np.sin(data_array)\n", "\n", "print('NumPy is on avarage %2.3f faster' % (t1.average/t2.average))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Reduction operations\n", "\n", "Functions that return the array of reduced size: **sum**, **min**,\n", "**max** , **mean**, **all**, **any**\n", "\n", "[https://numpy.org/doc/stable/reference/routines.math.html](https://numpy.org/doc/stable/reference/routines.math.html)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0 1 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n", "66\n", "[ 3 12 21 30]\n", "[[ 2]\n", " [ 5]\n", " [ 8]\n", " [11]]\n" ] } ], "source": [ "x = np.arange(12).reshape(4,3)\n", "print(x)\n", "print(np.sum(x))\n", "print(np.sum(x, axis=1))\n", "print(np.maximum.reduce(x,axis=1,keepdims=True))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[[ 0 1 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n", "\n", " [[12 13 14]\n", " [15 16 17]\n", " [18 19 20]\n", " [21 22 23]]]\n", "(2, 4)\n", "[[ 0 3 6 9]\n", " [12 15 18 21]]\n" ] } ], "source": [ "x = np.arange(24).reshape(2,4,3)\n", "print(x)\n", "y = np.min(x,axis=2)\n", "# y = np.mean(x,axis=(1,2))\n", "print(y.shape)\n", "print(y)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### References and mutability\n", "\n", "NumPy tries not to copy data in the arrays when not necessary\n", "\n", "- principle: whether it is possible to maintain simple pointer arithmetic \n", "- slices are generally not copied \n", "- numerical indexing/mask generally copied \n", "- **.flags** to check \n", "- **.copy** to make a true copy " ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0 1 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n", "[[ 1]\n", " [ 4]\n", " [ 7]\n", " [10]]\n", " C_CONTIGUOUS : False\n", " F_CONTIGUOUS : False\n", " OWNDATA : False\n", " WRITEABLE : True\n", " ALIGNED : True\n", " WRITEBACKIFCOPY : False\n", " UPDATEIFCOPY : False\n", "\n", "[[ 0 999 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n" ] } ], "source": [ "x = np.arange(12).reshape(4,3) # 2-dim array\n", "print(x)\n", "y = x[:,1:2]\n", "print(y)\n", "print(y.flags)\n", "y[0] = 999\n", "print(x)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 5]\n", " C_CONTIGUOUS : True\n", " F_CONTIGUOUS : True\n", " OWNDATA : True\n", " WRITEABLE : True\n", " ALIGNED : True\n", " WRITEBACKIFCOPY : False\n", " UPDATEIFCOPY : False\n", "\n", "[[ 0 999 2]\n", " [ 3 4 5]\n", " [ 6 7 8]\n", " [ 9 10 11]]\n" ] } ], "source": [ "y = x[[0,1],[0,2]]\n", "print(y)\n", "print(y.flags)\n", "y[0]=-100\n", "print(x)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "[https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flags.html](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flags.html)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Linear algebra with NumPy\n", "\n", "Submodule numpy.linalg\n", "\n", "- matrix decompositions \n", "- eigenvalues \n", "- determinant, rank \n", "- matrix inverse \n", "- linear systems of equations \n", "\n", "\n", "[https://numpy.org/doc/stable/reference/routines.linalg.html](https://numpy.org/doc/stable/reference/routines.linalg.html)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "------------------------------------------------------------\n", "Object report on object = \n", "Objec class : \n", "Parent classes : \n", "Occupied memory : 88 bytes\n", "PUBLIC METHODS\n", "LinAlgError \n", "cholesky \n", "cond \n", "det \n", "eig \n", "eigh \n", "eigvals \n", "eigvalsh \n", "inv \n", "lstsq \n", "matrix_power \n", "matrix_rank \n", "multi_dot \n", "norm \n", "pinv \n", "qr \n", "slogdet \n", "solve \n", "svd \n", "tensorinv \n", "tensorsolve \n", "test \n" ] } ], "source": [ "import numpy.linalg as linalg\n", "obj_explore(linalg,'public methods')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Matrix operations\n", "\n", "$$\n", "A=\n", "\\begin{pmatrix}\n", "1 & 2 & 0 & 5 \\\\\n", "4 & -2 & 1 & 1 \\\\\n", "0 & 0 & -2 & 7 \\\\\n", "3 & 1 & 4 & 0 \\\\\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "hide-output": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1 2 0 5]\n", " [ 4 -2 1 1]\n", " [ 0 0 -2 7]\n", " [ 3 1 4 0]]\n" ] } ], "source": [ "A = np.array([[1,2,0,5],[4,-2,1,1],[0,0,-2,7],[3,1,4,0]])\n", "print(A)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1 2 0 5]\n", " [ 4 -2 1 1]\n", " [ 0 0 -2 7]\n", " [ 3 1 4 0]]\n", "\n" ] } ], "source": [ "print(A,end='\\n\\n')\n", "# print( A.transpose() )\n", "# print( np.linalg.matrix_rank(A) )\n", "# print( np.linalg.inv(A) )\n", "# print( np.linalg.det(A) )\n", "B = A[:2,:]\n", "# print(B,end='\\n\\n')\n", "# print( B * B ) # element by element\n", "# print( B @ B ) # matrix multiplication\n", "# print( B @ B.T )" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Linear systems of equations\n", "\n", "$$\n", "\\begin{eqnarray*}\n", "1x_1+2x_2+5x_4&=&5\\\\\n", "4x_1-2x_2+x_3+x_4&=&5\\\\\n", "-2x_3+7x_4&=&0\\\\\n", "3x_1+x_2+4x_3&=&-3\\\\\n", "\\end{eqnarray*}\n", "$$\n", "\n", "In matrix notation $ Ax=b $ where\n", "\n", "$$\n", "A=\n", "\\begin{pmatrix}\n", "1 & 2 & 0 & 5 \\\\\n", "4 & -2 & 1 & 1 \\\\\n", "0 & 0 & -2 & 7 \\\\\n", "3 & 1 & 4 & 0 \\\\\n", "\\end{pmatrix},\\;\\;\n", "b=\\begin{pmatrix}\n", "5\\\\\n", "5\\\\\n", "0\\\\\n", "-3\n", "\\end{pmatrix},\\;\\;\n", "x=\\begin{pmatrix}\n", "x_1\\\\\n", "x_2\\\\\n", "x_3\\\\\n", "x_4\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution is array([ 4.95555556, 3.91111111, -5.44444444, -1.55555556])\n", "Check: max(Ax-b) = 8.88178e-16\n" ] } ], "source": [ "b = np.array([5,5,0,-3])\n", "x=np.linalg.solve(A, b)\n", "print('Solution is %r' % x)\n", "print('Check: max(Ax-b) = %1.5e' % np.max(A@x-b))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Overdetermined linear systems of equations\n", "\n", "$$\n", "\\begin{eqnarray*}\n", "1x_1+2x_2&=&5\\\\\n", "4x_1-2x_2+x_3&=&5\\\\\n", "-2x_3&=&0\\\\\n", "3x_1+x_2+4x_3&=&-3\\\\\n", "\\end{eqnarray*}\n", "$$" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1 2 0]\n", " [ 4 -2 1]\n", " [ 0 0 -2]\n", " [ 3 1 4]]\n", "Solution is array([ 1.75294118, 0.63529412, -1.72941176])\n", "Check: max(Ax-b) = 3.45882e+00\n" ] } ], "source": [ "A = np.array([[1,2,0],[4,-2,1],[0,0,-2],[3,1,4]])\n", "x=np.linalg.lstsq(A, b, rcond=None)\n", "\n", "A = np.array([[1,2,0],[4,-2,1],[0,0,-2],[3,1,4]])\n", "x,*_=np.linalg.lstsq(A, b, rcond=None) # ignore all outputs except the first\n", "print(A)\n", "print('Solution is %r' % x)\n", "print('Check: max(Ax-b) = %1.5e' % np.max(A@x-b))\n", "# help(np.linalg.lstsq)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### Market equilibrium in a linear model\n", "\n", "- Prices $ p $, quantities $ q $ \n", "- $ n $ goods in the economy \n", "- Linear demand $ D(p) = A p + d $, where $ A $ is n by n, and\n", " $ p $ is n by 1 \n", "- Linear supply $ S(p) = B p + s $, where $ B $ is n by n, and\n", " $ p $ is n by 1 \n", "- Market clearing prices: $ D(p)=S(p) $ " ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "hide-output": false, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Demand is given by Ap+d:\n", "A=array([[ -1.50293121, -3.83492595, -2.53004055],\n", " [ -3.83492595, -14.91324372, -12.76421734],\n", " [ -2.53004055, -12.76421734, -12.93230177]])\n", "d=array([100, 100, 100])\n", "Supply is given by Bp+s:\n", "B=array([[4.96246147, 2.14375478, 4.0190379 ],\n", " [2.14375478, 1.24921004, 1.0231901 ],\n", " [4.0190379 , 1.0231901 , 5.54449637]])\n", "s=array([0., 0., 0.])\n", "Equilibrium prices are p = [15.18283269 1.4987406 -1.08770522]\n", "Equilibrium quantities are q = [74.185626 33.30758277 56.52309892]\n" ] } ], "source": [ "def random_matrix(n,positive=True,maxeigen=10):\n", " '''Generates square random positive/negative semi-definite matrix'''\n", " e=np.random.uniform(0,maxeigen,n) # random eigenvalues\n", " r=np.random.uniform(0,1,n*n).reshape(n,n) # rotation\n", " e = e if positive else -e\n", " A = np.diag(e) # diagonal matrix with\n", " return r @ A @ r.T # positive/negative semi-definite\n", "n = 3 # number of products\n", "A = random_matrix(n,positive=False) # demand\n", "d = np.array([100,]*n)\n", "B = random_matrix(n) # supply\n", "s = np.zeros(n)\n", "p = np.linalg.solve(A-B, s-d) # solve for quilibrium\n", "q = A @ p + d # equilibrium quantities\n", "print('Demand is given by Ap+d:\\nA=%r\\nd=%r' % (A,d))\n", "print('Supply is given by Bp+s:\\nB=%r\\ns=%r' % (B,s))\n", "print('Equilibrium prices are p = {}'.format(p))\n", "print('Equilibrium quantities are q = {}'.format(q))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Further learning resources\n", "\n", "- Reference manual for Numpy\n", " [https://numpy.org/doc/stable/reference/](https://numpy.org/doc/stable/reference/) \n", "- 📖 Kevin Sheppard “Introduction to Python for Econometrics, Statistics\n", " and Data Analysis.” *Chapters: 4, 6, 7, 8* \n", "- SciPy 2017 Tutorial on NumPy (2.5h)\n", " [https://www.youtube.com/watch?v=lKcwuPnSHIQ&ab_channel=Enthought](https://www.youtube.com/watch?v=lKcwuPnSHIQ&ab_channel=Enthought) \n", "- Essence of linear algebra playlist by 3Blue1Brown\n", " [https://www.youtube.com/watch?v=fNk_zzaMoSs&list=PLZHQObOWTQDPD3MizzM2xVFitgF8hE_ab](https://www.youtube.com/watch?v=fNk_zzaMoSs&list=PLZHQObOWTQDPD3MizzM2xVFitgF8hE_ab) " ] } ], "metadata": { "celltoolbar": "Slideshow", "date": 1612589584.7823849, "download_nb": false, "filename": "14_numpy.rst", "filename_with_path": "14_numpy", "kernelspec": { "display_name": "Python", "language": "python3", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" }, "title": "Foundations of Computational Economics #14" }, "nbformat": 4, "nbformat_minor": 4 }