diff options
Diffstat (limited to 'lecture-notes/advanced-python/arrays.rst')
| -rw-r--r-- | lecture-notes/advanced-python/arrays.rst | 736 |
1 files changed, 736 insertions, 0 deletions
diff --git a/lecture-notes/advanced-python/arrays.rst b/lecture-notes/advanced-python/arrays.rst new file mode 100644 index 0000000..d5699ce --- /dev/null +++ b/lecture-notes/advanced-python/arrays.rst @@ -0,0 +1,736 @@ +Arrays +====== + +In this section, we shall learn about a very powerful data structure, +specially useful for scientific computations, the array. We shall look at +initialization of arrays, operations on arrays and a brief comparison with +lists. + +Arrays are homogeneous data structures, unlike lists which can be +heterogeneous. They can have only one type of data as their entries. + +Arrays of a given length are comparatively much faster in mathematical +operations than lists of the same length, because of the fact that they are +homogeneous. + +Now let us see how to create arrays. + +To create an array we will use the function ``array()``, + +:: + + a1 = array([1,2,3,4]) + +Notice that we created a one dimensional array here. Also notice that we +passed a list as an argument, to create the array. + +We create two dimensional array by converting a list of lists to an array + +:: + + a2 = array([[1,2,3,4],[5,6,7,8]]) + +A two dimensional array has been created by passing a list of lists to the +array function. + +Now let us use ``arange()`` function to create the same arrays as before. + +:: + + ar1 = arange(1, 5) + +The ``arange`` function is similar to the range function that we have already +looked at, in the basic Python section. + +To create the two dimensional array, we first obtain a 1D array with the +elements from 1 to 8 and then convert it to a 2D array. + +To obtain the 1D array, we use the ``arange`` function + +:: + ar2 = arange(1, 9) + print ar2 + +We use the shape method to change the shape of the array, into a 2D array. + +:: + + ar2.shape = 2, 4 + print ar2 + +This gives us the required array. + +Instead of passing a list to the array function, we could also pass a list +variable. For, instance + +:: + + l1 = [1,2,3,4] + a3 = array(l1) + +We used the shape method to change the shape of the array. But it can also +be used to find out the shape of an array that we have. + +:: + a1.shape + a2.shape + +As, already stated, arrays are homogeneous. Let us try to create a new +array with a mix of elements and see what happens. + +:: + + a4 = array([1,2,3,'a string']) + +We expect an error, but there wasn't one. What happened? Let's see what a4 +has. + +:: + + a4 + +There you go! The data type of the array has been set to a string of +length 8. All the elements have been implicitly type-casted as strings, +though our first three elements were meant to be integers. The dtype +specifies the data-type of the array. + +There some other special methods as well, to create special types of +arrays. We can create an 2 dimensional identity array, using the function +``identity()``. The function ``identity()`` takes an integer argument which +specifies the size of the diagonal of the array. + +:: + + identity(3) + +As you can see the identity function returned a three by three array, with +all the diagonal elements as ones and the rest of the elements as zeros. + +``zeros()`` function accepts a tuple, which is the shape of the array that +we want to create, and it generates an array with all elements as zeros. + +Let us creates an array of the order four by five with all the elements +zero. We can do it using the method zeros, :: + + zeros((4,5)) + +Notice that we passed a tuple to the function zeros. Similarly, ``ones()`` +gives us an array with all elements as ones. + +Also, ``zeros_like`` gives us an array with all elements as zeros, with a +shape similar to the array passed as the argument and the same data-type as +the argument array. + +:: + + a = zeros_like([1.5, 1, 2, 3] + print a, a.dtype + +Similarly, the function ``ones_like``. + +Let us try out some basic operations with arrays, before we end this +section. + +Let's check the value of a1, by typing it out on the interpreter. + +:: + + a1 + +``a1`` is a single dimensional array, array([1,2,3,4]). Now, try + +:: + + a1 * 2 + +We get back a new array, with all the elements multiplied by 2. Note that +the value of elements of ``a1`` remains the same. + +:: + + a1 + +Similarly we can perform addition or subtraction or division. + +:: + + a1 + 3 + a1 - 7 + a1 / 2.0 + +We can change the array, by simply assigning the newly returned array to +the old array. + +:: + + a1 = a1 + 2 + +Notice the change in elements of a1, + +:: + + a1 + +We could also do the augmented assignment, but there's a small nuance here. +For instance, + +:: + + a = arange(1, 5) + b = arange(1, 5) + + print a, a.dtype + print b, b.dtype + + a = a/2.0 + b /= 2.0 + + print a, a.dtype + print b, b.dtype + +As you can see, ``b`` doesn't have the expected result because the +augmented assignment that we are doing, is an inplace operation. Given that +arrays are optimized to be very fast, by fixing their datatype and hence +the amount of memory they can occupy, the division by a float, when +performed inplace, fails to change the dtype of the array. So, be cautious, +when using in-place assignments. + +Now let us try operations between two arrays. + +:: + + a1 + a1 + a1 * a2 + +Note that both the addition and the multiplication are element-wise. + +Accessing pieces of arrays +========================== + +Now, that we know how to create arrays, let's see how to access individual +elements of arrays, get rows and columns and other chunks of arrays using +slicing and striding. + +Let us have two arrays, A and C, as the sample arrays that we will use to +learn how to access parts of arrays. + +:: + + A = array([12, 23, 34, 45, 56]) + + C = array([[11, 12, 13, 14, 15], + [21, 22, 23, 24, 25], + [31, 32, 33, 34, 35], + [41, 42, 43, 44, 45], + [51, 52, 53, 54, 55]]) + +Let us begin with the most elementary thing, accessing individual elements. +Also, let us first do it with the one-dimensional array A, and then do the +same thing with the two-dimensional array. + +To access, the element 34 in A, we say, + +:: + + A[2] + +A of 2, note that we are using square brackets. + +Like lists, indexing starts from 0 in arrays, too. So, 34, the third +element has the index 2. + +Now, let us access the element 34 from C. To do this, we say + +:: + + C[2, 3] + +C of 2,3. + +34 is in the third row and the fourth column, and since indexing begins +from zero, the row index is 2 and column index is 3. + +Now, that we have accessed one element of the array, let us change it. We +shall change the 34 to -34 in both A and C. To do this, we simply assign +the new value after accessing the element. + +:: + + A[2] = -34 + C[2, 3] = -34 + +Now that we have accessed and changed a single element, let us access and +change more than one element at a time; first rows and then columns. + +Let us access one row of C, say the third row. We do it by saying, + +:: + + C[2] + +How do we access the last row of C? We could say, + +:: + + C[4] + +for the fifth row, or as with lists, use negative indexing and say + +:: + + C[-1] + +Now, we could change the last row into all zeros, using either + +:: + + C[-1] = [0, 0, 0, 0, 0] + +or + +:: + + C[-1] = 0 + +Now, how do we access one column of C? As with accessing individual +elements, the column is the second parameter to be specified (after the +comma). The first parameter, is replaced with a ``:``. This specifies that +we want all the elements of that dimension, instead of just one particular +element. We access the third column by + +:: + + C[:, 2] + +So, we can change the last column of C to zeroes, by + +:: + + C[:, -1] = 0 + +Since A is one dimensional, rows and columns of A don't make much sense. It +has just one row and + +:: + + A[:] + +gives the whole of A. + +Now, that we know how to access, rows and columns of an array, we shall +learn how to access other, larger pieces of an array. For this purpose, we +will be using image arrays. + +To read an image into an array, we use the ``imread`` command. We shall use +the image ``squares.png``. We shall first navigate to that path in the OS +and see what the image contains. + +Let us now read the data in ``squares.png`` into the array ``I``. + +:: + + I = imread('squares.png') + +We can see the contents of the image, using the command ``imshow``. We say, + +:: + + imshow(I) + +to see what has been read into ``I``. We do not see white and black +because, ``pylab`` has mapped white and black to different colors. This can +be changed by using a different colormap. + +To see that ``I`` is really, just an array, we say, + +:: + + I + +at the prompt, and see that an array is displayed. + +To check the dimensions of any array, we can use ``.shape``. We say + +:: + + I.shape + +to get the dimensions of the image. As we can see, ``squares.png`` has the +dimensions of 300x300. + +Our goal for now, is be to get the top-left quadrant of the image. To do +this, we need to access, a few of the rows and a few of the columns of the +array. + +To access, the third column of C, we said, ``C[:, 2]``. Essentially, we are +accessing all the rows in column three of C. Now, let us modify this to +access only the first three rows, of column three of C. +We say, + +:: + + C[0:3, 2] + +to get the elements of rows indexed from 0 to 3, 3 not included and column +indexed 2. Note that, the index before the colon is included and the index +after it is not included in the slice that we have obtained. This is very +similar to the ``range`` function, where ``range`` returns a list, in which +the upper limit or stop value is not included. + +Now, if we wish to access the elements of row with index 2, and in +columns indexed 0 to 2 (included), we say, + +:: + + C[2, 0:3] + +Note that when specifying ranges, if you are starting from the beginning or +going up-to the end, the corresponding element may be dropped. + +:: + + C[2, :3] + +also gives us the same elements, [31, 32, 33] + +Now, we are ready to obtain the top left quarter of the image. How do we go +about doing it? Since, we know the shape of the image to be 300, we know +that we need to get the first 150 rows and first 150 columns. + +:: + + I[:150, :150] + +gives us the top-left corner of the image. + +We use the ``imshow`` command to see the slice we obtained in the +form of an image and confirm. + +:: + + imshow(I[:150, :150]) + +How would we obtain the square in the center of the image? + +:: + + imshow(I[75:225, 75:225]) + +Our next goal is to compress the image, using a very simple technique to +reduce the space that the image takes on disk while not compromising too +heavily on the image quality. The idea is to drop alternate rows and +columns of the image and save it. This way we will be reducing the data to +a fourth of the original data but losing only so much of visual +information. + +We shall first learn the idea of striding using the smaller array C. +Suppose we wish to access only the odd rows and columns (first, third, +fifth). We do this by, + +:: + + C[0:5:2, 0:5:2] + +if we wish to be explicit, or simply, +:: + + C[::2, ::2] + +This is very similar to the step specified to the ``range`` function. It +specifies, the jump or step in which to move, while accessing the elements. +If no step is specified, a default value of 1 is assumed. + +:: + + C[1::2, ::2] + +gives the elements, [[21, 23, 0], [41, 43, 0]] + +Now, that we know how to stride over an array, we can drop alternate rows +and columns out of the image in I. + +:: + + I[::2, ::2] + +To see this image, we say, +:: + + imshow(I[::2, ::2]) + +This does not have much data to notice any real difference, but notice that +the scale has reduced to show that we have dropped alternate rows and +columns. If you notice carefully, you will be able to observe some blurring +near the edges. To notice this effect more clearly, increase the step to 4. + +:: + + imshow(I[::4, ::4]) + +That brings us to our discussion on accessing pieces of arrays. We shall +look at matrices, next. + +Matrices +======== + +In this course, we shall perform all matrix operations using the array +data-structure. We shall create a 2D array and treat it as a matrix. + +:: + + m1 = array([[1,2,3,4]]) + m1.shape + +As we already know, we can get a 2D array from a list of lists as well. + +:: + + l1 = [[1,2,3,4],[5,6,7,8]] + m2 = array(l1) + m3 = array([[5,6,7,8],[9,10,11,12]]) + +We can do matrix addition and subtraction as, + +:: + + m3 + m2 + +does element by element addition, thus matrix addition. + +Similarly, + +:: + + m3 - m2 + +it does matrix subtraction, that is element by element +subtraction. Now let us try, + +:: + + m3 * m2 + +Note that in arrays ``m3 * m2`` does element wise multiplication and not +matrix multiplication, + +And matrix multiplication in matrices are done using the function ``dot()`` + +:: + + dot(m3, m2) + +but for matrix multiplication, we need arrays of compatible sizes and hence +the multiplication fails, in this case. + +To see an example of matrix multiplication, we choose a proper pair. + +:: + + m1.shape + +m1 is of the shape one by four, let us create an array of the shape four by +two, + +:: + + m4 = array([[1,2],[3,4],[5,6],[7,8]]) + dot(m1, m4) + +Thus, the function ``dot()`` can be used for matrix multiplication. + +We have already seen the special functions like ``identity()``, +``zeros()``, ``ones()``, etc. to create special arrays. + +Let us now look at some Matrix specific operations. + +To find out the transpose, we use the ``.T`` method. + +:: + + print m4 + print m4.T + +Now let us try to find out the Frobenius norm of inverse of a 4 by 4 +matrix, the matrix being, + +:: + + m5 = arange(1,17).reshape(4,4) + print m5 + +The inverse of a matrix A, A raise to minus one is also called the +reciprocal matrix such that A multiplied by A inverse will give 1. + +:: + + im5 = inv(m5) + +The Frobenius norm of a matrix is defined as square root of sum of squares +of elements in the matrix. The Frobenius norm of ``im5`` can be found by, + +:: + + sqrt(sum(im5 * im5)) + +Now try to find out the infinity norm of the matrix im5. The infinity norm +is defined as the maximum value of sum of the absolute of elements in each +row. + +The solution for the problem is, + +:: + + max(sum(abs(im5), axis=1)) + +Well! to find the Frobenius norm and Infinity norm we have an even easier +method, and let us see that now. + +The norm of a matrix can be found out using the method ``norm()``. In order +to find out the Frobenius norm of the im5, we do, + +:: + + norm(im5) + +And to find out the Infinity norm of the matrix im5, we do, +:: + + norm(im5,ord=inf) + + +The determinant of a square matrix can be obtained using the function +``det()`` and the determinant of m5 by, + +:: + + det(m5) + +The eigen values and eigen vector of a square matrix can be computed +using the function ``eig()`` and ``eigvals()``. + +Let us find out the eigen values and eigen vectors of the matrix +m5. We can do it as, + +:: + + eig(m5) + + +Note that it returned a tuple of two matrices. The first element in +the tuple are the eigen values and the second element in the tuple are +the eigen vectors. Thus the eigen values are, +:: + + eig(m5)[0] + +and the eigen vectors are, + +:: + + eig(m5)[1] + +The eigen values can also be computed using the function ``eigvals()`` as, + +:: + + eigvals(m5) + + +We can also find the singular value decomposition or S V D of a matrix. + +The SVD of ``m5`` can be found by + +:: + + svd(m5) + +Notice that it returned a tuple of 3 elements. The first one U the +next one Sigma and the third one V star. + +That brings us to our discussion of matrices and operations on them. But we +shall continue to use them in the next section on Least square fit. + +Least square fit +================ + +First let us have a look at the problem. + +We have an input file generated from a simple pendulum experiment, which we +have already looked at. We shall find the least square fit of the plot +obtained by plotting l vs. t^2, where l is the length of the pendulum and t +is the corresponding time-period. + +:: + + l, t = loadtxt("/home/fossee/pendulum.txt", unpack=True) + l + t + +We can see that l and t are two sequences containing length and time values +correspondingly. + +Let us first plot l vs t^2. Type +:: + + tsq = t * t + plot(l, tsq, 'bo') + +We can see that there is a visible linear trend, but we do not get a +straight line connecting them. We shall, therefore, generate a least square +fit line. + +We are first going to generate the two matrices ``tsq`` and A, the vander +monde matrix. Then we are going to use the ``lstsq`` function to find the +values of m and c. + +Let us now generate the A matrix with l values. We shall first generate a 2 +x 90 matrix with the first row as l values and the second row as ones. Then +take the transpose of it. Type + +:: + + A = array((l, ones_like(l))) + A + +We see that we have intermediate matrix. Now we need the transpose. Type + +:: + + A = A.T + A + +Now we have both the matrices A and tsq. We only need to use the ``lstsq`` + +:: + + result = lstsq(A, tsq) + +The result is a sequence of values. The first item in this sequence, is the +matrix p i.e., the values of m and c. Hence, + +:: + + m, c = result[0] + m + c + +Now that we have m and c, we need to generate the fitted values of t^2. Type + +:: + + tsq_fit = m * l + c + plot(l, tsq, 'bo') + plot(l, tsq_fit, 'r') + +We get the least square fit of l vs t^2 + +That brings us to the end of our discussion on least square fit curve and +also our discussion of arrays. + +.. + Local Variables: + mode: rst + indent-tabs-mode: nil + sentence-end-double-space: nil + fill-column: 75 + End: |