Scientific Programming using Python
Python has excellent supprt for scientific computing. It offers a lot of packages for this purpose. There is also a significant amount of web support available for these packages. One of the most widely used packages is Scipy (pronounced "Sigh-Pie"). It has a number of core packages:
- Numpy (N-dimensional array package)
- Scipy Library (Fundamental library for scientific computing)
- Matplotlib (Comprehensive 2D Plotting)
- Ipython (Enhanced Interactive Console)
- Sympy (Symbolic mathematics)
- Pandas (Data structures & analysis)
- Nose (Framework for testing Python code)
Some of the packages built on Scipy:
- Chaco (Embedded, interactive plotting)
- Mayavi (3D visualization)
- Cython (build C/C++ libraries)
- Scikits (scikit-image for image processing, scikit-learn for machine learning)
If you are interested you should also take a look at the list of topical software available here.
In this lecture, we will cover the basics of Numpy, Scipy Library and Matplotlib. A list of some good tutorials on Scipy core/based packages are given below:
Creating Numpy Arrays¶
We can create numpy arrays as follows:
import numpy as np
A = np.array([1,2,3,4,5])
#Creates a list and then converts it to an array
print 'A =',A
print 'type:',type(A)
print 'data type:', A.dtype
print 'dimensions:', A.ndim
print 'Shape',A.shape
There are some functions to create arrays that are frequently required:
print "Zeros"
Z = np.zeros(2)
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "More Zeros"
Z = np.zeros((1,2))
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Ones"
Z = np.ones((3,2))
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Identity"
Z = np.eye(3,3)
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Integer Range"
Z = np.arange(0,10,2)
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Linear Spaced"
Z = np.linspace(0,1,10)
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Random floating point numbers"
Z = np.random.random((3,3))
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
print "Random Integers"
Z = np.random.randint(0,5,(3,3))
print 'Z =',Z
print 'type:',type(Z)
print 'data type:', Z.dtype
print 'dimensions:', Z.ndim
print 'Shape',Z.shape
print '\n'+'='*10+'\n'
Scalar Operations
Scalar operations work as expected!
A = np.zeros((3,4))#similarly 'ones'
A+=1 #Le Addition
print A
A = 2*A
print A
A = 1/A
print A
A = A/2
print A
print A**0.5
Array Operations
from numpy.random import random as rand
A = np.floor(10*rand((2,3)))
B = np.floor(10*rand(A.shape))
C = np.floor(10*rand(A.shape[::-1]))
D = np.array([[1,1],[2,3]])
print 'A:\n', A, '\n'
print 'B:\n', B, '\n'
print 'C:\n', C, '\n'
print 'D:\n', D, '\n'
print 'A+B:\n', A+B, '\n'
print 'A*B (Element wise):\n', A*B, '\n'
print 'Matrix product:\n', A.dot(C), '\n'
print 'Determinant (D):\n', np.linalg.det(D), '\n'
print 'Transpose (A):\n', A.T, '\n'
print 'Inverse (D):\n',np.linalg.inv(D), '\n'
print 'Pseudo Inverse (A):\n',np.linalg.pinv(A), '\n'
print 'Norm (D):\n',np.linalg.norm(D), '\n'
print 'Sum (A):\n',np.sum(A), '\n' #Similarly, mean, std, min, max, argmin, argmax
print 'Sum of rows (A):\n',np.sum(A,axis=0), '\n'
print 'Sum of columns (A):\n',np.sum(A,axis=1), '\n'
print 'Number of elements in A greater than 5:\n',np.sum(A>5), '\n'
Indexing, Slicing and Iterating
Numpy indexing allows more flexibility in comparison to indexing in Python lists. For example:
x = np.linspace(0,1,10)
y = np.sin(2*np.pi*x)
print 'y:\n',y
print '\n'+'='*10+'\n'
print 'y[2]:\n',y[2]
print 'y[-2]:\n',y[-2]
print 'y[:2]:\n',y[:2]
print 'y[2:]:\n',y[2:]
print 'y[2:-2]:\n',y[2:-2]
print 'y[2:-2:2]:\n',y[2:-2:2]
print 'y[::-1]:\n',y[::-1]
print 'y[[1,3,4]]:\n',y[[1,3,4]] #This is not possible directly with Python lists
#You can also assign values
xnew=x.copy()
xnew[1:3]=[0,1]
print 'xnew:\n',xnew
Below are some examples of indexing with multidimensional arrays:
X=np.random.random((4,5))
print 'X:\n',X
print 'X[0]:\n',X[0]
print 'X[0,:]:\n',X[0,:]
print 'X[:,0]:\n',X[:,0]
print 'X[...,0]:\n',X[...,0] #ellipsis work the same way as the colon
print 'X[2,0]:\n',X[2,0]
print 'X[2][0]:\n',X[2][0]
print 'X[1:3,0]:\n',X[1:3,0]
This demonstrates only basic indexing operations. For more elaborate indexing schemes (such as boolean indexing and ix_(), please see the indexing section in the Numpy tutorial.
You can also iterate over an array as shown below:
for row in X:
print X
for element in x:
print element
Array Shape Manipulation¶
Numpy provides a lot of functions for manipulating the shape of an array. np.reshape is a very useful function in this regard. It returns a new array reshaped to your liking.
print 'x:\n',x
print 'x.rehsape(2,5):\n',x.reshape(2,5) #the reshaped array must have the same number of elements as original array
print 'x.rehsape(2,5):\n',x.reshape(-1,2) #Automatically sets the unspecified dimension
xnew.resize(2,5) will resize the array in place.
print 'xnew:\n',xnew
xnew.resize(2,5) #resize doesn't allow negative index
print 'xnew resized:\n',xnew
print 'X:\n',X
print 'X.ravel():\n',X.ravel()
print 'X.flatten():\n',X.flatten()
Some other useful functions include: vstack, hstack and np.newaxis which adds a dimension to an array.
print 'vstack:\n',np.vstack((x,x))
print 'hstack:\n',np.hstack((x,x))
print 'x[:,np.newaxis]:\n',x[:,np.newaxis]
print 'hstack on x[:,np.newaxis]:\n',np.hstack((x[:,np.newaxis],x[:,np.newaxis]))
Exercise: Using Matrices in Encryption (Hill Cipher)
Hill Cipher is based on simple matrix operations. Here is how it works:
- A given input message (say
s = 'target neutralized') is first coded numerically to giveN(s)as:'116-97-114-103-101-116-32-110-101-117-116-114-97-108-105-122-101-100 A (typically large)
n x nsecret non-singular coding matrix is then used for encryption, for example, $C=\begin{pmatrix} 3 & 3 \\ 2 & 5 \end{pmatrix}$, as follows:Rolling: Create vectors of
ncontiguous characters fromN(s), e.g., $v_0 = \begin{pmatrix} 116 \\ 97 \end{pmatrix}$, $v_1 = \begin{pmatrix} 114 \\ 103 \end{pmatrix}$,.... If the length of the string not a multiple ofn, then empty spaces are appended to the end of the string to make it so. Construct the matrix: $V=\begin{pmatrix}v_0 & v_1 & \ldots\end{pmatrix}$.Compute the matrix: $E = (CV) \bmod 251$ to get the encrypted n-dimensional vectors. 251 is chosen because it is the biggest prime number less than 256.
Unrolling the numbers in $E$ and converting them back to string gives encrypted data.
- Decryption is carried out by using: $d_i = (DE)(\bmod251)$ and then the text is obtained through reverse numerical coding. Here, $D=d \times adj(C)$ with $d$ being a scalar such that $d \times |C| \equiv 1 (\bmod 251)$ and $adj(C)$ is the adjoint of C. $d$ can be found by looping through 0 to 251 until you encounter a value for which $(d \times |C|)\bmod 251)$ equals $1$. If no such value exists then this means that 251 and the determinant of the chosen matrix have common factors and it is not possible to use this matrix for decoding.
Implement a function encryptHill(s,C) which returns the encrypted version of s using Hill Cipher with C as the code matrix. Also implement decryptHill(e,C) which takes a list of encrypted numbers (output of encryptHill and returns the decrypred string.
Sample Solution¶
import numpy as np
import random
import string
def numerify(s,n):
return np.reshape([ord(c) for c in s],(n,-1))
def denumerify(lst):
return ''.join([chr(int(round(i))) for i in lst])
def hillcrypt(s,C,nmax):
n=C.shape[0]
P=numerify(s,n)
E=C.dot(P)%nmax
etext=denumerify(E.flatten())
return etext
def encryptHill(s,C,nmax):
s+=' '*int(n*np.ceil(len(s)/float(n))-len(s))
return hillcrypt(s,C,nmax)
def decryptHill(e,C,nmax):
Ci=moduloInv(C,nmax)
return hillcrypt(e,Ci,nmax)
def moduloInv(C,nmax):
dC=np.linalg.det(C)
di=None
for i in range(nmax):
if (round(dC*i)%nmax)==1:
di=i
break
Ci=dC*np.linalg.inv(C)*di # remember: adj(C) = |C| x inv(C)
return Ci%nmax
n=2
C=np.array([[3,3],[2,5]])
# for random C: np.random.randint(0,2,size=(n,n))
nmax=251
#random text
s= ''.join([random.choice(string.ascii_uppercase) for _ in range(500)])
# _ creates an unnamed variable. Useful when the index isn't going to be used.
etext = encryptHill(s,C,nmax)
dtext = decryptHill(etext,C,nmax)
print s,'\n'+'*'*10
print etext,'\n'+'*'*10
print dtext,'\n'+'*'*10
print 'Decryption Correct:',dtext.strip()==s.strip()
Just for the fun of it, let's encrypt the code for the function 'encryptHill' using itself:
import inspect
ptext = inspect.getsource(encryptHill)
print 'Plain text:\n'+ptext
enc = encryptHill(ptext,C,nmax)
print '\nEncrypted text:\n'+enc
dec = decryptHill(enc,C,nmax)
print '\nDecrypted text:\n'+dec
Extra credit for Home Work:¶
- Try the above solution with
Cbeing a random binary matrix withn = 3. Most of the time you will not be able to decrypt the text. Do you know why? Can you fix it? Try the above solution with
n = 100for a random binary matrix. Most of the time it will not decrypt the message properly. Do you know why? Can you fix it?Hill Cipher is vulnerable to known plain text attacks, i.e., you can find out what the secret matrix was if you are given enough plain text and its encrypted form. Implement
breakHill(s,e)which returns both the numerical coding and the encrypting matrix using the plain textsand its encrypted forme. For ease, you can initially assume that the numerical coding scheme is known. Here is an illustrative example.
(Basic) Plotting¶
Plotting is supported by matplotlib. It's real easy! Given two numpy arrays or python lists x and y, the plt.plot(x,y) function makes a plot of the x vs. y. What it essentially does is takes corresponding points in x and y and, by default, draws a line between them. You can refer to the matplotlib gallery for more elaborate graphing options.
import matplotlib.pyplot as plt
import numpy as np
t = np.linspace(0,1,100)
y = np.sin(2*np.pi*2*t)
plt.plot(t,y) # t and y can be python lists or numpy arrays
plt.xlabel('t')
plt.ylabel('y')
plt.show()
You can also plot a histogram
plt.hist(np.random.random(1000)); plt.show()
Class Exercise: Determinants of random binary matrices¶
Let's see how the determinant of random binary matrices changes in response to their size. We will take a number of random matrices of a certain size n and find the average of the log of the magnitude of the determinant and then plot this average for different values of n.
import numpy as np
from numpy.linalg import det #for determinant
from numpy.random import randint as randi #for random binary matrices
import matplotlib.pyplot as plt
ntrials = range(100) # number of matrices for a single n
nrange = range(2,50,2) # array sizes (n)
dRand=[]
for n in nrange:
d = [np.abs(det(randi(0,2,(n,n)))) for _ in ntrials]
dRand.append(np.mean(d)) # calc
plt.plot(nrange,dRand,'ro-')
plt.xlabel('n')
plt.ylabel('magnitude of determinant')
plt.yscale('log')
plt.grid()
plt.show()
As you can see above, the magnitude of the determinant rises very quickly as the order of the matrix increases. This is why numpy, apart from a providing the det function, numpy also provides a slogdet funtion which returns a tuple containing the sign (1 for positve, -1 for negative and 0 for zero) and the natural log of the determinant of the input matrix.