Advanced Python

Interactive Plotting

We have looked at the basics of Python. In this section we shall look at the basics of plotting.

Lets start our IPython interpreter, by

$ ipython -pylab

Pylab is a python library which provides plotting functionality. It also provides many other important mathematical and scientific functions. After running ipython -pylab in your shell if at the top of the output, you see something like

`ERROR: matplotlib could NOT be imported! Starting normal IPython.`

Then you have to install python-matplotlib and run this command again.

Let us get started directly and try plotting a cosine curve in the range -pi to pi.

p = linspace(-pi,pi,100)
plot(p, cos(p))

As you can see we have obtained a cosine curve in the required range. Let’s try and understand what we have done here.

The variable p has a hundred points in the range -pi to pi. To see this, check the documentation of the linspace command.

linspace?

As you can see, it returns num evenly spaced samples, calculated over the interval start and stop.

You can verify this, by looking at the first and last points of p and the length of p

print p[0], p[-1], len(p)

Now, let’s clear the plot we obtained and plot a sine curve in the same range.

clf()
plot(p, sin(p))

clf command clears the plot, and the plot command following it, plots the required curve.

Let’s briefly explore the plot window and look at what features it has.

Firstly, note that moving the mouse around gives us the point where mouse points at.

Also we have some buttons the right most among them is for saving the file.

Left to the save button is the slider button to specify the margins.

Left to this is zoom button to zoom into the plot. Just specify the region to zoom into.

The button left to it can be used to move the axes of the plot.

The next two buttons with a left and right arrow icons change the state of the plot and take it to the previous state it was in. It more or less acts like a back and forward button in the browser.

The last one is ‘home’, which returns the plot to the original view.

Embellishing Plots

Now that we know how to make simple plots, let us look at embellishing the plots. We shall look at how to modify the colour, thickness and line-style of a plot. We shall then learn how to add title to a plot and then look at adding labels to x and y axes. We shall also look at adding annotations to the plot and setting the limits on the axes.

Let us decorate the sine curve that we have already plotted. If you have closed that window or your IPython terminal, redo the plot.

p = linspace(-pi,pi,100)
plot(p, sin(p))

As we can see, the default colour and the default thickness of the line is as decided by pylab. Wouldn’t it be nice if we could control these parameters in the plot? This is possible by passing additional arguments to the plot command.

The additional argument that we shall be passing in here now is the color argument. We shall first clear the figure and plot the same in red color.

clf()
plot(p, sin(p), 'r')

As we can see we have the same plot but now in red color.

To alter the thickness of the line, we use the linewidth argument in the plot command.

plot(p, cos(p), linewidth=2)

produces a plot with a thicker line, to be more precise plot with line thickness 2.

Occasionally we would also want to alter the style of line. Sometimes all we want is just a bunch of points not joined. This is possible by passing the linestyle argument along with or instead of the colour argument.

clf()
plot(p, sin(p), '.')

produces a plot with only points.

To produce the same plot but now in blue colour, we do

clf()
plot(p, sin(p), 'b.')

Other available options can be seen in the documentation of plot.

plot?

Now that we know how to produce a bare minimum plot with colour, style and thickness of our interest, we shall look at decorating the plot.

Let us start with a plot of the function -x^2 + 4x - 5 in the range -2 to 4.

x = linspace(-2, 4, 50)
plot(x, -x*x + 4*x - 5, 'r', linewidth=2)

We now have the plot in a colour and linewidth of our interest. As you can see, the figure does not have any description describing the plot.

We will now add a title to the plot by using the title command.

title("Parabolic function -x^2+4x-5")

The figure now has a title which describes what the plot is. The title command as you can see, takes a string as an argument and sets the title accordingly.

The formatting in title is messed and it does not look clean. You can imagine what would be the situation if there were fractions and more complex functions like log and exp. Wouldn’t it be good if there was LaTeX like formatting?

That is also possible by adding a $ sign before and after the part of the string that should be in LaTeX style. For instance,

title("Parabolic function $-x^2+4x-5$")

gives us the polynomial formatted properly.

Although we have title, the plot is not complete without labelling x and y axes. Hence we shall label x-axis to “x” and y-axis to “f(x)”

xlabel("x")

As you can see, xlabel command takes a string as an argument, similar to the title command and sets it as the label to x-axis.

Similarly,

ylabel("f(x)")

sets the name of the y-axis as “f(x)”

Let us now see how to name or annotate a point on the plot. For example the point (2, -1) is the local maxima. We would like to name the point accordingly. We can do this by using

annotate("local maxima", xy=(2, -1))

As you can see, the first argument to annotate command is the name we would like to mark the point as and the second argument is the co-ordinates of the point at which the name should appear. It is a sequence containing two numbers. The first is x co-ordinate and second is y co-ordinate.

As we can see, every annotate command makes a new annotation on the figure.

Now we have everything we need to decorate a plot. But, it would be nice to change the limits of the plot area, so that it looks better.

We shall look at how to get and set them from the script. We shall first get the existing limits and then change them as required.

xlim()
ylim()

We see that xlim function returns the current x axis limits and ylim function returns the current y-axis limits.

Let us look at how to set the limits.

xlim(-4, 5)

We see the limits of x-axis are now set to -4 and 5.

Similarly

ylim(-15, 2)

sets the limits of y-axis appropriately.

Saving to Scripts

While we are at it, let’s look at how to save all the things we typed into our interpreter into scripts that can be used anytime we like.

Let’s now look at the history of the commands we have typed. The history can be retreived by using %hist command.

%hist

As you can see, it displays a list of recent commands that we typed. Every command has a number in front, to specify in which order and when it was typed.

Please note that there is a % sign before the hist command. This implies that %hist is a command that is specific to IPython and not available in the vanilla Python interpreter. These type of commands are called as magic commands.

Also note that, the %hist itself is a command and is displayed as the most recent command. We should note that anything we type in is stored as history, irrespective of whether it is a valid command or an error or an IPython magic command.

If we want only the recent 5 commands to be displayed, we pass the number as an argument to %hist command.

%hist 5

displays the recent 5 commands, inclusive of the %hist command. The default number is 40.

To list all the commands between 5 and 10, type

%hist 5-10

Now that we have the history, we would like to save the required lines of code from history to reproduce the plot of the parabolic function. This is possible by using the %save command.

Before we do that, let us first look at history and identify what lines of code we require.

%hist

The first command is linspace. But second command is a command that gave us an error. Hence we do not need second command. The commands from third to sixth and eighth are required.

%save plot_script.py 1 3-6 8

The command saves first and then third to sixth and eighth lines of code into the specified file.

The first argument to %save is the path of file to save the commands and the arguments there after are the commands to be saved in the given order.

Now that we have the required lines of code in a file, let us learn how to run the file as a python script. We use the IPython magic command %run to do this.

%run -i plot_script.py

The script runs but we do not see the plot. This happens because when we are running a script and we are not in interactive mode anymore.

Hence on your IPython terminal type

show()

to show the plot.

The reason for including a -i after run is to tell the interpreter that if any name is not found in script, search for it in the interpreter. Hence all these sin, plot, pi and show which are not available in script, are taken from the interpreter and used to run the script.

Saving Plots

Let us now learn to save the plots, from the command line, in different formats.

Let us plot a sine curve from minus 3pi to 3pi.

x = linspace(-3*pi,3*pi,100)
plot(x,sin(x))

As, you can see we now have a sine curve. Let’s now see how to save the plot.

For saving the plot, we will use savefig() function, and it has to be done with the plot window open. The statement is,

savefig('sine.png')

Notice that savefig function takes one argument which is the filename. The last 3 characters after the . in the filename is the extension or type of the file which determines the format in which you want to save.

Also, note that we gave the full path or the absolute path to which we want to save the file.

Here we have used an extension .png which means the plot gets saved as a PNG image.

You can check file which has been saved as sine.png

savefig can save the plot in many formats, such as pdf, ps, eps, svg and png.

Multiple Plots

Let us now learn, how to draw more than one plot, how to add legends to each plot to indicate what each plot represents. We will also learn how to switch between the plots and create multiple plots with different regular axes which are also called as subplots.

Let us first create set of points for our plot. For this we will use the command called linspace:

x = linspace(0, 50, 10)

As we know, x has 10 points in the interval between 0 and 50 both inclusive.

Now let us draw a plot simple sine plot using these points

plot(x, sin(x))

This should give us a nice sine plot.

Oh! wait! It isn’t as nice, as we expected. The curve isn’t very smooth, why? In the linspace command, we chose too few points in a large interval between 0 and 50 for the curve to be smooth. So now let us use linspace again to get 500 points between 0 and 50 and draw the sine plot

y = linspace(0, 50, 500)
plot(y, sin(y))

Now we see a smooth curve. We can also see that, we have two plots, one overlaid upon another. Pylab overlays plots by default.

Since, we have two plots now overlaid upon each other we would like to have a way to indicate what each plot represents to distinguish between them. This is accomplished using legends. Equivalently, the legend command does this for us

legend(['sine-10 points', 'sine-500 points'])

Now we can see the legends being displayed for the respective plots. The legend command takes a single list of parameters where each parameter is the text indicating the plots in the order of their serial number.

Additionally, we could give the legend command an additional argument to choose the location of the placement, manually. By default, pylab places it in the location it thinks is the best. The example below, places it in the center of the plot. Look at the doc-string of legend for other locations.

legend(['sine-10 points', 'sine-500 points'], loc='center')

We now know how to draw multiple plots and use legends to indicate which plot represents what function, but we would like to have more control over the plots we draw. Like plot them in different windows, switch between them, perform some operations or labeling on them individually and so on. Let us see how to accomplish this. Before we move on, let us clear our screen.

clf()

To accomplishing this, we use the figure command

x = linspace(0, 50, 500)
figure(1)
plot(x, sin(x), 'b')
figure(2)
plot(x, cos(x), 'g')

Now we have two plots, a sine plot and a cosine plot in two different figures.

The figure command takes an integer as an argument which is the serial number of the plot. This selects the corresponding plot. All the plot commands we run after this are applied to the selected plot. In this example figure 1 is the sine plot and figure 2 is the cosine plot. We can, for example, save each plot separately

savefig('cosine.png')
figure(1)
title('sin(y)')
savefig('sine.png')

We also titled our first plot as ‘sin(y)’ which we did not do for the second plot.

Let us now do the following. Draw a line of the form y = x as one figure and another line of the form y = 2x + 3. Switch back to the first figure, annotate the x and y intercepts. Now switch to the second figure and annotate its x and y intercepts. Save each of them.

To solve this problem we should first create the first figure using the figure command. Before that, let us first run clf command to make sure all the previous plots are cleared

clf()
figure(1)
x = linspace(-5, 5, 100)
plot(x, x)

Now we can use figure command to create second plotting area and plot the figure

figure(2)
plot(x, ((2 * x) + 3))

Now to switch between the figures we can use figure command. So let us switch to figure 1. We are asked to annotate x and y intercepts of the figure 1 but since figure 1 passes through origin we will have to annotate the origin. We will annotate the intercepts for the second figure and save them as follows

figure(1)
annotate('Origin', xy=(0.0, 0.0)
figure(2)
annotate('x-intercept', xy=(0, 3))
annotate('y-intercept', xy=(0, -1.5))
savefig('plot2.png')
figure(1)
savefig('plot1.png')

We can close the figures from the terminal by using the close() command.

close()
close()

The first close command closes figure 1, and the second one closes the figure 2. We could have also use the close command with an argument ‘all’, to close all the figures.

close('all')

At times we run into situations where we want to compare two plots and in such cases we want to draw both the plots in the same plotting area. The situation is such that the two plots have different regular axes which means we cannot draw overlaid plots. In such cases we can draw subplots.

We use subplot command to accomplish this

subplot(2, 1, 1)

subplot command takes three arguments, number of rows, number of columns and the plot number, which specifies what subplot must be created now in the order of the serial number. In this case we passed 1 as the argument, so the first subplot that is top half is created. If we execute the subplot command as

subplot(2, 1, 2)

the lower subplot is created. Now we can draw plots in each of the subplot area using the plot command.

x = linspace(0, 50, 500)
plot(x, cos(x))

subplot(2, 1, 1)
y = linspace(0, 5, 100)
plot(y, y ** 2)

This created two plots one in each of the subplot area. The top subplot holds a parabola and the bottom subplot holds a cosine curve.

As seen here we can use subplot command to switch between the subplot as well, but we have to use the same arguments as we used to create that subplot, otherwise the previous subplot at that place will be automatically erased. It is clear from the two subplots that both have different regular axes. For the cosine plot x-axis varies from 0 to 100 and y-axis varies from 0 to 1 where as for the parabolic plot the x-axis varies from 0 to 10 and y-axis varies from 0 to 100

Plotting Data

We often require to plot points obtained from experimental observations, instead of the analytic functions that we have been plotting, until now. We shall now learn to read data from files and read it into sequences that can later be used to plot.

We shall use the loadtxt command to load data from files. We will be looking at how to read a file with multiple columns of data and load each column of data into a sequence.

Now, Let us begin with reading the file primes.txt, which contains just a list of primes listed in a column, using the loadtxt command. The file, in our case, is present in primes.txt.

Now let us read this list into the variable primes.

primes = loadtxt('primes.txt')

primes is now a sequence of primes, that was listed in the file, primes.txt.

We now type, print primes to see the sequence printed.

We observe that all of the numbers end with a period. This is so, because these numbers are actually read as floats.

Now, let us use the loadtxt command to read a file that contains two columns of data, pendulum.txt. This file contains the length of the pendulum in the first column and the corresponding time period in the second. Note that loadtxt needs both the columns to have equal number of rows.

Let us, now, read the data into the variable pend. Again, it is assumed that the file is in our current working directory.

pend = loadtxt('pendulum.txt')

Let us now print the variable pend and see what’s in it.

print pend

Notice that pend is not a simple sequence like primes. It has a sequence of items in which each item contains two values. Let us use an additional argument of the loadtxt command, to read it into two separate, simple sequences. We add the argument unpack=True to the loadtxt command.

L, T = loadtxt('pendulum.txt', unpack=True)

Let us now, print the variables L and T, to see what they contain.

print L
print T

Notice, that L and T now contain the first and second columns of data from the data file, pendulum.txt, and they are both simple sequences. unpack=True has given us the two columns into two separate sequences instead of one complex sequence.

Now that we have the required data in sequences, let us see how to plot it.

Since we already have the values of L and T as two different sequences, we now need to calculate T squared. We shall be plotting L vs. T^2 values.

To obtain the square of sequence T we will use the function square

Tsq = square(T)

Now to plot L vs T^2 we will simply type

plot(L, Tsq, '.')

For any experimental data, there is always an error in measurements due to instrumental and human constraints. Now we shall try and take into account error into our plots.

We shall read the read the error values from the file pendulum_error.txt along with the L and T values.

L, T, L_err, T_err = loadtxt('pendulum_error.txt', unpack=True)

Now to plot L vs T^2 with an error bar we use the function errorbar()

errorbar(L, Tsq , xerr=L_err, yerr=T_err, fmt='b.')

This gives a plot with error bar for x and y axis. The dots are of blue color. The parameters xerr and yerr are error on x and y axis and fmt is the format of the plot.

You can explore other options of errorbar by looking at it’s documentation.

errorbar?

Other kinds of Plots

We shall now briefly look at a few other kinds of plots, namely, the scatter plot, the pie chart, the bar chart and the log-log plot.

Let us start with scatter plot.

In a scatter plot, the data is displayed as a collection of points, each having the value of one variable determining the position on the horizontal axis and the value of the other variable determining the position on the vertical axis. This kind of plot is also called a scatter chart, a scatter diagram or a scatter graph.

Now, let us plot a scatter plot showing the percentage profit of a company A from the year 2000-2010. The data for the same is available in the file company-a-data.txt.

The data file has two lines with a set of values in each line, the first line representing years and the second line representing the profit percentages.

To produce the scatter plot, we first need to load the data from the file using loadtxt.

year,profit = loadtxt('company-a-data.txt', dtype=int)

By default loadtxt converts the value to float. The dtype=type(int()) argument in loadtxt converts the value to integer as we require the data as integers.

In-order to generate the scatter graph we will use the function scatter()

scatter(year, profit)

Notice that we passed two arguments to scatter() function, first one the values in x-coordinate, year, and the other the values in y-coordinate, the profit percentage.

Now let plot a pie chart for the same data. A pie chart or a circle graph is a circular chart divided into sectors, illustrating proportion.

Plot a pie chart representing the profit percentage of company A, with the same data from file company-a-data.txt. We shall reuse the data we have already read from the file.

We can plot the pie chart using the function pie().

pie(profit, labels=year)

Notice that we passed two arguments to the function pie(). First one the values and the next one the set of labels to be used in the pie chart.

Now let us move on to the bar charts. A bar chart or bar graph is a chart with rectangular bars with lengths proportional to the values that they represent.

Plot a bar chart representing the profit percentage of company A, with the same data from file company-a-data.txt.

We can plot the bar chart using the function bar().

bar(year, profit)

Note that the function bar() needs at least two arguments one the values in x-coordinate and the other values in y-coordinate which is used to determine the height of the bars.

Now let us move on to the log-log plot. A log-log graph or a log-log plot is a two-dimensional graph of numerical data that uses logarithmic scales on both the horizontal and vertical axes.

Due to the nonlinear scaling of the axes, a function of the form y = ax^b will appear as a straight line on a log-log graph.

Plot a log-log chart of y=5*x^3 for x from 1-20.

Before we actually plot let us calculate the points needed for that.

x = linspace(1,20,100)
y = 5*x**3

Now we can plot the log-log chart using loglog() function,

loglog(x, y)

For the sake of clarity, let us make a linear plot of x vs. y.

plot(x, y)

Observing the two plots together should give you some clarity about the loglog plot.

Help about matplotlib can be obtained from http://matplotlib.sourceforge.net/contents.html

That brings us to the end of the discussion on plots and matplotlib. We have looked a making simple plots, embellishing plots, saving plots, making multiple plots, plotting data from files, and a few varieties of plots.

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.

In this section we shall look at using scipy to do various common scientific computation tasks. We shall be looking at solving equations (linear and non-linear) and solving ODEs. We shall also briefly look at FFTs.

Solving Equations

Let us begin with solving equations, specifically linear equations.

Consider the set of equations,

3x + 2y -z = 1,
2x-2y + 4z = -2,
-x+ 1/2 y-z = 0

We shall use the solve function, to solve the given system of linear equations. solve requires the coefficients and the constants to be in the form of matrices, of the form Ax = b to solve the system.

We begin by entering the coefficients and the constants as matrices.

A = array([[3,2,-1],
           [2,-2,4],
           [-1, 0.5, -1]])

A is a 3X3 matrix of the coefficients of x, y and z

b = array([1, -2, 0])

Now, we can use the solve function to solve the given system.

x = solve(A, b)

Type x, to look at the solution obtained.

The equation is of the form Ax = b, so we verify the solution by obtaining a matrix product of A and x, and comparing it with b. As we have seen earlier, we should use the dot function, for a matrix product and not the * operator.

Ax = dot(A, x)
Ax

The result Ax, doesn’t look exactly like b, but if we carefully observe, we will see that it is the same as b. To save ourself all this trouble, we can use the allclose function.

allclose checks if two matrices are close enough to each other (with-in the specified tolerance level). Here we shall use the default tolerance level of the function.

::
allclose(Ax, b)

The function returns True, which implies that the product of A & x is very close to the value of b. This validates our solution.

Let’s move to finding the roots of a polynomial. We shall use the roots function for this.

The function requires an array of the coefficients of the polynomial in the descending order of powers. Consider the polynomial x^2-5x+6 = 0

coeffs = [1, -5, 6]
roots(coeffs)

As we can see, roots returns the result in an array. It even works for polynomials with imaginary roots.

roots([1, 1, 1])

As you can see, the roots of that equation are complex.

What if, we want the solution of non linear equations?

For that we shall use the fsolve function. We shall use the equation sin(x)+cos^2(x), as an example. fsolve is not part of the pylab package which we imported at the beginning, so we will have to import it. It is part of scipy package. Let’s import it,

from scipy.optimize import fsolve

We shall look at the details of importing, later.

Now, let’s look at the documentation of fsolve by typing fsolve?

fsolve?

As mentioned in documentation the first argument, func, is a python function that takes atleast one argument. The second argument, x0, is the initial estimate of the roots of the function. Based on this initial guess, fsolve returns a root.

Let’s define a function called f, that returns the value of (sin(x)+cos^2(x)) evaluated at the input value x.

def f(x):
    return sin(x)+cos(x)*cos(x)

We can test our function, by calling it with an argument for which the output value is known, say x = 0. We can see that

Let’s check our function definition, by calling it with 0 as an argument.

f(0)

We can see that the output is 1 as expected, since sin(x) + cos^2(x) has a value of 1, when x = 0.

Now, that we have our function, we can use fsolve to obtain a root of the expression sin(x)+cos^2(x). Let’s use 0 as our initial guess.

fsolve(f, 0)

fsolve has returned a root of sin(x)+cos^2(x) that is close to 0.

That brings us to the end of our discussion on solving equations. We discussed solution of linear equations, finding roots of polynomials and non-linear equations.

ODEs

Let’s see how to solve Ordinary Differential Equations (ODEs), using Python. Let’s consider the classic problem of the spread of an epidemic in a population, as an example.

This is given by the ordinary differential equation dy/dt = ky(L-y) where L is the total population and k is an arbitrary constant.

For our problem, let us use L=250000, k=0.00003. Let the boundary condition be y(0)=250.

We shall use the odeint function to solve this ODE. As before, this function resides in a submodule of SciPy and doesn’t come along with Pylab. We import it,

from scipy.integrate import odeint

We can represent the given ODE as a Python function. This function takes the dependent variable y and the independent variable t as arguments and returns the ODE.

def epid(y, t):
    k = 0.00003
    L = 250000
    return k*y*(L-y)

Independent variable t can be assigned the values in the interval of 0 and 12 with 60 points using linspace:

In []: t = linspace(0, 12, 60)

Now obtaining the solution of the ode we defined, is as simple as calling the Python’s odeint function which we just imported

y = odeint(epid, 250, t)

We can plot the the values of y against t to get a graphical picture our ODE.

plot(y, t)

Let’s now try and solve an ODE of second order. Let’s take the example of a simple pendulum.

The equations can be written as a system of two first order ODEs

d(theta)/dt = omega

and

d(omega)/dt = - g/L sin(theta)

Let us define the boundary conditions as: at t = 0, theta = theta-naught = 10 degrees and omega = 0

Let us first define our system of equations as a Python function, pend_int. As in the earlier case of single ODE we shall use odeint function to solve this system of equations.

def pend_int(initial, t):
    theta = initial[0]
    omega = initial[1]
    g = 9.81
    L = 0.2
    f=[omega, -(g/L)*sin(theta)]
    return f

It takes two arguments. The first argument itself containing two dependent variables in the system, theta and omega. The second argument is the independent variable t.

In the function we assign the first and second values of the initial argument to theta and omega respectively. Acceleration due to gravity, as we know is 9.81 meter per second sqaure. Let the length of the the pendulum be 0.2 meter.

We create a list, f, of two equations which corresponds to our two ODEs, that is d(theta)/dt = omega and d(omega)/dt = - g/L sin(theta). We return this list of equations f.

Now we can create a set of values for our time variable t over which we need to integrate our system of ODEs. Let us say,

t = linspace(0, 20, 100)

We shall assign the boundary conditions to the variable initial.

initial = [10*2*pi/360, 0]

Now solving this system is just a matter of calling the odeint function with the correct arguments.

pend_sol = odeint(pend_int, initial,t)

plot(pend_sol)

This gives us 2 plots, omega vs t and theta vs t.

That brings us to the end of our discussion on ODEs and solving them in Python.

Using Python modules

We shall, in this section, see how to run Python scripts from the command line, and more details about importing modules.

Let us create a simple python script to print hello world. Open your text editor and type the following, and save the script as hello.py.

print "Hello world!"

Until now we have been running scripts from inside IPython, using

%run -i hello.py

But, as we know, IPython is just an advanced interpreter and it is not mandatory that every one who wants to run your program has IPython installed.

So, the right method to run the script, would be to run it using the Python interpreter. Open the terminal and navigate to the directory where hello.py is saved.

Then run the script by saying,

python hello.py

The script is executed and the string Hello World! has been printed.

Now let us run a script that plots a simple sine curve in the range -2pi to 2pi, using Python, instead of running it from within IPython.

Type the following lines and save it in sine_plot.py file.

x = linspace(-2*pi, 2*pi, 100)
plot(x, sin(x))
show()

Now let us run sine_plot.py as a python script.

python sine_plot.py

Do we get the plot? No! All we get is error messages. Python complains that linspace isn’t defined. What happened? Let us try to run the same script from within IPython. Start IPython, with the -pylab argument and run the script. It works!

What is going on, here? IPython when started with the -pylab option, is importing a whole set of functionality for us to work with, without bothering about importing etc.

Let us now try to fix the problem by importing linspace.

Let’s add the following line as the first line in the script.

from scipy import *

Now let us run the script again,

python sine_plot.py

Now we get an error, that says plot is not defined. Let’s edit the file and import this from pylab. Let’s add the following line, as the second line of our script.

from pylab import *

And run the script,

python sine_plot.py

Yes! it worked. So what did we do?

We actually imported the required functions and keywords, using import. By using the * , we are asking python to import everything from the scipy and pylab modules. This isn’t a good practice, as 1. it imports a lot of unnecessary things 2. the two modules may have functions with the same name, which would

cause a conflict.

One of the ways out is to import only the stuff that we need, explicitly.

Let us modify sine_plot.py by replacing the import statements with the following lines.

from scipy import linspace, pi, sin
from pylab import plot, show

Now let us try running the code again as,

python show_plot.py

As expected, this works. But, once the number of functions that you need to import increases, this is slightly inconvenient. Also, you cannot use functions with the same name, from different modules. To overcome this, there’s another method.

We could just import the modules as it is, and use the functions or other objects as a part of those modules. Change the import lines, to the following.

import scipy
import pylab

Now, replace pi with scipy.pi. Similarly, for sin and linspace. Replace plot and show with pylab.plot and pylab.show, respectively.

Now, run the file and see that we get the output, as expected.

We have learned how to import from modules, but what exactly is a module? How is one written?

A module is simply a Python script containing the definitions and statements, which can be imported. So, it is very easy to create your own modules. You just need to stick in all your definitions into your python file and put the file in your current directory.

When importing, Python searches the locations present in a variable called the Python path. So, our module file should be present in one of the locations in the Python path. The first location to be searched is the current directory of the script being run. So, having our module file in the current working directory, is the simplest way to allow it to be used as a module and import things from it.

Python has a very rich standard library of modules. It is very extensive, offering a wide range of facilities. Some of the standard modules are,

for Math: math, random for Internet access: urllib2, smtplib for System, Command line arguments: sys for Operating system interface: os for regular expressions: re for compression: gzip, zipfile, tarfile And there are lot more.

Find more information at Python Library reference, http://docs.python.org/library/

There are a lot of other modules like pylab, scipy, Mayavi, etc which are not part of the standard Python library.

This brings us to the end of our discussion on modules and running scripts from the command line.

Writing modules

In this section we shall look at writing modules, in some more detail. Often we will have to reuse the code that we have previously written. We do that by writing functions. Functions can then be put into modules, and imported as and when required.

Let us first write a function that computes the gcd of two numbers and save it in a script.

def gcd(a, b):
    while b:
        a, b = b, a%b
    return a

Now, we shall write a test function in the script that tests the gcd function, to see if it works.

if gcd(40, 12) == 4 and gcd(12, 13) == 1:
    print "Everything OK"
else:
    print "The GCD function is wrong"

Let us save the file as gcd_script.py in /home/fossee/gcd_script.py and run it.

$ python /home/fossee/gcd_script.py

We can see that the script is executed and everything is fine.

What if we want to use the gcd function in some of our other scripts. This is also possible since every python file can be used as a module.

But first, we shall understand what happens when you import a module.

Open IPython and type

import sys
sys.path

This is a list of locations where python searches for a module when it encounters an import statement. Hence, when we just did import sys, python searches for a file named sys.py or a folder named sys in all these locations one by one, until it finds one. We can place our script in any one of these locations and import it.

The first item in the list is an empty string which means the current working directory is also searched.

Alternatively, we can also import the module if we are working in same directory where the script exists.

Since we are in /home/fossee, we can simply do

import gcd_script

We can see that the gcd_script is imported. But the test code that we added at the end of the file is also executed.

But we want the test code to be executed only when the file is run as a python script and not when it is imported.

This is possible by using __name__ variable.

Go to the file and add

if __name__ == "__main__":

before the test code and indent the test code.

Let us first run the code.

$ python gcd_script.py

We can see that the test runs successfully.

Now we shall import the file

import gcd_script

We see that now the test code is not executed.

The __name__ variable is local to every module and it is equal to __main__ only when the file is run as a script. Hence, all the code that goes in to the if block, if __name__ == "__main__": is executed only when the file is run as a python script.

Object Oriented Programming

At the end of this section, you will be able to -

  • Understand the differences between Object Oriented Programming and Procedural Programming
  • Appreciate the need for Object Oriented Programming
  • Read and understand Object Oriented Programs
  • Write simple Object Oriented Programs

Suppose we have a list of talks, to be presented at a conference. How would we store the details of the talks? We could possibly have a dictionary for each talk, that contains keys and values for Speaker, Title and Tags of the talk. Also, say we would like to get the first name of the speaker and the tags of the talk as a list, when required. We could do it, as below.

talk = {'Speaker': 'Guido van Rossum',
        'Title': 'The History of Python'
        'Tags': 'python,history,C,advanced'}

def get_first_name(talk):
    return talk['Speaker'].split()[0]

def get_tags(talk):
    return talk['Tags'].split(',')

This is fine, when we have a small number of talks and a small number of operations that we wish to perform. But, as the number of talks increases, this gets inconvenient. Say, you are writing another function in some other module, that uses this talk dictionary, you will also need to pass the functions that act on talk to that function. This gets quite messy, when you have a lot of functions and objects.

This is where Objects come in handy. Objects, essentially, group data with the methods that act on the data into a single block/object.

Objects and Methods

The idea of objects and object oriented programming is being introduced, now, but this doesn’t mean that we haven’t come across objects. Everything in Python is an object. We have been dealing with objects all the while. Strings, lists, functions and even modules are objects in Python. As we have seen, objects combine data along with functions that act upon the data and we have been seeing this since the beginning. The functions that are tied to an object are called methods. We have seen various methods, until now.

s = "Hello World"
s.lower()

l = [1, 2, 3, 4, 5]
l.append(6)

lower is a string method and is being called upon s, which is a string object. Similarly, append is a list method, which is being called on the list object l.

Functions are also objects and they can be passed to and returned from functions, as we have seen in the SciPy section.

Objects are also useful, because the provide a similar interface, without us needing to bother about which exact type of object we are dealing with. For example, we can iterate over the items in a sequence, as shown below without really worrying about whether it’s a list or a dictionary or a file-object.

for element in (1, 2, 3):
    print element
for key in {'one':1, 'two':2}:
    print key
for char in "123":
    print char
for line in open("myfile.txt"):
    print line
for line in urllib2.urlopen('http://site.com'):
    print line

All objects providing a similar inteface can be used the same way.

Classes

When we created a string s, we obtained an object of type string.

s = "Hello World"
type(s)

s already comes with all the methods of strings. So, it suggests that there should be some template, based on which the object s is built. This template or blueprint to build an object is the Class. The class definition gives the blueprint for building objects of that kind, and each object is an instance of that class.

As you would’ve expected, we can define our own classes in Python. Let us define a simple Talk class for the example that we started this section with.

class Talk:
    """A class for the Talks."""

    def __init__(self, speaker, title, tags):
        self.speaker = speaker
        self.title = title
        self.tags = tags

    def get_speaker_firstname(self):
        return self.speaker.split()[0]

    def get_tags(self):
        return self.tags.split(',')

The above example introduces a lot of new things. Let us look at it, piece by piece.

A class is defined using a class block – the keyword class followed by the class name, in turn followed by a semicolon. All the statements within the class are enclosed in it’s block. Here, we have defined a class named Talk.

Our class has the same two functions that we had defined before, to get the speaker firstname and the tags. But along with that, we have a new function __init__. We will see, what it is, in a short while. By the way, the functions inside a class are called methods, as you already know. By design, each method of a class requires to have the same instance of the class, (i.e., the object) from which it was called as the first argument. This argument is generally defined using self.

self.speaker, self.title and self.tags are variables that contain the respective data. So, as you can see, we have combined the data and the methods operating on it, into a single entity, an object.

Let’s now initialize a Talk which is equivalent to the example of the talk, that we started with. Initializing an object is similar to calling a function.

bdfl = Talk('Guido van Rossum',
            'The History of Python',
            'python,history,C,advanced')

We pass the arguments of the __init__ function to the class name. We are creating an object bdfl, that is an instance of the class Talk and represents the talk by Guido van Rossum on the History of Python. We can now use the methods of the class, using the dot notation, that we have been doing all the while.

bdfl.get_tags()
bdfl.get_speaker_firstname()

The __init__ method is a special method, that is called, each time an object is created from a class, i.e., an instance of a class is created.

print bdfl.speaker
print bdfl.tags
print bdfl.title

As you can see, the __init__ method was called and the variables of the bdfl object have been set. object have been set. Also notice that, the __init__ function takes 4 arguments, but we have passed only three. The first argument self as we have already seen, is a reference to the object itself.

Inheritance

Now assume that we have a different category for Tutorials. They are almost like talks, except that they can be hands-on or not. Now, we do not wish to re-write the whole code that we wrote for the Talk class. Here, the idea of inheritance comes in handy. We “inherit” the Talk class and modify it to suit our needs.

class Tutorial(Talk):
    """A class for the tutorials."""

    def __init__(self, speaker, title, tags, handson=True):
        Talk.__init__(self, speaker, title, tags)
        self.handson = handson

    def is_handson(self):
        return self.handson

We have now derived the Tutorial class from the Talk class. The Tutorial class, has a different __init__ method, and a new is_handson method. But, since it is derived from the Talk method it also has the methods, get_tags and get_speaker_firstname. This concept of inheriting methods and values is called inheritance.

numpy = Tutorial('Travis Oliphant', 'Numpy Basics', 'numpy,python,beginner')
numpy.is_handson()
numpy.get_speaker_firstname()

As you can see, it has saved a lot of code duplication and effort.

That brings us to the end of the section on Object Oriented Programming. In this section we have learnt,

  • the fundamental difference in paradigm, between Object Oriented Programming and Procedural Programming
  • to write our own classes
  • to write new classes that inherit from existing classes

Functional programming

In this section, we shall look at a different style of programming called Functional programming, which is supported by Python.

The fundamental idea behind functional programming is that the output of a function, should not depend on the state of the program. It should only depend on the inputs to the program and should return the same output, whenever the function is called with the same arguments. Essentially, the functions of our code, behave like mathematical functions, where the output is dependent only on the variables on which the function depends. There is no “state” associated with the program.

Let’s look at some of the functionality that Python provides for writing code in the Functional approach.

List Comprehensions

Let’s take a very simple problem to illustrate how list comprehensions work. Let’s say, we are given a list of weights of people, and we need to calculate the corresponding weights on the moon. Essentially, return a new list, with each of the values divided by 6.0.

It’s a simple problem, and let’s first solve it using a for loop. We shall initialize a new empty, list and keep appending the new weights calculated

weights = [30, 45.5, 78, 81, 55.5, 62.5]

weights_moon = []

for w in weights:
    weights_moon.append(w/6.0)

print weights_moon

We had to initialize an empty list, and write a for loop that appends each of the values to the new list. List comprehensions provide a way of writing a one liner, that does the same thing, without resorting to initializing a new empty list. Here’s how we would do it, using list comprehensions

[ w/6.0 for w in weights ]

As we can see, we have the same list, that we obtained previously, using the for loop.

Let’s now look at a slightly more complex example, where we wish to calculate the weights on the moon, only if the weight on the earth is greater than 50.

weights = [30, 45.5, 78, 81, 55.5, 62.5]
weights_moon = []
for w in weights:
    if w > 50:
        weights_moon.append(w/6.0)

print weights_moon

The for loop above can be translated to a list comprehension as shown

[ w/6.0 for w in weights if w>50 ]

Now, let’s change the problem again. Say, we wish to give the weight on the moon (divide by 6), when the weight is greater than 50, otherwise triple the weight.

weights = [30, 45.5, 78, 81, 55.5, 62.5]
weights_migrate = []
for w in weights:
    if w > 50:
        weights_migrate.append(w/6.0)
    else:
        weights_migrate.append(w * 3.0)

print weights_migrate

This problem, however, cannot be translated into a list comprehension. We shall use the map built-in, to solve the problem in a functional style.

map

But before getting to the more complex problem, let’s solve the easier problem of returning the weight on moon for all the weights.

The map function essentially allows us to apply a function to all the elements of a list, and returns a new list that consists of the outputs from each of the call. So, to use map we need to define a function, that takes an input and returns a sixth of it.

def moonize(weight):
    return weight/6.0

Now, that we have our moonize function, we can pass the function and the list of weights as an argument to map and get the required list.

map(moonize, weights)

Let’s define a new function, that compares the weight value with 50 and returns a new value based on the condition.

def migrate(weight):
    if weight < 50:
        return weight*3.0
    else:
        return weight/6.0

Now, we can use map to obtain the new weight values

map(migrate, weights)

As you can see, we obtained the result, that we obtained before.

But, defining such functions each time, is slightly inconvenient. We are not going to use these functions at any later point, in our code. So, it is slightly wasteful to define functions for this. Wouldn’t it be nice to write them, without actually defining functions? Enter lambda!

lambda

Essentially, lambda allows us to define small functions anonymously. The first example that uses the moonize function can now be re-written as below, using the lambda function.

map(lambda x: x/6.0, weights)

lambda above is defining a function, which takes one argument x and returns a sixth of that argument. The lambda actually returns a function object which we could in fact assign to a name and use later.

l_moonize = lambda x: x/6.0

The l_moonize function, now behaves similarly to the moonize function.

The migrate function can be written using a lambda function as below

l_migrate = lambda x: x*3.0 if x < 50 else x/6.0

If you observed, we have carefully avoided the discussion of the example where only the weights that were above 50 were calculated and returned. This is because, this cannot be done using map. We may return None instead of calculating a sixth of the element, but we cannot ensure that the element is not present in the new list.

This can be done using filter and map in conjunction.

filter

The filter function, like the map takes two arguments, a function and a sequence and calls the function for each element of the sequence. The output of filter is a sequence consisting of elements for which the function returned True

The problem of returning a sixth of only those weights which are more than 50, can be solved as below

map(lambda x: x/6.0, filter(lambda x: x > 50, weights))

The filter function returns a list containing only the values which are greater than 50.

filter(lambda x: x > 50, weights)

reduce

As, the name suggests, reduce reduces a sequence to a single object. reduce takes two arguments, a function and a sequence, but the function should take two arguments as input.

reduce initially passes the first two elements as arguments, and continues iterating of the sequence and passes the output of the previous call with the current element, as the arguments to the function. The final result therefore, is just a single element.

reduce(lambda x,y: x*y, [1, 2, 3, 4])

multiplies all the elements of the list and returns 24.

Sage notebook

In this section, we will learn what Sage is, what is Sage notebook, how to start and use the sage notebook. In the notebook we will be specifically learning how to execute our code, how to write annotations and other content, typesetting the content and how to use the offline help available.

To start with, What is Sage? Sage is a free, open-source mathematical software. Sage can do a lot of math stuff for you including, but not limited to, algebra, calculus, geometry, cryptography, graph theory among other things. It can also be used as aid in teaching and research in any of the areas that Sage supports. So let us start Sage now

We are assuming that you have Sage installed on your computer now. If not please visit the page http://sagemath.org/doc/tutorial/introduction.html#installation for the tutorial on how to install Sage.

Let us now learn how to start Sage. On the terminal type

sage

This should start a new Sage shell with the prompt sage: which looks like this

sage:

So now we can type all the commands that Sage supports here. But Sage comes bundled with a much more elegant tool called Sage Notebook. It provides a web based user interface to use Sage. So once we have a Sage notebook server up and running, all we need is a browser to access the Sage functionality.

We could also use a server run by somebody, else like the official instance of the Sage server at http://sagenb.org. So all we need is just a browser, a modern browser, to use Sage and nothing else! The Sage notebook also provides a convenient way of sharing and publishing our work, which is very handy for research and teaching.

However we can also run our own instances of Sage notebook servers on all the computers we have a local installation of Sage. To start the notebook server just type

notebook()

on the Sage prompt. This will start the Sage Notebook server. If we are starting the notebook server for the first time, we are prompted to enter the password for the admin. Type the password and make a note of it. After this Sage automatically starts a browser page, with the notebook opened.

If it doesn’t automatically start a browser page check if the Notebook server started and there were no problems. If so open your browser and in the address bar type the URL shown in the instructions upon running the notebook command on the sage prompt.

If you are not logged in yet, it shows the Notebook home page and textboxes to type the username and the password. You can use the username ‘admin’ and the password you gave while starting the notebook server for the first time. There are also links to recover forgotten password and to create new accounts.

Once we are logged in with the admin account we can see the notebook admin page. A notebook can contain a collection of Sage Notebook worksheets.

Worksheet is basically a working area. This is where we enter all the Sage commands on the notebook.

The admin page lists all the worksheets created. On the topmost part of this page we have the links to various pages.

The home link takes us to the admin home page. The published link takes us to the page which lists all the published worksheets. The log link has the complete log of all the actions we did on the notebook. We have the settings link where we can configure our notebook, the notebook server, create and mangage accounts. We have a link to help upon clicking opens a new window with the complete help of Sage. The entire documentation of Sage is supplied with Sage for offline reference and the help link is the way to get into it. Then we can report bugs about Sage by clicking on Report a Problem link and there is a link to sign out of the notebook.

We can create a new worksheet by clicking New Worksheet link

Sage prompts you for a name for the worksheet. Let us name the worksheet as nbtutorial. Now we have our first worksheet which is empty.

A worksheet will contain a collection of cells. Every Sage command must be entered in this cell. Cell is equivalent to the prompt on console. When we create a new worksheet, to start with we will have one empty cell. Let us try out some math here

2 + 2
57.1 ^ 100

The hat operator is used for exponentiation. We can observe that only the output of the last command is displayed. By default each cell displays the result of only the last operation. We have to use print statement to display all the results we want to be displayed.

Now to perform more operations we want more cells. So how do we create a new cell? It is very simple. As we hover our mouse above or below the existing cells we see a blue line, by clicking on this new line we can create a new cell.

We have a cell, we have typed some commands in it, but how do we evaluate that cell? Pressing Shift along with Enter evaluates the cell. Alternatively we can also click on the evaluate link to evaluate the cell

matrix([[1,2], [3,4]])^(-1)

After we create many cells, we may want to move between the cells. To move between the cells use Up and Down arrow keys. Also clicking on the cell will let you edit that particular cell.

To delete a cell, clear the contents of the cell and hit backspace.

If you want to add annotations in the worksheet itself on the blue line that appears on hovering the mouse around the cell, Hold Shift and click on the line. This creates a What You See Is What You Get cell.

We can make our text here rich text. We can make it bold, Italics, we can create bulleted and enumerated lists in this area

This text contains both the **bold** text and also *italicised*
text.
It also contains bulleted list:
* Item 1
* Item 2
It also contains enumerate list:
1. Item 1
2. Item 2

In the same cell we can display typeset math using the LaTeX like syntax

$\int_0^\infty e^{-x} \, dx$

We enclose the math to be typeset within $ and $ or $$ and $$ as in LaTeX.

We can also obtain help for a particular Sage command or function within the worksheet itself by using a question mark following the command

sin?

Evaluating this cell gives me the entire help for the sin function inline on the worksheet itself. Similarly we can also look at the source code of each command or function using double question mark

matrix??

Sage notebook also provides the feature for autocompletion. To auto-complete a command type first few unique characters and hit tab key

sudo<tab>

To see all the commands starting with a specific name type those characters and hit tab

plo<tab>

To list all the methods that are available for a certain variable or a data-type we can use the variable name followed by the dot to access the methods available on it and then hit tab

s = 'Hello'
s.rep<tab>

The output produced by each cell can be one of the three states. It can be either the full output, or truncated output or hidden output. The output area will display the error if the Sage code we wrote in the cell did not successfully execute

a, b = 10

The default output we obtained now is a truncated output. Clicking at the left of the output area when the mouse pointer turns to hand gives us the full output, clicking again makes the output hidden and it cycles.

Lastly, Sage supports a variety of languages and each cell on the worksheet can contain code written in a specific language. It is possible to instruct Sage to interpret the code in the language we have written. This can be done by putting percentage sign(%) followed by the name of the language. For example, to interpret the cell as Python code we put

%python

as the first line in the cell. Similarly we have: %sh for shell scripting, %fortran for Fortran, %gap for GAP and so on. Let us see how this works. Say we have an integer. The type of the integer in default Sage mode is

a = 1
type(a)

Output: <type 'sage.rings.integer.Integer'>

We see that Integers are Sage Integers. Now let us put %python as the first line of the cell and execute the same code snippet

%python
a = 1
type(a)

Output: <type 'int'>

Now we see that the integer is a Python integer. Why? Because now we instructed Sage to interpret that cell as Python code.

This brings us to the end of the section on using Sage.

Symbolics

In this section, we shall learn to define symbolic expressions in Sage, use built-in constants and functions, perform integration and differentiation using Sage, define matrices and symbolic functions and simplify and solve them.

In addtion to a lot of other things, Sage can do Symbolic Math and we shall start with defining symbolic expressions in Sage.

On the sage notebook type

sin(y)

It raises a name error saying that y is not defined. We need to declare y as a symbol. We do it using the var function.

var('y')

Now if you type

sin(y)

Sage simply returns the expression.

Sage treats sin(y) as a symbolic expression. We can use this to do symbolic math using Sage’s built-in constants and expressions.

Let us try out a few examples.

var('x, alpha, y, beta')
(x^2/alpha^2)+(y^2/beta^2)

We have defined 4 variables, x, y, alpha and beta and have defined a symbolic expression using them.

Here is an expression in theta

var('theta')
sin(theta)*sin(theta)+cos(theta)*cos(theta)

Sage also provides built-in constants which are commonly used in mathematics, for instance pi, e, infinity. The function n gives the numerical values of all these constants.

n(pi)
n(e)
n(oo)

If you look into the documentation of function n by doing

n?

You will see what all arguments it takes and what it returns. It will be very helpful if you look at the documentation of all functions introduced, in this section.

Also we can define the number of digits we wish to have in the constants. For this we have to pass an argument – digits.

n(pi, digits = 10)

Apart from the constants Sage also has a lot of built-in functions like sin, cos, log, factorial, gamma, exp, arcsin etc ...

Lets try some of them out on the Sage notebook.

sin(pi/2)

arctan(oo)

log(e,e)

Given that we have defined variables like x, y etc., we can define an arbitrary function with desired name

var('x')
function('f',x)

Here f is the name of the function and x is the independent variable . Now we can define f(x) to be

f(x) = x/2 + sin(x)

Evaluating this function f for the value x=pi returns pi/2.

f(pi)

We can also define functions that are not continuous but defined piece-wise. Let us define a function which is a parabola between 0 to 1 and a constant from 1 to 2 .

var('x')
h(x)=x^2
g(x)=1

f=Piecewise([[(0,1),h(x)],[(1,2),g(x)]],x)
f

We can also define functions convergent series and other series.

We first define a function f(n) in the way discussed above.

var('n')
function('f', n)

To sum the function for a range of discrete values of n, we use the sage function sum.

For a convergent series , f(n)=1/n^2 we can say

var('n')
function('f', n)
f(n) = 1/n^2
sum(f(n), n, 1, oo)

Lets us now try another series

f(n) = (-1)^(n-1)*1/(2*n - 1)
sum(f(n), n, 1, oo)

This series converges to pi/4.

Let’s now look at some calculus.

diff(x**2 + sin(x), x)

The diff function differentiates an expression or a function. It’s first argument is expression or function and second argument is the independent variable.

We have already tried an expression now lets try a function

f = exp(x^2) + arcsin(x)
diff(f(x), x)

To get a higher order differential we need to add an extra third argument for order

diff(f(x), x, 3)

in this case it is 3.

Just like differentiation of expression you can also integrate them

x = var('x')
s = integral(1/(1 + (tan(x))**2),x)
s

Many a times we need to find factors of an expression, we can use the “factor” function

y = (x^100 - x^70)*(cos(x)^2 + cos(x)^2*tan(x)^2)
f = factor(y)

One can simplify complicated expression

f.simplify_full()

This simplifies the expression fully. We can also do simplification of just the algebraic part and the trigonometric part

f.simplify_exp()
f.simplify_trig()

One can also find roots of an equation by using find_root function

phi = var('phi')
find_root(cos(phi)==sin(phi),0,pi/2)

Let’s substitute this solution into the equation and see we were correct

var('phi')
f(phi)=cos(phi)-sin(phi)
root=find_root(f(phi)==0,0,pi/2)
f.substitute(phi=root)

as we can see when we substitute the value the answer is almost = 0 showing the solution we got was correct.

Lets us now try some matrix algebra symbolically

var('a,b,c,d')
A=matrix([[a,1,0],[0,b,0],[0,c,d]])
A

Now lets do some of the matrix operations on this matrix

::
A.det() A.inverse()

That brings us to the end of our discussion on symbolics.

Plotting using Sage

In this section we shall look at

  • 2D plotting in SAGE
  • 3D plotting in SAGE

We shall first create a symbolic variable x

x = var('x')

We shall plot the function sin(x) - cos(x) ^ 2 in the range (-5, 5).

plot(sin(x) - cos(x) ^ 2, (x, -5, 5))

As we can see, the plot is shown.

plot command takes the symbolic function as the first argument and the range as the second argument.

We have seen that plot command plots the given function on a linear range.

What if the x and y values are functions of another variable. For instance, lets plot the trajectory of a projectile.

A projectile was thrown at 50 m/s^2 and at an angle of 45 degrees from the ground. We shall plot the trajectory of the particle for 5 seconds.

These types of plots can be drawn using the parametric_plot function. We first define the time variable.

t = var('t')

Then we define the x and y as functions of t.

f_x = 50 * cos(pi/4)
f_y = 50 * sin(pi/4) * t - 1/2 * 9.81 * t^2 )

We then call the parametric_plot function as

parametric_plot((f_x, f_y), (t, 0, 5))

And we can see the trajectory of the projectile.

The parametric_plot funciton takes a tuple of two functions as the first argument and the range over which the independent variable varies as the second argument.

Now we shall look at how to plot a set of points.

We have the line function to achieve this.

We shall plot sin(x) at few points and join them.

First we need the set of points.

points = [ (x, sin(x)) for x in srange(-2*float(pi), 2*float(pi), 0.75) ]

srange takes a start, a stop and a step argument and returns a list of point. We generate list of tuples in which the first value is x and second is sin(x).

line(points)

plots the points and joins them with a line.

The line function behaves like the plot command in matplotlib. The difference is that plot command takes two sequences while line command expects a sequence of co-ordinates.

As we can see, the axes limits are set by SAGE. Often we would want to set them ourselves. Moreover, the plot is shown here since the last command that is executed produces a plot.

Let us try this example

plot(cos(x), (x,0,2*pi))
# Does the plot show up??

As we can see here, the plot is not shown since the last command does not produce a plot.

The actual way of showing a plot is to use the show command.

p1 = plot(cos(x), (x,0,2*pi))
show(p1)
# What happens now??

As we can see the plot is shown since we used it with show command.

show command is also used set the axes limits.

p1 = plot(cos(x), (x,0,2*pi))
show(p1, xmin=0, xmax=2*pi, ymin=-1.2, ymax=1.2)

As we can see, we just have to pass the right keyword arguments to the show command to set the axes limits.

The show command can also be used to show multiple plots.

p1 = plot(cos(x), (x, 0, 2*pi))
p2 = plot(sin(x), (x, 0, 2*pi))
show(p1+p2)

As we can see, we can add the plots and use them in the show command.

Now we shall look at 3D plotting in SAGE.

We have the plot3d function that takes a function in terms of two independent variables and the range over which they vary.

x, y = var('x y')
plot3d(x^2 + y^2, (x, 0, 2), (y, 0, 2))

We get a 3D plot which can be rotated and zoomed using the mouse.

parametric_plot3d function plots the surface in which x, y and z are functions of another variable.

u, v = var("u v")
f_x = u
f_y = v
f_z = u^2 + v^2
parametric_plot3d((f_x, f_y, f_z), (u, 0, 2), (v, 0, 2))

Using Sage

In this section we shall quickly look at a few examples of using Sage for Linear Algebra, Calculus, Graph Theory and Number theory.

Let us begin with Calculus. We shall be looking at limits, differentiation, integration, and Taylor polynomial.

To find the limit of the function x*sin(1/x), at x=0, we say

lim(x*sin(1/x), x=0)

We get the limit to be 0, as expected.

It is also possible to the limit at a point from one direction. For example, let us find the limit of 1/x at x=0, when approaching from the positive side.

lim(1/x, x=0, dir='above')

To find the limit from the negative side, we say,

lim(1/x, x=0, dir='below')

Let us now see how to differentiate, using Sage. We shall find the differential of the expression exp(sin(x^2))/x w.r.t x. We shall first define the expression, and then use the diff function to obtain the differential of the expression.

var('x')
f = exp(sin(x^2))/x

diff(f, x)

We can also obtain the partial differentiation of an expression w.r.t one of the variables. Let us differentiate the expression exp(sin(y - x^2))/x w.r.t x and y.

var('x y')
f = exp(sin(y - x^2))/x

diff(f, x)

diff(f, y)

Now, let us look at integration. We shall use the expression obtained from the differentiation that we did before, diff(f, y)e^(sin(-x^2 + y))*cos(-x^2 + y)/x. The integrate command is used to obtain the integral of an expression or function.

integrate(e^(sin(-x^2 + y))*cos(-x^2 + y)/x, y)

We get back the correct expression. The minus sign being inside or outside the sin function doesn’t change much.

Now, let us find the value of the integral between the limits 0 and pi/2.

integral(e^(sin(-x^2 + y))*cos(-x^2 + y)/x, y, 0, pi/2)

Let us now see how to obtain the Taylor expansion of an expression using sage. Let us obtain the Taylor expansion of (x + 1)^n up to degree 4 about 0.

var('x n')
taylor((x+1)^n, x, 0, 4)

This brings us to the end of the features of Sage for Calculus, that we will be looking at. For more, look at the Calculus quick-ref from the Sage Wiki.

Next let us move on to Matrix Algebra.

Let us begin with solving the equation Ax = v, where A is the matrix matrix([[1,2],[3,4]]) and v is the vector vector([1,2]).

To solve the equation, Ax = v we simply say

x = solve_right(A, v)

To solve the equation, xA = v we simply say

x = solve_left(A, v)

The left and right here, denote the position of A, relative to x.

Now, let us look at Graph Theory in Sage.

We shall look at some ways to create graphs and some of the graph families available in Sage.

The simplest way to define an arbitrary graph is to use a dictionary of lists. We create a simple graph by

G = Graph({0:[1,2,3], 2:[4]})

We say

G.show()

to view the visualization of the graph.

Similarly, we can obtain a directed graph using the DiGraph function.

G = DiGraph({0:[1,2,3], 2:[4]})

Sage also provides a lot of graph families which can be viewed by typing graph.<tab>. Let us obtain a complete graph with 5 vertices and then show the graph.

G = graphs.CompleteGraph(5)

G.show()

Sage provides other functions for Number theory and Combinatorics. Let’s have a glimpse of a few of them.

prime_range(100, 200)

gives primes in the range 100 to 200.

is_prime(1999)

checks if 1999 is a prime number or not.

factor(2001)

gives the factorized form of 2001.

C = Permutations([1, 2, 3, 4])
C.list()

gives the permutations of [1, 2, 3, 4]

C = Combinations([1, 2, 3, 4])
C.list()

gives all the combinations of [1, 2, 3, 4]

That brings us to the end of this session showing various features available in Sage.

Using Sage to Teach

In this section, we shall look at

  • How to use the “@interact” feature of SAGE for better demonstration
  • How to use SAGE for collaborative learning

Let us look at a typical example of demonstrating a damped oscillation.

t = var('t')
p1 = plot( e^(-t) * sin(2*t), (t, 0, 15))
show(p1)

Now let us reduce the damping factor

t = var('t')
p1 = plot(e^(-t/2) * sin(2*t), (t, 0, 15))
show(p1)

Now if we want to reduce the damping factor even more, we would be using e^(-t/3). We can observe that every time we have to change, all we do is change something very small and re evaluate the cell.

This process can be simplified, using the @interact feature of SAGE.

@interact
def plot_damped(n=1):
    t = var('t')
    p1 = plot( e^(-t/n) * sin(2*t), (t, 0, 20))
    show(p1)

We can see that the function is evaluated and the plot is shown. We can also see that there is a field to enter the value of n and it is currently set to 1. Let us change it to 2 and hit enter.

We see that the new plot with reduced damping factor is shown. Similarly we can change n to any desired value and hit enter and the function will be evaluated.

This is a very handy tool while demonstrating or teaching.

Often we would want to vary a parameter over range instead of taking it as an input from the user. For instance we do not want the user to give n as 0 for the damping oscillation we discussed. In such cases we use a range of values as the default argument.

@interact
def plot_damped(n=(1..10)):
    t = var('t')
    p1 = plot( e^(-t/n) * sin(2*t), (t, 0, 20))
    show(p1)

We get similar plot but the only difference is the input widget. Here it is a slider unlike an input field. We can see that as the slider is moved, the function is evaluated and plotted accordingly.

Sometimes we want the user to have only a given set of options. We use a list of items as the default argument in such situations.

@interact
def str_shift(s="STRING", shift=(0..8), direction=["Left", "Right"]):
    shift_len = shift % len(s)
    chars = list(s)
    if direction == "Right":
        shifted_chars = chars[-shift_len:] + chars[:-shift_len]
    else:
        shifted_chars = chars[shift_len:] + chars[:shift_len]
    print "Actual String:", s
    print "Shifted String:", "".join(shifted_chars)

We can see that buttons are displayed which enables us to select from a given set of options.

We have learnt how to use the @interact feature of SAGE for better demonstration. We shall look at how to use SAGE worksheets for collaborative learning.

The first feature we shall see is the publish feature. Open a worksheet and in the top right, we can see a button called publish. Click on that and we get a confirmation page with an option for re publishing.

For now lets forget that option and simply publish by clicking yes. The worksheet is now published.

Now lets sign out and go to the sage notebook home. We see link to browse published worksheets. Lets click on it and we can see the worksheet. This does not require login and anyone can view the worksheet.

Alternatively, if one wants to edit the sheet, there is a link on top left corner that enables the user to download a copy of the sheet onto their home. This way they can edit a copy of the worksheet.

We have learnt how to publish the worksheets to enable users to edit a copy. Next, we shall look at how to enable users to edit the actual worksheet itself.

Let us open the worksheet and we see a link called share on the top right corner of the worksheet. Click the link and we get a box where we can type the usernames of users whom we want to share the worksheet with. We can even specify multiple users by seperating their names using commas. Once we have shared the worksheet, the worksheet appears on the home of shared users.

ReST

The Pythonic way to Document

What is ReST?

  • Sage Docs - Linear Algebra | Sources
  • ReST is a lightweight markup language.
  • Developed by the Python community
  • Heavily used in documentation of Python code

Why ReST?

  • Highly readable source format.
  • Easy to learn and write.
  • Simple yet Powerful markup that can produce html/LaTeX output.

Tools used

  • You will need python-docutils to compile your ReST documents
  sudo apt-get install python-docutils

+ To build websites, look at Sphinx, rest2web, etc.
+ rst2html doesn't support math. We shall use Sphinx to see how math works
sudo apt-get install python-sphinx

Generating output

  • html:

    rst2html source.rst destination.html
  • LaTeX

    rst2latex source.rst destination.tex
  • Slide shows

    rst2s5 source.rst destination.html

Paragraph

  • The most basic structural element

  • Just a chunk of text separated by blank lines.

  • Must begin at the left edge margin.

  • Indented paragraphs are output as quoted text.

  • Example

    Sage provides standard constructions from linear algebra,
    e.g., the characteristic polynomial, echelon form, trace,
    decomposition, etc., of a matrix.

Section Headings

  • Single line of text with adornment

  • Adornment is underline alone or an over- and underline

  • - = ` :  ' " ~  ^ _ *  + # <  > can be used for adornment

  • All headings with same adornment are at same level

  • NOTE- The over/under line should be at least as long as the heading

    Linear Algebra
    ==============
    
    Matrix spaces
    -------------
    
    Sparse Linear Algebra
    ---------------------

Text Styles

Markup text Resulting text
*italics* italics
**strong** strong
``inline literal`` inline literal

Code Samples

  • To include code, end the prior paragraph with =::=

  • Code needs to be indented one level

  • The code block ends, when text falls back to the previous indentation

  • For instance

    For example, each of the following is legal ::
    
      plot(x, y)         # plot x and y using default line style and color
      plot(x, y, 'bo')   # plot x and y using blue circle markers
      plot(y)            # plot y using x as index array 0..N-1
      plot(y, 'r+')      # ditto, but with red plusses

Math

  • ReST in itself doesn’t support math

  • Sphinx has support for math using ~jsmath~ or ~pngmath~
    :math: `3 \times 3`
    
    .. math::
    
    \sum_{n=0}^N x_n = y

Lists

  • Three flavors - Enumerated, Bulleted, Definition

  • Always start as a new paragraph — preceeded by a new line

  • Enumerated

    Number or Letter followed by a =.=, =)= or surrounded by =( )=.

    1) Numbers
    #) auto numbered
    A. Upper case letters
    a) lower case letters
    i) roman numerals
    (I) more roman numerals

Lists ...

  • Bulleted lists

    Start a line with -, + or *

* a bullet point using "*"

  - a sub-list using "-"

    + yet another sub-list

  - another item

Lists ...

  • Definition Lists
    • Consist of Term, and it’s definition.
    • Term is one line phrase; Definition is one or more paragraphs
    • Definition is indented relative to the term
    • Blank lines are not allowed between term and it’s definition
what
Definition lists associate a term with a definition.

Tables

  • Simple Tables
    • Each line is a row.
    • The table ends with ~=~
    • Column Header is specified by using ~=~
    • Cells may span columns; ~-~ is used to specify cells spanning columns.
============ ============ ===========
 Header 1     Header 2     Header 3
============ ============ ===========
 body row 1   column 2     column 3
 body row 2   Cells may span columns.
------------ ------------------------
 body row 3   column 2     column 3
============ ============ ===========

Tables...

Grid Tables

+------------+------------+-----------+
| Header 1   | Header 2   | Header 3  |
+============+============+===========+
  body row 1   column 2     column 3
+------------+------------+-----------+
  body row 2   Cells may span columns.
+------------+------------+-----------+
  body row 3   Cells may    - Cells
+------------+ span rows.   - contain
  body row 4                - blocks.
+------------+------------+-----------+

Footnotes

This[#]_ gives auto-numbered[#]_ footnotes.

This[*]_ gives auto-symbol footnotes[*]_.

.. [#] First auto numbered footnote
.. [#] Second auto numbered footnote
.. [*] First auto symbol footnote
.. [*] Second auto symbol footnote

References

  • An Introduction to reStructured Text – David Goodger
  • Quick reStructuredText
  • reStructuredText– Bits and Pieces – Christoph Reller

Exercises

  1. Consider the iteration x n+1=f(xn) where f(x)=kx(1-x) Plot the successive iterates of this process.
  2. Plot this using a cobweb plot as follows:
    1. Start at (x0,0)
    2. Draw line to (xi,f(xi));
    3. Set x i+1=f(xi)
    4. Draw line to (xi,xi)
    5. Repeat from 2 for as long as you want
  3. Plot the Koch snowflake. Write a function to generate the necessary points given the two points constituting a line.
    1. Split the line into 4 segments.
    2. The first and last segments are trivial.
    3. To rotate the point you can use complex numbers, recall that ze jθ rotates a point z in 2D by θ.
    4. Do this for all line segments till everything is done.
  4. Show rate of convergence for a first and second order finite difference of sin(x)
  5. Given, the position of a projectile in in pos.txt plot it’s trajectory.
    • Label both the axes.
    • What kind of motion is this?
    • Title the graph accordingly.
    • Annotate the position where vertical velocity is zero.
  6. Write a Program that plots a regular n-gon(Let n = 5).
  7. Create a sequence of images in which the damped oscillator ( e -x/10sin(x)) slowly evolves over time.
  8. Given a list of numbers, find all the indices at which 1 is present. numbers = [1, 1, 3, 4, 3, 6, 7, 8, 1, 2, 4, 1]
  9. Given a list of numbers, find all the indices at which 1 is present. numbers = [1, 1, 3, 4, 3, 6, 7, 8, 1, 2, 4, 1]. Solve the problem using a functional approach.