{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Matrices, Inverses, Equations, Solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Today I'll be using lists as I am still not sure if numpy and scipy are installed in the lab."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Matrices, reading from a file."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It will be painful to input a 16x16 matrix whenever we need to. So we shall load matrices from files. We shall use the split function (*method*) we learnt before. We shall write a function which takes the name of file as input and returns a list of lists of numbers."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def read_matrix(filename) :\n",
      "    fl = open(filename, 'r')\n",
      "    matrix = []\n",
      "    for line in fl :\n",
      "        row = []\n",
      "        words = line.split()\n",
      "        for word in words :\n",
      "            row.append(float(word))\n",
      "        matrix.append(row)\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = read_matrix(\"files/A.txt\")\n",
      "print A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[1.0, 2.0, 3.0], [2.0, 1.0, 2.0], [3.0, 2.0, 1.0]]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The output isn't really very pretty. Let us make it pretty by"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def print_matrix(matrix, myformat) :\n",
      "    for row in matrix :\n",
      "        for no in row :\n",
      "            print myformat % no,\n",
      "        print"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(A, \"%10.0f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "         1          2          3\n",
        "         2          1          2\n",
        "         3          2          1\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Coding elementary row operations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Row operation 1 : Interchanging two rows\n",
      "As one can guess, this will be a function which takes a matrix as an input, and row numbers of two rows; like `elem1(A, 0, 2)`. However, before we do that, it might be useful to define a function called `copy_row(source)`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def copy_row(src) :\n",
      "    dest = []\n",
      "    for entry in src :\n",
      "        dest.append(entry)\n",
      "    return dest\n",
      "\n",
      "def elem1(matrix, row1, row2) :\n",
      "    if row1 < 0 or row2 < 0 or row1 >= len(matrix) or row2 >= len(matrix) :\n",
      "        print \"Row out of range.\"\n",
      "    else :\n",
      "        temp_row = copy_row(matrix[row1])\n",
      "        matrix[row1] = copy_row(matrix[row2])\n",
      "        matrix[row2] = copy_row(temp_row)\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note lists are also passed by reference. So modifying matrix actually changes the actual matrix. You can see it below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(elem1(A, 1, 2), \"%2.0f\")\n",
      "print \"-------------------------\"\n",
      "print_matrix(A, \"%2.0f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1  2  3\n",
        " 3  2  1\n",
        " 2  1  2\n",
        "-------------------------\n",
        " 1  2  3\n",
        " 3  2  1\n",
        " 2  1  2\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Row operation 2 : multiplying a row by a non-zero scalar\n",
      "The functionw will take a matrix, and row number and a scalar as input. We'll also check if the scalar is non-zero."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def elem2(matrix, row, scalar) :\n",
      "    if row < 0 or row >= len(matrix) :\n",
      "        print \"row %d out of range.\" % row\n",
      "    elif scalar == 0 :\n",
      "        print \"scalar has to be non-zero.\"\n",
      "    else :\n",
      "        for i in range(len(matrix[row])) :\n",
      "            matrix[row][i] *= scalar\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(A, \"%2.0f\")\n",
      "print \"--------------------------\"\n",
      "print_matrix(elem2(A, 0, 4), \"%2.0f\")\n",
      "print \"-------------------------\"\n",
      "print_matrix(A, \"%2.0f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1  2  3\n",
        " 3  2  1\n",
        " 2  1  2\n",
        "--------------------------\n",
        " 4  8 12\n",
        " 3  2  1\n",
        " 2  1  2\n",
        "-------------------------\n",
        " 4  8 12\n",
        " 3  2  1\n",
        " 2  1  2\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####row operation 3 : replacing a row by the sum of that row and a multiple of another row.\n",
      "Here we need 4 inputs to the function:\n",
      "+ The matrix\n",
      "+ The row to be changed\n",
      "+ The row whose multiple will be added to the first row,\n",
      "+ a scalar : the multiplication factor.\n",
      "\n",
      "Again for this we use some helper functions which are readily available for vectors."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def add_vects(lst1, lst2) :\n",
      "    return [a + b for (a, b) in zip(lst1, lst2)]\n",
      "\n",
      "def scalar_mult(a, lst) :\n",
      "    return [a * no for no in lst]\n",
      "\n",
      "\n",
      "def elem3(matrix, row_2b_changed, row_used_4_change, scalar) :\n",
      "    temp_row = copy_row(scalar_mult(scalar, matrix[row_used_4_change]))\n",
      "    matrix[row_2b_changed] = add_vects(matrix[row_2b_changed], temp_row)\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(A, \"%2.0f\")\n",
      "print \"-\"*20\n",
      "print_matrix(elem3(A, 0, 1, -2), \"%2.0f\")\n",
      "print \"-\"*21\n",
      "print_matrix(A, \"%2.0f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 4  8 12\n",
        " 3  2  1\n",
        " 2  1  2\n",
        "--------------------\n",
        "-2  4 10\n",
        " 3  2  1\n",
        " 2  1  2\n",
        "---------------------\n",
        "-2  4 10\n",
        " 3  2  1\n",
        " 2  1  2\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Sweeping a column"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us read another matrix $B$. This is the same one as in the class last week."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "B = read_matrix(\"files/B.txt\")\n",
      "print_matrix(B, \"%2.0f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0  2  3  2\n",
        " 1  4  5  5\n",
        " 3  0  2  6\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we recall, sweeping a column requires a matrix and an entry (position in the matrix) which we call pivot. We shall define the pivot as the 2-tuple (row_no, col_no)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sweep(matrix, pivot) :\n",
      "    (r, c) = pivot\n",
      "    if r < 0 or r >= len(matrix) :\n",
      "        print \"row out of range\"\n",
      "    elif c < 0 or c >= len(matrix[r]) :\n",
      "        print \"column out of range\"\n",
      "    elif matrix[r][c] == 0:\n",
      "        # pivot cannot be zero.\n",
      "        print(\"pivot cannot be zero.\")\n",
      "    else :\n",
      "        # print_matrix(matrix, \"%5.2f\")\n",
      "        pivot = matrix[r][c]\n",
      "        \n",
      "        # Step 1 : multiply the row containing the pivot by 1/pivot to make the pivot 1.\n",
      "        matrix = elem2(matrix, r, 1.0/pivot)\n",
      "        # print_matrix(matrix, \"%5.2f\")\n",
      "        \n",
      "        # Step 2 : for row != that of pivot, subtract matrix[row][c] times the row containing pivot\n",
      "        for i in range(len(matrix)) :\n",
      "            if i != r :\n",
      "                matrix = elem3(matrix, i, r, -matrix[i][c])\n",
      "                # print_matrix(matrix, \"%5.2f\")\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print_matrix(sweep(B, (1, 2)), \"%6.2f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " -0.60  -0.40   0.00  -1.00\n",
        "  0.20   0.80   1.00   1.00\n",
        "  2.60  -1.60   0.00   4.00\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Gaussian elimination"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Gaussian elemination is an algorithm to reduce a matrix to its reduced echelon form. Here *reduction* means performing row operations till the final matrix satisfies the definition of a reduced echelon form. For your convenience let us recall the definition.\n",
      "\n",
      "An $m \\times n$ matrix is said to be in echelon form, if it has $r$, $0 \\leq r \\leq m$ non-zero rows and\n",
      "1. All the non-zero rows are on the top.\n",
      "1. For $1 \\leq i \\leq r$, $p_i$ denotes the column containing the *first non-zero entry* of the $i$-th row, then\n",
      "   $p_1 < p_2 < \\dotsb < p_r$.\n",
      "1. $a_{i p_i} = 1$.\n",
      "The matrix is said to be in reduced echelon form if in addition to being in the echelon form, the $p_i$'th columns have all but one zeroes. The non-zero entry is $a_{i p_i}$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####Example :\n",
      "The matrix\n",
      "\n",
      "0 1 4 5\n",
      "\n",
      "1 4 3 3\n",
      "\n",
      "0 0 0 5\n",
      "\n",
      "is not even in echelon form as $p_1 = 2$ and $p_2 = 1$ and hence $p_1 < p_2$ is not satisfied. Note, for the sake of giving an example,  here $p_3 = 4$.\n",
      "\n",
      "2 1 4 5\n",
      "\n",
      "1 4 3 3\n",
      "\n",
      "0 0 0 5\n",
      "\n",
      "is also not in echelon form because of the same reason.\n",
      "\n",
      "1 4 5 3\n",
      "\n",
      "0 3 8 2\n",
      "\n",
      "0 0 5 6\n",
      "\n",
      "does satisy the conditions on $p_i$s but $a_{2 p_2} \\neq 1$ and hence is not in echelon form.\n",
      "\n",
      "1 4 0 3\n",
      "\n",
      "0 1 0 2\n",
      "\n",
      "0 0 1 6\n",
      "\n",
      "is in echelon form but not in reduced echelon form as $p_2$-th column has extra non-zero entries (namely 4).\n",
      "\n",
      "1 0 0 3\n",
      "\n",
      "0 1 0 2\n",
      "\n",
      "0 0 1 6\n",
      "\n",
      "is in reduced echelon form.\n",
      "\n",
      "\n",
      "#### The algorithm to reduce to echelon form.\n",
      "1. First find the *first* non-zero column. Suppose the the column number is $p_1$ and the first non-zero entry is in the $i$-th row.\n",
      "1. Switch the 1st row with the $i$-th row.\n",
      "1. Sweep the $p_1$-th column\n",
      "1. Repeat this for the smaller matrix spanning row 2 - last row and column $p_1 + 1$ to the last column."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### First let us write a function to find the first column with a non-zero entry.\n",
      "It will return a tuple (row, column) for the non-zero entry. If the matrix is zero. It will return (-1, -1)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def first_non_zero(matrix) :\n",
      "    found_entry = False\n",
      "    column_scanning = 0\n",
      "    if matrix == [] or matrix == [[]] :\n",
      "        print \"The matrix is empty.\"\n",
      "        row = -1\n",
      "        col = -1\n",
      "    else :\n",
      "        while found_entry == False and column_scanning < len(matrix[0]) :\n",
      "            # len(matrix[0]) = number of columns of the matrix.\n",
      "            row_scanning = 0\n",
      "            while found_entry == False and row_scanning < len(matrix) :\n",
      "                #print \"r, c : \", row_scanning, column_scanning\n",
      "                if matrix[row_scanning][column_scanning] != 0 :\n",
      "                    found_entry = True\n",
      "                    row = row_scanning\n",
      "                    col = column_scanning\n",
      "                else :\n",
      "                    row_scanning += 1\n",
      "            column_scanning += 1\n",
      "        if found_entry == False :\n",
      "            row = -1\n",
      "            col = -1\n",
      "            print \"The matrix is the zero matrix.\"\n",
      "    return (row, col)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# To test it out :\n",
      "print first_non_zero([[0, 0, 0, 0], [0, 2, 1, 0], [0, 0, 5, 1], [0, 5, 5, 6]])\n",
      "print first_non_zero([[0, 0, 0, 0], [0, 0, 0, 0]])\n",
      "print first_non_zero([[], []])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(1, 1)\n",
        "The matrix is the zero matrix.\n",
        "(-1, -1)\n",
        "The matrix is the zero matrix.\n",
        "(-1, -1)\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we are ready to do Gaussian elimination"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gauss_elim_prelim(matrix) :\n",
      "    row = 0\n",
      "    row_offset = 0\n",
      "    col_offset = 0\n",
      "    \n",
      "    (r, c) = first_non_zero(matrix)\n",
      "    if r == -1 :\n",
      "        print \"Nothing to do: The matrix is zero.\"\n",
      "    else :\n",
      "        matrix = elem1(matrix, 0, r)\n",
      "        matrix = sweep(matrix, (0, c))\n",
      "        print_matrix(matrix, \"%6.2f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gauss_elim_prelim([[0,0,2,2],\\\n",
      "            [1,4,4,4],\\\n",
      "            [0,1,3,3]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  1.00   4.00   4.00   4.00\n",
        "  0.00   0.00   2.00   2.00\n",
        "  0.00   1.00   3.00   3.00\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gauss_elim(matrix) :\n",
      "    row = 0\n",
      "    row_offset = 0\n",
      "    col_offset = 0\n",
      "    \n",
      "    while row < len(matrix) :\n",
      "        submatrix = []\n",
      "        for i in range(len(matrix) - row_offset) :\n",
      "            row_submatrix = []\n",
      "            for j in range(len(matrix[0]) - col_offset) :\n",
      "                row_submatrix.append(matrix[i + row_offset][j + col_offset])\n",
      "            submatrix.append(row_submatrix)\n",
      "        #print_matrix(submatrix, \"%6.2f\")\n",
      "        print row, row_offset, col_offset, submatrix, matrix\n",
      "        print \"-\"*20\n",
      "                                     \n",
      "        (r, c) = first_non_zero(submatrix)\n",
      "        if r == -1 :\n",
      "            print \"Nothing to do: The (sub)matrix is zero.\"\n",
      "        else :\n",
      "            r += row_offset\n",
      "            c += col_offset\n",
      "            matrix = elem1(matrix, row_offset, r)\n",
      "            matrix = sweep(matrix, (row_offset, c))\n",
      "            row_offset += 1\n",
      "            col_offset = c+1\n",
      "        row += 1\n",
      "    return matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "read_matrix(\"files/B.txt\")\n",
      "print_matrix(gauss_elim(B), \"%6.2f\")\n",
      "print_matrix(gauss_elim([[0,0,0,0], [0,1,1,0], [0,0,1,0]]), \"%6.2f\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0 0 0 [[-0.6000000000000001, -0.40000000000000036, 0.0, -1.0], [0.2, 0.8, 1.0, 1.0], [2.6, -1.6, 0.0, 4.0]] [[-0.6000000000000001, -0.40000000000000036, 0.0, -1.0], [0.2, 0.8, 1.0, 1.0], [2.6, -1.6, 0.0, 4.0]]\n",
        "--------------------\n",
        "1 1 1 [[0.6666666666666666, 1.0, 0.6666666666666667], [-3.333333333333335, 0.0, -0.33333333333333304]] [[1.0, 0.6666666666666672, -0.0, 1.6666666666666665], [0.0, 0.6666666666666666, 1.0, 0.6666666666666667], [0.0, -3.333333333333335, 0.0, -0.33333333333333304]]\n",
        "--------------------\n",
        "2 2 2 [[5.000000000000002, 3.0000000000000018]] [[1.0, 0.0, -1.0000000000000009, 0.9999999999999993], [0.0, 1.0, 1.5, 1.0], [0.0, 0.0, 5.000000000000002, 3.0000000000000018]]\n",
        "--------------------\n",
        "  1.00   0.00   0.00   1.60\n",
        "  0.00   1.00   0.00   0.10\n",
        "  0.00   0.00   1.00   0.60\n",
        "0 0 0 [[0, 0, 0, 0], [0, 1, 1, 0], [0, 0, 1, 0]] [[0, 0, 0, 0], [0, 1, 1, 0], [0, 0, 1, 0]]\n",
        "--------------------\n",
        "1 1 2 [[0.0, 0.0], [1.0, 0.0]] [[0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]\n",
        "--------------------\n",
        "2 2 3 [[0.0]] [[0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0]]\n",
        "--------------------\n",
        "The matrix is the zero matrix.\n",
        "Nothing to do: The (sub)matrix is zero.\n",
        "  0.00   1.00   0.00   0.00\n",
        "  0.00   0.00   1.00   0.00\n",
        "  0.00   0.00   0.00   0.00\n"
       ]
      }
     ],
     "prompt_number": 19
    }
   ],
   "metadata": {}
  }
 ]
}