{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy Array Operations and Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** We will see more tools for linear algebra in the section on [Scipy Tools](scipy-tools-for-linear-algebra.ipynb), which introduces the package *Scipy*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As of version 3.5, Python can handle matrix multiplication on Numpy arrays, using the at sign \"@\":\n",
    "\n",
    "    C = A @ B\n",
    "    \n",
    "Why this strange notation? Because\n",
    "\n",
    "    D = A * B\n",
    "\n",
    "is instead \"point-wise\" multiplication of arrays, with `D[i,j] = A[i,j] * B[i,j]` instead of\n",
    "\n",
    "$$C_{i,j}= \\sum_{k} a_{ik} b_{kj}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Aside 1: Numpy `matrix` is now redundant!**\n",
    "\n",
    "Previously, Numpy supported matrix calculations withe variables of type `matrix`, a special sub-tpye of `array`.\n",
    "However, with the addition of \"@\" for matrix multiplication of arrays, that option is now un-needed, and I advise you to avoid it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Aside 2: Floating point numbers vs integers**\n",
    "\n",
    "Though Python is generally good at understanding when an integer like `7` is to be used as a floating point (real) number, it is sometimes best to make this distinction explicitly when working with module `numpy`;\n",
    "otherwise sometimes division done within numpy functions returns an integer answer, like `7/2 = 3`.\n",
    "\n",
    "Thus from now on, when I mean an integer to be used as a floating point number, I give it a decimal point: `7./2.` will reliably be `3.5`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we explore some basic features offered by Numpy; then we will look at some more advanced tools from package Scipy, and specifically its module `scipy.linalg` for linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "B is\n",
      "[[1 2]\n",
      " [3 4]\n",
      " [5 6]]\n",
      "C is\n",
      "[[10  9  8]\n",
      " [ 7  6  5]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3],[4, 5, 6]])\n",
    "B = np.array([[1, 2], [3,4], [5, 6]])\n",
    "C = np.array([[10, 9, 8],[7, 6, 5]])\n",
    "\n",
    "print(f\"A is\\n{A}\")\n",
    "print(f\"B is\\n{B}\")\n",
    "print(f\"C is\\n{C}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The matrix product A times B is:\n",
      " [[22 28]\n",
      " [49 64]]\n",
      "The matrix product B times A is:\n",
      " [[ 9 12 15]\n",
      " [19 26 33]\n",
      " [29 40 51]]\n"
     ]
    }
   ],
   "source": [
    "print(f\"The matrix product A times B is:\\n {A @ B}\")\n",
    "print(f\"The matrix product B times A is:\\n {B @ A}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (2,3) (3,2) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-3bb26336abf9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"The array product A * B fails:\\n{A * B}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (2,3) (3,2) "
     ]
    }
   ],
   "source": [
    "print(f\"The array product A * B fails:\\n{A * B}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Array product A times C is\n",
      "[[10 18 24]\n",
      " [28 30 30]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f\"Array product A times C is\\n{A * C}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-55b834c682c2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"Matrix product A times C fails:\\n{A @ C}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)"
     ]
    }
   ],
   "source": [
    "print(f\"Matrix product A times C fails:\\n{A @ C}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The matrix transpose $A^T$ is given by `A.T`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A-transpose is\n",
      "[[1 4]\n",
      " [2 5]\n",
      " [3 6]]\n"
     ]
    }
   ],
   "source": [
    "print(f\"A-transpose is\\n{A.T}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Slicing*: Extracting rows, columns, and other rectangular chunks from matrices\n",
    "\n",
    "This works with Python lists and Numpy arrays, and we have seen some of it before;\n",
    "I review it here because it will help with doing the row operations of linear algebra.\n",
    "\n",
    "### Index notation for slicing\n",
    "\n",
    "For an index with n possible values, from 0 to n-1:\n",
    "\n",
    "- `a:b` means indices i in th usual semi-open interval $a \\leq i < b$ or $[a,b)$\n",
    "\n",
    "- `a:` is short for `a:n`, so indices $a \\leq i$, all the way to the maximum index value \n",
    "- `:b` is short for `0:b`, so all indices $i < b$\n",
    "- `:` combines both of the above, so is short for `0:n`; all possible indices\n",
    "- index value `-1` refers the last entry; the same as index n-1, but without needing to know n.\n",
    "- index value `-m` refers the \"m-th last\" entry; the same as index n-m, again without needing to know n."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "The column of index 1 (presented as a row vector): [2 5]\n",
      "The row of index 1: [4 5 6]\n",
      "The first 2 elements of the row of index 1: [4 5]\n",
      "Another way to say the above: [4 5]\n",
      "The bottom-right element: 6\n",
      "The 2x2 sub-matrix in the bottom-right corner:\n",
      "[[2 3]\n",
      " [5 6]]\n"
     ]
    }
   ],
   "source": [
    "print(f\"A:\\n{A}\")\n",
    "print(f\"The column of index 1 (presented as a row vector): {A[:,1]}\")\n",
    "print(f\"The row of index 1: {A[1,:]}\")\n",
    "print(f\"The first 2 elements of the row of index 1: {A[1,:2]}\")\n",
    "print(f\"Another way to say the above: {A[1,0:2]}\")\n",
    "print(f\"The bottom-right element: {A[-1,-1]}\")\n",
    "print(f\"The 2x2 sub-matrix in the bottom-right corner:\\n{A[-2:,-2:]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synonyms for array names with the equal sign (not copying!)\n",
    "\n",
    "If we use the equal sign between two array names, that makes them synonyms, referring to the same values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "Anickname is\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "Anickname is now\n",
      "[[12  2  3]\n",
      " [ 4  5  6]]\n",
      "and so is A!:\n",
      "[[12  2  3]\n",
      " [ 4  5  6]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3],[4, 5, 6]])\n",
    "print(f\"A is\\n{A}\")\n",
    "Anickname = A\n",
    "print(f\"Anickname is\\n{Anickname}\")\n",
    "Anickname[0,0] = 12\n",
    "print(f\"Anickname is now\\n{Anickname}\")\n",
    "print(f\"and so is A!:\\n{A}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copying arrays with method .copy()\n",
    "\n",
    "Thus if we want a separate new array or list with the same elements initially, we must make a copy with the method `.copy()`, not the equal sign:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "Acopy is\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n",
      "Acopy is now\n",
      "[[54  2  3]\n",
      " [ 4  5  6]]\n",
      "A is still\n",
      "[[1 2 3]\n",
      " [4 5 6]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2, 3],[4, 5, 6]])\n",
    "print(f\"A is\\n{A}\")\n",
    "Acopy = A.copy()\n",
    "print(f\"Acopy is\\n{Acopy}\")\n",
    "Acopy[0,0] = 54\n",
    "print(f\"Acopy is now\\n{Acopy}\")\n",
    "print(f\"A is still\\n{A}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=Exercise-A></a>\n",
    "#### Exercise A\n",
    "\n",
    "Create a Numpy array (not a Numpy *matrix*; those are now mostly obsolete!) containing the matrix\n",
    "\n",
    "$$A = \\left[ \\begin{array}{ccc} 4. & 2. & 1. \\\\ 9. & 3. & 1. \\\\ 25. & 5. & 1. \\end{array} \\right]$$\n",
    "\n",
    "and one containing the vector\n",
    "\n",
    "$$b = [ 0.693147, 1.098612, 1.386294]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=Exercise-B></a>\n",
    "#### Exercise B\n",
    "\n",
    "Create arrays $c$ and $d$ containing respectively the last row of $A$ and the middle column of $A$.\n",
    "\n",
    "**Note:** Do this by manipulating the array `A` with indexing and slicing operations, without entering any numerical values for array entries explicity."
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}