14_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #14\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Vectors and matrixes (Numpy)\n",
    "\n",
    "<img src=\"_static/img/lecture.png\" style=\"width:64px;\">"
   ]
  },
  {
   "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",
    "<img src=\"_static/img/numpy_logo.png\" style=\"width:512px;\">\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": null,
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "obj_explore(a,'public')\n",
    "# help(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Memory footprint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array([1,3,5.0,7])\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a = np.eye(5) # identity matrix\n",
    "b = np.ones((2,3))\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(b)\n",
    "print(b*b) # element by element\n",
    "print(c*c) # matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "z = np.linspace(0, 2, 15)\n",
    "z = np.reshape(z,[5,3])\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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",
    "<img src=\"_static/img/broadcasting.png\" style=\"width:800px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Operation broadcasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": 1597824775.264069,
  "filename": "14_numpy.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Foundations of Computational Economics #14"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}