Lecture 13  Numerical Python (numpy)  ipython environment  Plotting using matplotlib

Matrix-Matrix Multiply  In matrix computations, AB is the matrix product of matrix A

with B (NOT element-wise multiply)  If we multiply the following 2x2 matrices for example, a

b

e

f

c

d

g

h

=

ae+bg

af+bh

ce+dg

cf+dh

 Product of two matrices, A (pxq) and B(qxr), can be obtained

using dot function from numpy. This will result in a pxr matrix  How can we create these matrices? A

C

B

1

2

3

4

5

6

7

8

9

.

1

1

1

1

1

1

1

1

1

=

6

6

6

15

15

15

24

24

24

How to create the input matrices? C

B

A 1

2

3

4

5

6

7

8

9

.

1

1

1

1

1

1

1

1

1

=

6

6

6

15

15

15

24

24

24

>>> import numpy as np  Create the 2D array explicitly to create the A matrix

>>> A = np.array([[1,2,3],[4,5,6],[7,8,9]]);  Or we can create the 1D array and reshape it to a 3x3 2D array

>>> A = np.arange(1,10); #This creates 1D array >>> A = np.reshape(A,(3,3))

How do we create the above 3x3 B matrix? >>> B = np.array([[1,1,1],[1,1,1],[1,1,1]]); # OR >>> B = np.ones((3,3),dtype=int)  For matrix-multiplication: >>> C = np.dot(A,B);

random and fft modules from NumPy  We will use the random module from numpy, i.e., np.random

(Note: There is also a random module in standard Python) >>> dir(np.random) (or) >>> help(np.random)  Set the seed of the random number generator manually (this will

generate the same sequence of random numbers every time) >>> np.random.seed(12)  Generate n uniform numbers on [0.0,1.0)

>>> u = np.random.random(n)  Generate n uniform numbers on [low, high)

>>> u = np.random.uniform(-1, 1, n)  The random module provides more general distributions – normal, binomial, chi-square, exponential etc.,  The fft module from numpy is accessed using, i.e., np.fft >>> np.fft.fft(x) # FFT >>> np.fft.ifft(x) # Inverse FFT

linalg module from NumPy  Sub-package of NumPy (useful but limited, SciPy has more

functionalities) >>> import numpy as np  linalg module has core linear algebra functions  Vector or Matrix Norm, Matrix Inverse, Solve a linear system of Equations, Determinant, Cholesky Decomposition etc. >>> dir(np.linalg) >>> help(np.linalg) >>> help(np.linalg.norm) Example: Let’s consider the linear system of equations: 3x + y = 9 ---- (1) x + 2y = 8 ---- (2)

3

x

9 =

1 We will set this up as Ax = b

1 2

y

8

x A b Matrix of known Vector of Vector of coefficients unknowns constants

linalg module from NumPy >>> a = np.array([[3, 1], [4, 2]]) >>> b = np.array([9,8]) >>> x = np.linalg.solve(a,b) >>> x  Check that the solution is correct (we use the method all() to check if

all elements are True in the array): >>> np.dot(a,x) == b >>> (np.dot(a,x) == b).all()  Can we do this on floating-point? (No, we want to compare the

absolute value to a threshold)  To obtain the inverse, use the inv function:

>>> ainv = np.linalg.inv(a) >>> np.dot(a,ainv)

# and this should equal to what?

NumPy’s matrix object  Create a matrix object (i.e., an object of matrix class)

>>> help(np.matrix) >>> x1 = np.array([[3, 1], [1,2]])  Convert the array object to the matrix object >>> x2 = np.matrix(x1) >>> type(x2)  Transpose

>>> x3 = x2.T  Inverse

>>> x3 = x2.I  Matrix-Matrix Multiply: (* when used on matrix objects means

matrix-matrix multiplication) >>> x1*x1 # ndarray object >>> x2*x2 # matrix object

NumPy Data Types Some common data-types: boolean – True or False (stored as a byte) np.int – (Platform dependent integer either int32 or int64) np.int8 – (Byte, signed -128 to 127) np.uint8 – (Byte, unsigned, 0 to 255) np.int16 – 16-bit integer (16-bit, -32768 to 32767) np.uint16 – 16-bit integer (16-bit, 0 to 65535) np.int32 – 32-bit integer np.float32 – (single-precision float, sign bit, 8 bit exponent, 23 bits mantissa) np.float64 – (double-precision float, sign bit, 11 bits exponent, 52bit mantissa) Documentation: http://docs.scipy.org/doc/numpy/user/basics.types.html

matplotlib Plotting Library  Python 2D Plotting Library  pyplot provides a MATLAB like framework  To use we import numpy and pyplot import numpy as np import matplotlib.pyplot as plt

 Interactive matplotlib use with X-server (we will use Xming)  Non-interactive matplotlib without X-server (use a backend to render

the plot and save it as an image file) Documentation: http://matplotlib.sourceforge.net/contents.html http://matplotlib.sourceforge.net/users/pyplot_tutorial.html Book uses Easyviz which is an interface to other plotting packages http://code.google.com/p/scitools/wiki/EasyvizDocumentation

Simple Line Plot (example1.py)  Let’s start with a simple line plot, first import numpy and

matplotlib.pyplot import numpy as np import matplotlib.pyplot as plt  Create an array using the linspace function from numpy

(this will create a vector (1D array) with discrete data points)

y = np.linspace(0, 10, 11);  A line plot connects these data points. We will use the plot and show

functions to plot and show the plot (x will be the index array 0 to 10)

plt.plot(y); plt.show();

Line Plot – Completed Script (example1.py) # Program Name: example1.py # Import the necessary modules import numpy as np import matplotlib.pyplot as plt # Create the array y = np.linspace(0, 10, 11); # Generate the plot plt.plot(y); plt.show()

Generated Line Plot

Plotting Capabilities  What are some things we might want to do while plotting ?  Decorating the plot (adding title, axes labels, legend etc)  Exporting the plot to a portable format, e.g. ps, pdf, or png  Multiple plots on the same figure  Creating multiple figures from the same program  Generating subplots (plot separately, with more than one plot per

window)  Creating different kinds of plots (histograms, 3D plots, etc.)

Plot the sine function (example2.py)  Let’s plot y = sin(x) where x takes values in 0, 4π  We will use π and sine function from numpy

 Import numpy and matplotlib (using local names)

import numpy as np import matplotlib.pyplot as plt  Create the input array, x using linspace and use vectorization to

compute y a=0 b = 4*np.pi; dx = 0.1; x = np.linspace(a, b, (b-a)/dx+1); y = np.sin(x)  We will plot y versus x using the plot function from plt module

plt.plot(x, y)

Decorating the Plot (example2.py)  Add labels to the x- and y-axes

plt.xlabel(‘x’) plt.ylabel(‘y’)  Add a title to the plot

plt.title(‘Plot of sine function’)  Display legend in the figure (i.e., we use labels to distinguish multiple

datasets) plt.legend((’y=sin(x)',))  We can display a grid if we want

plt.grid(True);  Finally we will display the plot

plt.show()

Completed Script (example2.py) # Program Name: example2.py # Import the necessary modules import numpy as np import matplotlib.pyplot as plt # Create the x array and evaluate y a=0 b = 4*np.pi; dx = 0.1; x = np.linspace(a, b, (b-a)/dx+1); y = np.sin(x) # Generate the plot plt.plot(x,y); plt.xlabel(‘x’) plt.ylabel(‘y’) plt.title(‘Plot of sine function’) plt.legend((’y=sin(x)',)) plt.grid(True); plt.show()

Line Plot for y=sin(x)

Other options: colors, markers (example2.py)  Line plot using red solid line:

plt.plot(x, y, ‘r’)

 Blue dashed lines:

plt.plot(x, y, ‘--b’)

 Green dashed lines with dotted markers (o):

plt.plot(x, y, ‘--go’)  Plot just using green square (s) markers:

plt.plot(x, y, ‘gs’)

 Plotting discrete data using stem plot

(A stem plot displays data as lines that end with a marker symbol at each data value): plt.stem(x, y)

Multiple plots on the same figure (example3.py)  Copy example2.py to example3.py  Let’s plot y = sin(x) and y = cos(x)where x takes values in 0, 4π  We will use π and sin and cos function from numpy

 Import numpy and matplotlib (using local names)

import numpy as np import matplotlib.pyplot as plt  Create an array using linspace and use vectorization to compute y

a=0 b = 4*np.pi; dx = 0.1; x = np.linspace(a, b, (b-a)/dx+1); y1 = np.sin(x) y2 = np.cos(x)

Multiple plots on the same figure (example3.py)  We will plot y versus x using the plot function from plt module

plt.plot(x, y1, ‘-rs’, x, y2, ‘-bd’)  Add labels to the x- and y-axes

plt.xlabel(‘x’) plt.ylabel(‘y’)  Add a title to the plot

plt.title(‘Plot of sine and cosine function’)  Display legend in the figure (i.e., we use labels to distinguish multiple

datasets) plt.legend((’y=sin(x)',’y=cos(x)’))  We can display a grid if we want

plt.grid(True);  Finally we will display the plot

plt.show()

Multiple Figures (example4.py)  Copy example3.py to example4.py  Let’s plot y = sin(x) and y = xsin(x)

where x takes values in -4π, 4π  Import numpy and matplotlib (using local names) import numpy as np import matplotlib.pyplot as plt  Create an array using linspace and use vectorization to compute y1

and y2 a = -4*np.pi; b = 4*np.pi; dx = 0.1; x = np.linspace(a, b, (b-a)/dx+1); y1 = np.sin(x) y2 = x*np.sin(x)

Example 4 (example4.py)  We will plot y versus x using the plot function from plt module

plt.figure(1); plt.plot(x, y1, ‘r’) plt.figure(2); plt.plot(x, y2, ‘b’)  We can display a grid if we want

plt.grid(True);  Finally we will display the plot

plt.show()

Generating subplots (example5.py)  Use the subplot command:

subplot(rows, columns, figureindex)  Let’s plot y = sin(x) and y = xsin(x) where x takes values in -4π, 4π  Import numpy and matplotlib (using local names) import numpy as np import matplotlib.pyplot as plt  Create an array using linspace and use vectorization to compute y1

and y2 a = -4*np.pi; b = 4*np.pi; dx = 0.1; x = np.linspace(a, b, (b-a)/dx+1); y1 = np.sin(x) y2 = x*np.sin(x)

Example 5 (example5.py)  We want to create a subplot with 2 rows, 1 column  We will plot y versus x using the plot function from plt module

plt.subplot(2,1,1) plt.plot(x, y1, ‘r’) plt.subplot(2,1,2) plt.plot(x, y2, ‘b’)  We display the plot

plt.show()

Other options  If we want to save it to a file:

plt.savefig(“sine.png”)  We can adjust axes manually using

plt.axis([x_min, x_max, y_min, y_max])  We can add other options to the axis labels

plt.xlabel(‘x’, fontsize = 14, color=‘red’) plt.plot(x, y, linewidth=2.0)  We can also use Tex

plt.title(r’$\sin(2 \omega t’) Refer: http://matplotlib.sourceforge.net/users/mathtext.html

3D surface plot example using mplot3d +

 mplot3d module enables 3D plots using matplotlib (Mayavi is another

package that provides 3D plotting capabilities) import numpy as np import matplotlib.pyplot as plt  Additionally import mplot3d

import mpl_toolkits.mplot3d as p3  Create two 1D arrays that represent the x- and y- co-ordinates of a grid

x = np.linspace(-2, 2, 500); y = np.linspace(-2, 2, 500);  meshgrid creates a grid over which we can evaluate a function, i.e. it uses the 1D arrays, x, and y to constructs matrices, X and Y (lengthY, lengthX) X, Y = np.meshgrid(x,y);  Create a function, Z Z = X**2 + Y**2; Z = np.exp(-Z); +http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/mlab.html

3D Plot (an example)  Create a new figure (this returns an instance of the class,

matplotlib.figure.Figure) fig = plt.figure();  This creates a 3D axes object (we are plotting 3D objects on a 2D

matplotlib figure) ax = p3.Axes3D(fig);  Create a surface plot, using the plot_surface command (by default, it

plots using shades of a solid color, color mapping using cmap argument) ax.plot_surface(X, Y, Z);  Save the figure

plt.savefig("test.png"); plt.show()