{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\re}{\\mathbb{R}}$\n",
    "Numpy provides many useful facilities to perform numerical computations including vectors, matrices and linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.141592653589793\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "print(numpy.pi)\n",
    "print(numpy.sin(numpy.pi/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " It is common practice to import numpy like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that ```np``` acts as an alias or short-hand for ```numpy```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.141592653589793\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "print(np.pi)\n",
    "print(np.sin(np.pi/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1-d Arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an array of zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "x = np.zeros(5)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These one dimensional arrays are of type ```ndarray```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an array of ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones(5)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create and add two arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5. 7. 9.]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1.0, 2.0, 3.0])\n",
    "y = np.array([4.0, 5.0, 6.0])\n",
    "z = x + y\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the size of array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "3\n"
     ]
    }
   ],
   "source": [
    "print(len(x))\n",
    "print(x.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the shape of an array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,)\n"
     ]
    }
   ],
   "source": [
    "print(x.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that these are arrays of reals by default. We can specify the type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "a = np.zeros(5, dtype=int)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## linspace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This behaves same way as Matlab's [linspace](https://mathworks.com/help/matlab/ref/linspace.html) function.\n",
    "\n",
    "Generate 10 uniformly spaced numbers in [1,10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(1,10,10)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this includes the end points 1 and 10. The output of linspace is an ```ndarray``` of floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`x = linspace(a,b,n)` is such that `x` is an array of `n` elements\n",
    "```\n",
    "x[i] = a + i*h,   i=0,1,2,...,n-1,    h = (b-a)/(n-1)\n",
    "```\n",
    "so that\n",
    "```\n",
    "x[0] = a,   x[-1] = b\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## arange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4 5 6 7 8 9]\n",
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(1,10)\n",
    "print(x)\n",
    "print(type(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $m < n$, `arange(m,n)` returns\n",
    "$$\n",
    "m, m+1, \\ldots, n-1\n",
    "$$\n",
    "Note the last value $n$ is not included. \n",
    "\n",
    "We can specify a step size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 5 7 9]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(1,10,2)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the arguments are float, the returned value is also float."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 2. 3. 4. 5. 6. 7. 8. 9.]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(1.0,10.0)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(0,1,0.1)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beware of pitfalls - 1\n",
    "\n",
    "Create an array of ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones(10)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we want set all elements to zero, so we might try this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "x = 0.0\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```x``` has changed from an array to a scalar !!! The correct way is this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones(10)\n",
    "print(x)\n",
    "x[:] = 0.0\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beware of pitfalls - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0.]\n",
      "[0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones(5)\n",
    "y = x\n",
    "x[:] = 0.0\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why did ```y``` change ? This happened because when we do\n",
    "```\n",
    "y = x\n",
    "```\n",
    "then `y` is just a pointer to `x`, so that changing `x` changes `y` also. If we want ```y``` to be an independent copy of ```x``` then do this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0.]\n",
      "[1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones(5)\n",
    "y = x.copy()     # or y = np.copy(x)\n",
    "x[:] = 0.0\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2-d Arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2-d arrays can be considered as matrices, though Numpy has a separate [matrix](https://numpy.org/doc/stable/reference/generated/numpy.matrix.html) class.\n",
    "\n",
    "Create an array of zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.zeros((5,5))\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the elements by two indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "A[1,1] = 1.0\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an array of ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 1.]\n",
      " [1. 1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.ones((2,3))\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create identity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "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"
     ]
    }
   ],
   "source": [
    "A = np.eye(5)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an array by specifying its elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2.]\n",
      " [3. 4.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1.0, 2.0], \n",
    "              [3.0, 4.0]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a random array and inspect its shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.56171896 0.66235847 0.78542121]\n",
      " [0.40789654 0.07574282 0.7198766 ]]\n",
      "A.size = 6\n",
      "A.shape = (2, 3)\n",
      "A.shape[0] = 2  (number of rows)\n",
      "A.shape[1] = 3  (number of columns)\n"
     ]
    }
   ],
   "source": [
    "m, n = 2, 3\n",
    "A = np.random.rand(m,n)\n",
    "print(A)\n",
    "print('A.size =',A.size)\n",
    "print('A.shape =',A.shape)\n",
    "print('A.shape[0] =',A.shape[0],' (number of rows)')\n",
    "print('A.shape[1] =',A.shape[1],' (number of columns)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To print the elements of an array, we need two for loops, one over rows and another over the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0 0.5617189572734113\n",
      "0 1 0.6623584720511232\n",
      "0 2 0.7854212087229677\n",
      "1 0 0.40789654277553433\n",
      "1 1 0.07574281632417834\n",
      "1 2 0.7198765978001941\n"
     ]
    }
   ],
   "source": [
    "for i in range(m): # loop over rows\n",
    "    for j in range(n): # loop over columns\n",
    "        print(i,j,A[i,j])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To transpose a 2-d array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "[[1 2]\n",
      " [3 4]]\n",
      "B =\n",
      "[[1 3]\n",
      " [2 4]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2],[3,4]])\n",
    "print(\"A =\"); print(A)\n",
    "B = A.T\n",
    "print(\"B =\"); print(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagonal matrix creation\n",
    "\n",
    "Create\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 2 & 0 & 0 \\\\\n",
    "0 & 0 & 3 & 0 \\\\\n",
    "0 & 0 & 0 & 4 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 0 0 0]\n",
      " [0 2 0 0]\n",
      " [0 0 3 0]\n",
      " [0 0 0 4]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1,2,3,4])    # diagonal elements\n",
    "A = np.diag(a)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create tri-diagonal matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4 -1  0  0]\n",
      " [ 1  5 -2  0]\n",
      " [ 0  2  6 -3]\n",
      " [ 0  0  3  7]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1,2,3])    # sub-diagonal  : length = n-1\n",
    "b = np.array([4,5,6,7])  # main diagonal : length = n\n",
    "c = np.array([-1,-2,-3]) # super-diagonal: length = n-1\n",
    "A = np.diag(a,-1) + np.diag(b,0) + np.diag(c,+1)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## empty_like, etc.\n",
    "\n",
    "Sometimes we have an existing array and we want to create a new array of the same shape and type, but without initializing its elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 3) (2, 3)\n"
     ]
    }
   ],
   "source": [
    "A = np.random.rand(2,3)\n",
    "B = np.empty_like(A)\n",
    "print(A.shape,B.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "C = np.zeros_like(A)\n",
    "print(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 1.]\n",
      " [1. 1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "D = np.ones_like(A)\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1-D array is neither row nor column vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,) [1 2 3]\n",
      "(3,) [1 2 3]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1,2,3])\n",
    "print(x.shape, x)\n",
    "y = x.T\n",
    "print(y.shape, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a row vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x.shape = (1, 3)\n",
      "[[1 2 3]]\n",
      "y.shape = (3, 1)\n",
      "[[1]\n",
      " [2]\n",
      " [3]]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([[1,2,3]]) # row vector\n",
    "print('x.shape =',x.shape)\n",
    "print(x)\n",
    "y = x.T\n",
    "print('y.shape =',y.shape)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a column vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x.shape = (3, 1)\n",
      "[[1]\n",
      " [2]\n",
      " [3]]\n",
      "y.shape = (1, 3)\n",
      "[[1 2 3]]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([[1],[2],[3]]) # column vector\n",
    "print('x.shape =',x.shape)\n",
    "print(x)\n",
    "y = x.T\n",
    "print('y.shape =',y.shape)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions like `zeros` and `ones` can return row/column vector rather than array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =\n",
      "[[1.]\n",
      " [1.]\n",
      " [1.]]\n",
      "y =\n",
      "[[1. 1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "x = np.ones((3,1)) # column vector\n",
    "print('x ='); print(x)\n",
    "y = np.ones((1,3)) # row vector\n",
    "print('y ='); print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Row/column vectors are like row/column matrix. We have to use two indices to access the elements of such row/column vectors.\n",
    "\n",
    "Here we access elements of a column vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "print(x[0][0])\n",
    "print(x[1][0])\n",
    "print(x[2][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`x[0]` gives an array for the first row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First row of x       = [1.] <class 'numpy.ndarray'>\n",
      "Element in first row = 1.0 <class 'numpy.float64'>\n"
     ]
    }
   ],
   "source": [
    "print('First row of x       =',x[0], type(x[0]))\n",
    "print('Element in first row =',x[0][0], type(x[0][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing portions of arrays\n",
    "Array of 10 elements\n",
    "\n",
    "| x[0] | x[1] | x[2] | x[3] | x[4] | x[5] | x[6] | x[7] | x[8] | x[9] |\n",
    "|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|\n",
    "| 0    | 1    | 2    | 3    | 4    | 5    | 6    | 7    | 8    | 9    |\n",
    "|x[-10]| x[-9]| x[-8]| x[-7]| x[-6]| x[-5]| x[-4]| x[-3]| x[-2]|x[-1] |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(0,9,10)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get elements ```x[2],...,x[5]```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2. 3. 4. 5.]\n"
     ]
    }
   ],
   "source": [
    "print(x[2:6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence ```x[m:n]``` gives the elements ```x[m],x[m+1],...,x[n-1]```.\n",
    "\n",
    "Get elements from ```x[5]``` upto the last element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5. 6. 7. 8. 9.]\n"
     ]
    }
   ],
   "source": [
    "print(x[5:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the last element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.0\n"
     ]
    }
   ],
   "source": [
    "print(x[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get element ```x[5]``` upto last but one element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5. 6. 7. 8.]\n"
     ]
    }
   ],
   "source": [
    "print(x[5:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Access every alternate element of array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 2. 4. 6. 8.]\n"
     ]
    }
   ],
   "source": [
    "print(x[0::2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 3. 5. 7. 9.]\n"
     ]
    }
   ],
   "source": [
    "print(x[1::2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These operations work on multi dimensional arrays also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[5.67717102e-01 2.94552656e-01 1.76235501e-04 3.49796707e-01]\n",
      " [5.60834819e-01 7.28843246e-04 8.59005230e-01 5.02693996e-01]\n",
      " [6.32336995e-01 6.88768612e-01 6.63854136e-01 5.80712788e-01]]\n"
     ]
    }
   ],
   "source": [
    "A = np.random.rand(3,4)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.67717102e-01 2.94552656e-01 1.76235501e-04 3.49796707e-01]\n"
     ]
    }
   ],
   "source": [
    "print(A[0,:]) # 0'th row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5677171  0.56083482 0.632337  ]\n"
     ]
    }
   ],
   "source": [
    "print(A[:,0]) # 0'th column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[5.67717102e-01 2.94552656e-01 1.76235501e-04]\n",
      " [5.60834819e-01 7.28843246e-04 8.59005230e-01]]\n"
     ]
    }
   ],
   "source": [
    "print(A[0:2,0:3]) # print submatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
      " [5.60834819e-01 7.28843246e-04 8.59005230e-01 5.02693996e-01]\n",
      " [6.32336995e-01 6.88768612e-01 6.63854136e-01 5.80712788e-01]]\n"
     ]
    }
   ],
   "source": [
    "A[0,:] = 0.0 # zero out zeroth row\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arithmetic operations on arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arithmetic operations act element-wise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4. 10. 18.]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1.0, 2.0, 3.0])\n",
    "y = np.array([4.0, 5.0, 6.0])\n",
    "print(x*y)  # multiply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.25 0.4  0.5 ]\n"
     ]
    }
   ],
   "source": [
    "print(x/y)  # divide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  4.  25. 216.]\n"
     ]
    }
   ],
   "source": [
    "print(y**x) # exponentiation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2. 3.]\n",
      " [1. 2. 3.]\n",
      " [1. 2. 3.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.ones((3,3))\n",
    "print(A*x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If ```A``` and ```x``` are arrays, then ```A*x``` does not give matrix-vector product. For that use ```dot```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6. 6. 6.]\n"
     ]
    }
   ],
   "source": [
    "print(A.dot(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or equivalently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6. 6. 6.]\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(A,x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Python3 versions, we can use `@` to achieve matrix operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6. 6. 6.]\n"
     ]
    }
   ],
   "source": [
    "print(A@x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can of course do matrix-matrix products using ```dot``` or ```@```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      " [[1. 1. 1.]\n",
      " [1. 1. 1.]\n",
      " [1. 1. 1.]]\n",
      "B =\n",
      " [[2. 2. 2.]\n",
      " [2. 2. 2.]\n",
      " [2. 2. 2.]]\n",
      "A*B =\n",
      " [[6. 6. 6.]\n",
      " [6. 6. 6.]\n",
      " [6. 6. 6.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.ones((3,3))\n",
    "B = 2*A\n",
    "print('A =\\n',A)\n",
    "print('B =\\n',B)\n",
    "print('A*B =\\n',A@B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On vectors, `@` performs dot product, i.e., `x@y = dot(x,y)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1,1,1])\n",
    "y = np.array([2,3,4])\n",
    "print(x@y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some matrix/vector functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min     =  -3\n",
      "max     =  3\n",
      "abs min =  0\n",
      "abs max =  3\n",
      "sum     =  0\n",
      "dot     =  28\n",
      "dot     =  28\n"
     ]
    }
   ],
   "source": [
    "x = np.array([-3,-2,-1,0,1,2,3])\n",
    "print('min     = ',np.min(x))\n",
    "print('max     = ',np.max(x))\n",
    "print('abs min = ',np.abs(x).min()) # np.min(np.abs(x))\n",
    "print('abs max = ',np.abs(x).max()) # np.max(np.abs(x))\n",
    "print('sum     = ',np.sum(x))\n",
    "print('dot     = ',np.dot(x,x)) \n",
    "print('dot     = ',np.sum(x*x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute vector norms using [numpy.linalg.norm](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html) (also see [scipy.linalg.norm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.norm.html))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.291502622129181\n",
      "5.291502622129181\n",
      "12.0\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "from numpy.linalg import norm\n",
    "print(norm(x))        # L2 norm\n",
    "print(norm(x,2))      # L2 norm\n",
    "print(norm(x,1))      # L1 norm\n",
    "print(norm(x,np.inf)) # Linf norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or import the `linalg` module and use methods inside it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.291502622129181\n",
      "5.291502622129181\n",
      "12.0\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "import numpy.linalg as la\n",
    "print(la.norm(x))        # L2 norm\n",
    "print(la.norm(x,2))      # L2 norm\n",
    "print(la.norm(x,1))      # L1 norm\n",
    "print(la.norm(x,np.inf)) # Linf norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These methods work on 2-d arrays which are square."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.40159839 -0.8591952   0.16223542]\n",
      " [-0.9098026  -0.51844933 -0.04381662]\n",
      " [-0.22818608 -0.20678772 -0.81851499]]\n",
      "min       =  -0.9098026020175445\n",
      "max       =  0.4015983909742702\n",
      "abs min   =  0.04381661776429091\n",
      "abs max   =  0.9098026020175445\n",
      "sum       =  -3.020918714630895\n",
      "Frob norm =  1.6700494585979953\n",
      "L1 norm   =  1.5844322417075245\n",
      "L2 norm   =  1.1504001789154477\n",
      "Linf norm =  1.4720685462097143\n"
     ]
    }
   ],
   "source": [
    "A = 2*np.random.rand(3,3) - 1\n",
    "print(A)\n",
    "print('min       = ',np.min(A))\n",
    "print('max       = ',np.max(A))\n",
    "print('abs min   = ',np.abs(A).min()) # np.min(np.abs(x))\n",
    "print('abs max   = ',np.abs(A).max()) # np.max(np.abs(x))\n",
    "print('sum       = ',np.sum(A))\n",
    "print('Frob norm = ',la.norm(A))        # Frobenius norm\n",
    "print('L1 norm   = ',la.norm(A,1))      # L1 norm\n",
    "print('L2 norm   = ',la.norm(A,2))      # L2 norm\n",
    "print('Linf norm = ',la.norm(A,np.inf)) # Linf norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Matrix-vector product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let\n",
    "$$\n",
    "x, y \\in \\re^n, \\qquad A \\in \\re^{n \\times n}\n",
    "$$\n",
    "To compute the matrix-vector product $y=Ax$, we can do it element-wise\n",
    "$$\n",
    "y_i = \\sum_{j=0}^{n-1} A_{ij} x_j, \\qquad 0 \\le i \\le n-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "x = np.random.rand(n)\n",
    "A = np.random.rand(n,n)\n",
    "y = np.zeros(n)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        y[i] += A[i,j]*x[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify that our result is correct by computing $\\| y - A x\\|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.108895957933346e-16\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.norm(y-A@x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute the product column-wise. Let\n",
    "$$\n",
    "A_{:,j} = \\textrm{j'th column of A}\n",
    "$$\n",
    "Then the matrix-vector product can also be written as\n",
    "$$\n",
    "y = \\sum_{j=0}^{n-1} A_{:,j} x_j\n",
    "$$\n",
    "**Warning**: This may have inefficient memory access since by default, numpy arrays have column-major ordering, see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.108895957933346e-16\n"
     ]
    }
   ],
   "source": [
    "y[:] = 0.0\n",
    "for j in range(n):\n",
    "    y += A[:,j]*x[j]\n",
    "\n",
    "# Now check the result\n",
    "print(np.linalg.norm(y-A@x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Matrix-Matrix product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $A \\in \\re^{m\\times n}$ and $B \\in \\re^{n \\times p}$ then $C = AB \\in \\re^{m \\times p}$ is given by\n",
    "$$\n",
    "C_{ij} = \\sum_{k=0}^{n-1} A_{ik} B_{kj}, \\qquad 0 \\le i \\le m-1, \\quad 0 \\le j \\le p-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n,p = 10,8,6\n",
    "A = np.random.rand(m,n)\n",
    "B = np.random.rand(n,p)\n",
    "C = np.zeros((m,p))\n",
    "for i in range(m):\n",
    "    for j in range(p):\n",
    "        for k in range(n):\n",
    "            C[i,j] += A[i,k]*B[k,j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us verify the result is correct by computing the Frobenius norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4770556661825108e-15\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.norm(C - A@B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another view-point is the following\n",
    "$$\n",
    "C_{ij} = (\\textrm{i'th row of A}) \\cdot (\\textrm{j'th column of B})\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "for i in range(m):\n",
    "    for j in range(p):\n",
    "        C[i,j] = A[i,:] @ B[:,j]\n",
    "\n",
    "# Now check the result\n",
    "print(np.linalg.norm(C - A@B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outer product of two vectors\n",
    "\n",
    "Given $m$-vector $a$ and $n$-vector $b$, both considered as column vectors, their outer product\n",
    "\n",
    "$$\n",
    "A = a b^\\top = \\begin{bmatrix} \n",
    "|     &  |    &        & | \\\\\n",
    "a b_0 & a b_1 & \\ldots & a b_{n-1} \\\\\n",
    "|     &  |    &        & | \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "is an $m \\times n$ matrix with elements\n",
    "\n",
    "$$\n",
    "A_{ij} = a_j b_j, \\qquad 0 \\le i \\le m-1, \\quad 0 \\le n \\le n-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 2)\n",
      "[[1 2]\n",
      " [2 4]\n",
      " [3 6]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = np.array([1,2])\n",
    "A = np.outer(a,b)\n",
    "print(A.shape)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Math functions\n",
    "\n",
    "Numpy provides standard functions like `sin`, `cos`, `log`, etc. which can act on arrays in an element-by-element manner. This is not the case for functions in `math` module, which can only take scalar arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [0.         1.57079633 3.14159265 4.71238898 6.28318531]\n",
      "y = [ 0.0000000e+00  1.0000000e+00  1.2246468e-16 -1.0000000e+00\n",
      " -2.4492936e-16]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(0.0, 2.0*np.pi, 5)\n",
    "y = np.sin(x)\n",
    "print('x =',x)\n",
    "print('y =',y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      " [[0.69799961 0.44193454 0.42009594 0.59430268 0.44412881]\n",
      " [0.88453597 0.60609966 0.57272725 0.68715342 0.16661875]\n",
      " [0.5920093  0.63078135 0.78454055 0.81219262 0.73284029]\n",
      " [0.27186644 0.51657847 0.28403543 0.18333706 0.2893531 ]\n",
      " [0.59918434 0.9679843  0.18048533 0.82641006 0.26114357]]\n",
      "Y =\n",
      " [[0.64268642 0.42768894 0.40784805 0.55993113 0.42967137]\n",
      " [0.77362104 0.5696662  0.54192611 0.63433919 0.16584888]\n",
      " [0.5580295  0.58977593 0.7065001  0.72579724 0.66898346]\n",
      " [0.26852979 0.49390796 0.28023166 0.18231171 0.28533228]\n",
      " [0.56396909 0.82374457 0.17950705 0.73550386 0.25818551]]\n"
     ]
    }
   ],
   "source": [
    "A = np.random.rand(5,5) # 5x5 matrix\n",
    "Y = np.sin(A)\n",
    "print('A =\\n',A)\n",
    "print('Y =\\n',Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Memory ordering in arrays*\n",
    "By default, the ordering is same as in C/C++, the inner-most index is the fastest running one. For example, if we have an array of size (2,3), they are stored in memory in this order\n",
    "```\n",
    "a[0,0], a[0,1], a[0,2], a[1,0], a[1,1], a[1,2]\n",
    "```\n",
    "i.e.,\n",
    "```\n",
    "a[0,0] --> a[0,1] --> a[0,2]\n",
    "                         |\n",
    "   |----------<----------|\n",
    "   |\n",
    "a[1,0] --> a[1,1] --> a[1,2]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Row contiguous = True\n",
      "Col contiguous = False\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  C_CONTIGUOUS : True\n",
       "  F_CONTIGUOUS : False\n",
       "  OWNDATA : True\n",
       "  WRITEABLE : True\n",
       "  ALIGNED : True\n",
       "  WRITEBACKIFCOPY : False"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([[1,2,3], [4,5,6]])\n",
    "print('Row contiguous =', a[0,:].data.contiguous)\n",
    "print('Col contiguous =', a[:,0].data.contiguous)\n",
    "a.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get fortran style ordering, where the outer-most index is the fastest running one, which corresponds to the following layout in memory\n",
    "```\n",
    "a[0,0], a[1,0], a[0,1], a[1,1], a[0,2], a[1,2]\n",
    "```\n",
    "i.e.,\n",
    "```\n",
    "a[0,0]    -- a[0,1]    -- a[0,2]\n",
    "   |     |     |      |      |  \n",
    "   v     ^     v      ^      v\n",
    "   |     |     |      |      |\n",
    "a[1,0] --    a[1,1] --    a[1,2]\n",
    "```\n",
    "create like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Row contiguous = False\n",
      "Col contiguous = True\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  C_CONTIGUOUS : False\n",
       "  F_CONTIGUOUS : True\n",
       "  OWNDATA : True\n",
       "  WRITEABLE : True\n",
       "  ALIGNED : True\n",
       "  WRITEBACKIFCOPY : False"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array([[1,2,3], [4,5,6]], order='F')\n",
    "print('Row contiguous =', b[0,:].data.contiguous)\n",
    "print('Col contiguous =', b[:,0].data.contiguous)\n",
    "b.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tensor product array: meshgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "len(x) =  4\n",
      "len(y) =  3\n",
      "shape X=  (3, 4)\n",
      "shape Y=  (3, 4)\n",
      "x =  [0. 1. 2. 3.]\n",
      "y =  [0. 1. 2.]\n",
      "X = \n",
      "[[0. 1. 2. 3.]\n",
      " [0. 1. 2. 3.]\n",
      " [0. 1. 2. 3.]]\n",
      "Y = \n",
      "[[0. 0. 0. 0.]\n",
      " [1. 1. 1. 1.]\n",
      " [2. 2. 2. 2.]]\n"
     ]
    }
   ],
   "source": [
    "x = np.linspace(0,3,4)\n",
    "y = np.linspace(0,2,3)\n",
    "X, Y = np.meshgrid(x,y)\n",
    "print('len(x) = ',len(x))\n",
    "print('len(y) = ',len(y))\n",
    "print('shape X= ',X.shape)\n",
    "print('shape Y= ',Y.shape)\n",
    "print('x = ', x)\n",
    "print('y = ', y)\n",
    "print('X = '); print(X)\n",
    "print('Y = '); print(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $x \\in \\re^m$ and $y \\in \\re^n$, then the output is arranged like this\n",
    "$$\n",
    "X[i,j] = x[j], \\qquad\n",
    "Y[i,j] = y[i], \\qquad 0 \\le i \\le n-1, \\quad 0 \\le j \\le m-1\n",
    "$$\n",
    "If we want the following arrangement\n",
    "$$\n",
    "X[i,j] = x[i], \\qquad\n",
    "Y[i,j] = y[j], \\qquad 0 \\le i \\le m-1, \\quad 0 \\le j \\le n-1\n",
    "$$\n",
    "we have to do the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "len(x) =  4\n",
      "len(y) =  3\n",
      "shape X=  (4, 3)\n",
      "shape Y=  (4, 3)\n"
     ]
    }
   ],
   "source": [
    "X, Y = np.meshgrid(x,y,indexing='ij')\n",
    "print('len(x) = ',len(x))\n",
    "print('len(y) = ',len(y))\n",
    "print('shape X= ',X.shape)\n",
    "print('shape Y= ',Y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second form is useful when working with finite difference schemes on Cartesian grids, where we want to use i index running along x-axis and j index running along y-axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshaping arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2,3],\n",
    "              [4,5,6]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4 5 6]\n"
     ]
    }
   ],
   "source": [
    "B = np.reshape(A,2*3,order='C')\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "A1 = np.reshape(B,(2,3),order='C')\n",
    "print(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 4 2 5 3 6]\n"
     ]
    }
   ],
   "source": [
    "C = np.reshape(A,2*3,order='F')\n",
    "print(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "A2 = np.reshape(C,(2,3),order='F')\n",
    "print(A2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing and reading files\n",
    "Write an array to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.23800817 0.96267486 0.42996244]\n",
      " [0.31093536 0.61795882 0.65229253]\n",
      " [0.44645543 0.3532086  0.88249008]]\n"
     ]
    }
   ],
   "source": [
    "A = np.random.rand(3,3)\n",
    "print(A)\n",
    "np.savetxt('A.txt',A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the contents of the file in your terminal\n",
    "```\n",
    "cat A.txt\n",
    "```\n",
    "We can also print it from within the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.380081667929159206e-01 9.626748557062282385e-01 4.299624366157210886e-01\n",
      "3.109353645036705416e-01 6.179588184832025544e-01 6.522925255033611425e-01\n",
      "4.464554342116026087e-01 3.532086017204271178e-01 8.824900753021190924e-01\n"
     ]
    }
   ],
   "source": [
    "!cat A.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write two 1-D arrays as columns into file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000000000000000000e+00 2.000000000000000000e+00\n",
      "2.000000000000000000e+00 4.000000000000000000e+00\n",
      "3.000000000000000000e+00 6.000000000000000000e+00\n",
      "4.000000000000000000e+00 8.000000000000000000e+00\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1.0,2.0,3.0,4.0])\n",
    "y = np.array([2.0,4.0,6.0,8.0])\n",
    "np.savetxt('data.txt',np.column_stack([x,y]))\n",
    "!cat data.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can control the number of decimals saved, and use scientific notation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000e+00 2.0000e+00\n",
      "2.0000e+00 4.0000e+00\n",
      "3.0000e+00 6.0000e+00\n",
      "4.0000e+00 8.0000e+00\n"
     ]
    }
   ],
   "source": [
    "np.savetxt('data.txt',np.column_stack([x,y]),fmt='%8.4e')\n",
    "!cat data.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can read an existing file like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape d = (4, 2)\n",
      "x = [1. 2. 3. 4.]\n",
      "y = [2. 4. 6. 8.]\n"
     ]
    }
   ],
   "source": [
    "d = np.loadtxt('data.txt')\n",
    "print('Shape d =', d.shape)\n",
    "x1 = d[:,0]\n",
    "y1 = d[:,1]\n",
    "print('x =',x1)\n",
    "print('y =',y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measuring time\n",
    "\n",
    "Loops are slow in Python and it is good to avoid them. The following example show how to time code execution.\n",
    "\n",
    "Create some random matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 1000\n",
    "A = np.random.random((m,m))\n",
    "B = np.random.random((m,m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we add two matrices with an explicit for loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "230 ms ± 5.1 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r10 -n10\n",
    "C = np.empty_like(A)\n",
    "for i in range(m):\n",
    "    for j in range(m):\n",
    "        C[i,j] = A[i,j] + B[i,j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use the + operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1 ms ± 362 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r10 -n10\n",
    "C = A + B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loop version is significantly slower."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "c6e4e9f98eb68ad3b7c296f83d20e6de614cb42e90992a65aa266555a3137d0d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}