01.6.6 Programming without indices (axpy) - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Programming without indices (axpy)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Preliminaries"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<p><font color=red> Again, copy this notebook so that you don't corrupt the original!  Then you can \"play\" with the copy of the notebook all you want! </font> </p>\n",
      "\n",
      "<p>In this notebook, we show how the FLAME notation (the notation in which vectors and/or matrices are partitioned into regions) can be leveraged to implement linear algebra operations without using indices (which are the root of all evil in programming...).</p>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's start by importing numpy and creating vectors $ x = \\left( \\begin{array}{r} 1 \\\\ 2 \\\\ 3 \\end{array} \\right) $, $ y = \\left( \\begin{array}{r} -1 \\\\ 0 \\\\ -2 \\end{array} \\right) $ and a scalar $ \\alpha = 2.5 $.  \n",
      "\n",
      "Execute the code in the box by clicking in the box and then on \"Cell -> Run\".  Alternatively, click on the box and push \"Shift\" and \"Return\" (or \"Enter\") together."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np    # This imports a package called \"numpy\" that will make working with matrices \n",
      "                      # simpler\n",
      "\n",
      "import laff\n",
      "# create two two-dimensional matrices of only one column each.  \n",
      "x = np.matrix( '1.;2.;3.' )\n",
      "print( 'x = ' )\n",
      "print( x )\n",
      "\n",
      "y = np.matrix( '-1.;0.;-2.' )\n",
      "print( 'y = ' )\n",
      "print( y )\n",
      "\n",
      "alpha = 2.5\n",
      "print( 'alpha = ' )\n",
      "print( alpha )\n",
      "\n",
      "yold = np.matrix( np.zeros( [3,1] ) )\n",
      "laff.copy( y, yold )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Axpy as a simple routine"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<p>\n",
      "Here is a simple routine for computing $ y := \\alpha x + y $:\n",
      "</p>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def axpy( alpha, x, y ):\n",
      "\n",
      "    m, n = np.shape( x )\n",
      "    \n",
      "    for i in range( m ):\n",
      "        y[ i, 0 ] = alpha * x[ i, 0 ] + y[ i, 0 ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Be sure the run the above box, or this notebook won't know about the routine!!!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, execute"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "laff.copy( yold, y )\n",
      "\n",
      "print( 'y before axpy:')\n",
      "\n",
      "print( y )\n",
      "\n",
      "axpy( alpha, x, y )\n",
      "\n",
      "print( 'y after axpy: ' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare new y to  alpha * x + yold:' )\n",
      "print( y - ( alpha * x + yold ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "An implementation with the FLAMEPy Application Programming Interface (API)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now show how to implement this same routine using the FLAMEPy API.\n",
      "\n",
      "Start by visiting the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage </a>.  Follow along with the video and paste the resulting code below.  Then follow along with the video and add the appropriate commands.\n",
      "\n",
      "Here is the algorithm as presented in Unit 1.6.5.  \n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/FLAME_Algorithm_Axpy.png\" alt=\"some_text\">\n",
      "\n",
      "In the video for Unit 1.6.6, you are encouraged to us the Spark webpage and the above algorithm to implement the axpy operation.\n",
      "\n",
      "<b> Note: </b>\n",
      "<ul>\n",
      "<li> Make sure that the name of the routine is <code> axpy_unb </code>.\n",
      "<li> Make sure you order the parameters <code> alpha, x, y </code>.\n",
      "<li> Make sure you work through the vectors in the right direction.\n",
      "</ul>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff\n",
      "\n",
      "def axpy_unb(alpha, x, y):\n",
      "\n",
      "    xT, \\\n",
      "    xB  = flame.part_2x1(x, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    yT, \\\n",
      "    yB  = flame.part_2x1(y, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    while xT.shape[0] < x.shape[0]:\n",
      "\n",
      "        x0,   \\\n",
      "        chi1, \\\n",
      "        x2    = flame.repart_2x1_to_3x1(xT, \\\n",
      "                                        xB, \\\n",
      "                                        1, 'BOTTOM')\n",
      "\n",
      "        y0,   \\\n",
      "        psi1, \\\n",
      "        y2    = flame.repart_2x1_to_3x1(yT, \\\n",
      "                                        yB, \\\n",
      "                                        1, 'BOTTOM')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        psi1 = alpha * chi1 + psi1\n",
      "        \n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        xT, \\\n",
      "        xB  = flame.cont_with_3x1_to_2x1(x0,   \\\n",
      "                                         chi1, \\\n",
      "                                         x2,   \\\n",
      "                                         'TOP')\n",
      "\n",
      "        yT, \\\n",
      "        yB  = flame.cont_with_3x1_to_2x1(y0,   \\\n",
      "                                         psi1, \\\n",
      "                                         y2,   \\\n",
      "                                         'TOP')\n",
      "\n",
      "    flame.merge_2x1(yT, \\\n",
      "                    yB, y)\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import laff\n",
      "\n",
      "laff.copy( yold, y )\n",
      "\n",
      "alpha = np.matrix( '-2.0' ) # the way we are going to program, scalars, vectors, and matrices are \n",
      "                           # all just matrices.  So, alpha here is a 1 x 1 matrix, which we \n",
      "                           # initialize to some number, in this case -2.0.\n",
      "\n",
      "axpy_unb( alpha, x, y )    # Takes alpha, x, and y, and then updates y with the \n",
      "                          # result of axpy(alpha, x, y ).  Notice that the contents of variable y\n",
      "                          # are updated.  This only works if y is passed in as an array \n",
      "                          # (a matrix in our case)\n",
      "\n",
      "print( 'y' )\n",
      "print( y )\n",
      "\n",
      "print( 'compare updated y to  alpha * x + yold:' )\n",
      "\n",
      "print( y - ( alpha[ 0, 0 ] * x + yold )  )  # we have to use alpha[ 0, 0 ] because otherwise numpy \n",
      "                                            # thinks alpha is a matrix, and then the dimensions for \n",
      "                                            # matrix-matrix multiplication don't match up.\n",
      "                                            # More on this, later."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The output should be:\n",
      "<code>\n",
      "y\n",
      "[[-3.]\n",
      " [-4.]\n",
      " [-8.]]\n",
      "compare updated y to  alpha * x + yold:\n",
      "[[ 0.]\n",
      " [ 0.]\n",
      " [ 0.]]\n",
      "</code>"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}