{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# \"Lab:\" *Solving Linear Systems Using Row Reductions*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This class requires NO CODING, i.e. this kind of lab nonsense will not show up on the exams.  However, LIFE as an applied mathematician, engineer, computer scientist, etc... tends to require coding skills.\n",
    "\n",
    "Hence, the **Purpose** of these \"labs\" are\n",
    "* To **provide Linear Algebra examples** that are modifiable, and potentially, interactive\n",
    "* To **expose** those who are interested to some python code examples;\n",
    "\n",
    "The **\"Rules\"**\n",
    "* The labs are are highly experimental and are for *Adults and Entertainment Purposes Only*\n",
    "* The code is sometimes (always?) crude, and should not be taken an example of Good Coding Standards[TM]\n",
    "* The labs \"live\" at http://terminus.sdsu.edu/SDSU/Math254/PythonLabs/\n",
    "* The labs are developed for *Math 254 : Introduction to Linear Algebra* at San Diego State University\n",
    "* The labs are provided free of cost, with a full money-back guarantee\n",
    "\n",
    "**Copyright, and License**\n",
    "* Copyright (C) 2019 : Peter Blomgren <blomgren.peter@gmail.com>.\n",
    "* License: [GNU General Public License, version 3](https://en.wikipedia.org/wiki/GNU_General_Public_License)\n",
    "\n",
    "**Interactive or Static Mode???**\n",
    "* If you are using the static link which means you are using https://nbviewer.jupyter.org/, you are looking at the notebook Inputs and Outputs, and you cannot modify or interact in any way with the notebook.\n",
    "* If you have downloaded the notebook and are viewing it some other way --- $(\\infty-1)$ possibilities, you can step through each code-block one-step-at-a-time by pressing [shift]+[enter], if you want to do this without installing anything on your computer one option is to use [mybinder](https://mybinder.org/v2/gh/ipython/ipython-in-depth/master?filepath=binder/Index.ipynb) (surely, this will break at some point), or [Google/colab](https://colab.research.google.com)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *(In interactive mode)* Press [shift]+[enter] to step through the blocks one at a time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import basic python packages for numerical computation (numpy),\n",
    "# symbolic computations (sympy), and plotting (matplotlib)\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# cm are colormaps -- for plotting\n",
    "from matplotlib import cm\n",
    "\n",
    "# mplot3d is needed for 3D-plotting\n",
    "from mpl_toolkits import mplot3d\n",
    "\n",
    "#(not needed in this lab) from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
    "#(not needed in this lab) from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *\"What, is like, your problem... maaaan?!?\"*\n",
    "Our problem of interest is\n",
    "$\\begin{bmatrix}\n",
    "2x &+& 8y &+& 4z \\\\\n",
    "2x &+& 5y &+& z \\\\\n",
    "4x &+& 10y &=& z\n",
    "\\end{bmatrix} = \\begin{bmatrix} 2 \\\\ 5 \\\\ 1\\end{bmatrix}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 8 & 4\\\\2 & 5 & 1\\\\4 & 10 & -1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2,  8,  4],\n",
       "[2,  5,  1],\n",
       "[4, 10, -1]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We can use \"Symbolic\" matrices by invoking the Matrix command from sympy\n",
    "# \"Symbolic\" matrices are computationally inefficient, but are useful for\n",
    "# teaching and learning linear algebra on small matrices.\n",
    "\n",
    "# First, we let A hold the coeffient matrix\n",
    "A = sp.Matrix([[2, 8, 4], [2, 5, 1], [4, 10, -1]])\n",
    "\n",
    "# ... and look at how pretty the output is, without any extra effort\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2\\\\5\\\\1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2],\n",
       "[5],\n",
       "[1]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Next, we define the right-hand-side of the equation\n",
    "#\n",
    "# When you say \"linear system\", linear algebraists think \"Ax=b\",\n",
    "# hence\n",
    "b = sp.Matrix([[2],[5],[1]])\n",
    "b\n",
    "# ... and, yeah, a column vector is really just a \"skinny matrix\" (at least in Python)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 8 & 4\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[2, 8, 4]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Note: like many computing environments, Python counts from 0\n",
    "#\n",
    "# This is the row of A\n",
    "A[0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}8\\\\5\\\\10\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ 8],\n",
       "[ 5],\n",
       "[10]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and this is the 2nd row\n",
    "A[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 8 & 4 & 2\\\\2 & 5 & 1 & 5\\\\4 & 10 & -1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2,  8,  4, 2],\n",
       "[2,  5,  1, 5],\n",
       "[4, 10, -1, 1]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Lets create a of A, and then augment it with the right-hand-side\n",
    "B = A\n",
    "B = B.row_join(b)\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4 & 2 & 1\\\\2 & 5 & 1 & 5\\\\4 & 10 & -1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,  4,  2, 1],\n",
       "[2,  5,  1, 5],\n",
       "[4, 10, -1, 1]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# First we create a leading 1, in the first-row first-column\n",
    "B[0,:] /= B[0,0]\n",
    "# Since B[0,0] started out with the value 2, this exactly the operation preformed in the [Notes 1.2, slide 14]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4 & 2 & 1\\\\0 & -3 & -3 & 3\\\\4 & 10 & -1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,  4,  2, 1],\n",
       "[0, -3, -3, 3],\n",
       "[4, 10, -1, 1]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Next, we subtract appropriate multipliers of the first row from the second and third:\n",
    "B[1,:] = B[1,:] - B[1,0] * B[0,:]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4 & 2 & 1\\\\0 & -3 & -3 & 3\\\\0 & -6 & -9 & -3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,  4,  2,  1],\n",
       "[0, -3, -3,  3],\n",
       "[0, -6, -9, -3]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Next, we subtract appropriate multipliers of the first row from the second and third:\n",
    "B[2,:] = B[2,:] - B[2,0] * B[0,:]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 4 & 2 & 1\\\\0 & 1 & 1 & -1\\\\0 & -6 & -9 & -3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1,  4,  2,  1],\n",
       "[0,  1,  1, -1],\n",
       "[0, -6, -9, -3]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We now move to the second row, and make the leading coefficient 1\n",
    "B[1,:] /= B[1,1]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & -2 & 5\\\\0 & 1 & 1 & -1\\\\0 & 0 & -3 & -9\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1, 0, -2,  5],\n",
       "[0, 1,  1, -1],\n",
       "[0, 0, -3, -9]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's eliminate \"up\" and \"down\" in the 2nd column\n",
    "B[0,:] = B[0,:] - B[0,1] * B[1,:]\n",
    "B[2,:] = B[2,:] - B[2,1] * B[1,:]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & -2 & 5\\\\0 & 1 & 1 & -1\\\\0 & 0 & 1 & 3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1, 0, -2,  5],\n",
       "[0, 1,  1, -1],\n",
       "[0, 0,  1,  3]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We move to the 3rd row, and the strategy is the same: make a 1, then make 0s...\n",
    "B[2,:] /= B[2,2]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 11\\\\0 & 1 & 0 & -4\\\\0 & 0 & 1 & 3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1, 0, 0, 11],\n",
       "[0, 1, 0, -4],\n",
       "[0, 0, 1,  3]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[0,:] = B[0,:] - B[0,2] * B[2,:]\n",
    "B[1,:] = B[1,:] - B[1,2] * B[2,:]\n",
    "B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can identify the solution: $$\\begin{bmatrix} x \\\\ y \\\\ z\\end{bmatrix} = \\begin{bmatrix} 11 \\\\ -4 \\\\ 3 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The steps are so predictable we can write ourselves a little function to perform all the steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WARNING: This function will run into trouble if at any point we have\n",
    "# WARNING: a zero in the current diagonal element M[i,i], since \n",
    "# WARNING: division by zero is a Very Bad Idea!!!\n",
    "#\n",
    "# WARNING: Division-by-zero is NOT the only problem this can run into.\n",
    "#\n",
    "def my_rref( M ):\n",
    "    rows, cols = M.shape\n",
    "    for i in range(0,rows):\n",
    "        M[i,:] /= M[i,i] \n",
    "        print(M)\n",
    "        for j in range(0,rows):\n",
    "            if i != j:\n",
    "                M[j,:] = M[j,:] - M[j,i] * M[i,:]\n",
    "                print(M)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 8 & 4 & 2\\\\2 & 5 & 1 & 5\\\\4 & 10 & -1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[2,  8,  4, 2],\n",
       "[2,  5,  1, 5],\n",
       "[4, 10, -1, 1]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We reconstruct the non-reduced augmented matrix\n",
    "B = A\n",
    "B = B.row_join(b)\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[1, 4, 2, 1], [2, 5, 1, 5], [4, 10, -1, 1]])\n",
      "Matrix([[1, 4, 2, 1], [0, -3, -3, 3], [4, 10, -1, 1]])\n",
      "Matrix([[1, 4, 2, 1], [0, -3, -3, 3], [0, -6, -9, -3]])\n",
      "Matrix([[1, 4, 2, 1], [0, 1, 1, -1], [0, -6, -9, -3]])\n",
      "Matrix([[1, 0, -2, 5], [0, 1, 1, -1], [0, -6, -9, -3]])\n",
      "Matrix([[1, 0, -2, 5], [0, 1, 1, -1], [0, 0, -3, -9]])\n",
      "Matrix([[1, 0, -2, 5], [0, 1, 1, -1], [0, 0, 1, 3]])\n",
      "Matrix([[1, 0, 0, 11], [0, 1, 1, -1], [0, 0, 1, 3]])\n",
      "Matrix([[1, 0, 0, 11], [0, 1, 0, -4], [0, 0, 1, 3]])\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 11\\\\0 & 1 & 0 & -4\\\\0 & 0 & 1 & 3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1, 0, 0, 11],\n",
       "[0, 1, 0, -4],\n",
       "[0, 0, 1,  3]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Then we invoke our magic, uh, function:\n",
    "rref_B = my_rref(B)\n",
    "rref_B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is pretty nice!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}