py41. Linear Rotation

We can use the Rotation matrix to rotate any line.

Here a line (y=2x+1) is rotated by 30 deg counter clockwise.

Since we don't specify the number of points in linspace, it uses the default num of 50.

# ex41.py
from __future__ import division, print_function
import numpy as np
from numpy import sin, cos, pi
import matplotlib.pyplot as plt
print('Original line: y=5x+1')
print('Rotation matrix: 30 deg ccw (+pi/6)')
ang = pi/6
R = np.array([[cos(ang),-sin(ang)],[sin(ang),cos(ang)]])
print('R (pi/6) =\n',R)
x = np.linspace(-2,2)
xyR = [R.dot([xv,5*xv+1]) for xv in x]
xR = [xv[0] for xv in xyR]
yR = [xv[1] for xv in xyR]
plt.plot(x,5*x+1,'r')
plt.plot(xR,yR,'b')
plt.title("x,y - Red, x',y' - Blue for R(pi/6)")
plt.xlabel("x,x'")
plt.ylabel("y,y'")
plt.xlim((-2,2))
plt.ylim((-2,2))
plt.show()

#    Original line: y=5x+1
#    Rotation matrix: 30 deg ccw (+pi/6)
#    R (pi/6) =
#     [[ 0.8660254 -0.5      ]
#     [ 0.5        0.8660254]]

Output:

py40. Matrix operations in Python

We can use numpy.matrix to create matrices, rather than numpy.array. With numpy.array, the * operation is element-wise multiplication. However, with numpy.matrix, the * operation is real matrix multiplication. With numpy.array, we can always do real matrix multiplication with numpy.dot function.

The function numpy.trace finds the trace of a matrix, which is the sum of the diagonal.

# ex40.py
from __future__ import division, print_function
import numpy as np
A = np.matrix([[1,5,1],[2,-1,6],[1,0,3]])
print('A = \n',A)
B = np.matrix([[2,3,0],[3,-1,7],[4,8,9]])
print('B = \n',B)
print('5*A-10*B+3*A*B =\n',5*A-10*B+3*A*B)
print('trace(A*B)=',np.trace(A*B))
print('trace(B*A)=',np.trace(B*A))

#    A = 
#     [[ 1  5  1]
#     [ 2 -1  6]
#     [ 1  0  3]]
#    B = 
#     [[ 2  3  0]
#     [ 3 -1  7]
#     [ 4  8  9]]
#    5*A-10*B+3*A*B =
#     [[ 48  13 137]
#     [ 55 170 101]
#     [  7   1   6]]
#    trace(A*B)= 103
#    trace(B*A)= 103

py39. Matrix inversion in Python

The function numpy.linalg.inv can be used to find the inverse of matrix.

We find x in the linear equation (Ax=b) using inversion x=inv(A)*b and also with numpy.linalg.solve.

# ex39.py
from __future__ import print_function, division
import numpy as np
A = np.array([[1, 1, 1],[1, -1, 1], [3, 5, -1] ])
print('A =\n',A)
print('shape of A =',A.shape)
inverse = np.linalg.inv(A)
print('inverse of A\n', inverse)
b = np.array([1,0,2])
print('b =\n',b)
print('For Ax=b')
x = np.linalg.solve(A,b)
print('x = ',x)
print('inv(A)*b =',inverse.dot(b))

#A =
# [[ 1  1  1]
# [ 1 -1  1]
# [ 3  5 -1]]
#shape of A = (3, 3)
#inverse of A
# [[-0.5   0.75  0.25]
# [ 0.5  -0.5  -0.  ]
# [ 1.   -0.25 -0.25]]
#b =
# [1 0 2]
#For Ax=b
#x =  [ 0.   0.5  0.5]
#inv(A)*b = [ 0.   0.5  0.5]

py38. Numpy arrays

Numpy arrays contain data of one kind, specified by the dtype attribute, which can be changed with astype(), with the proper string identifier.

We can select values by slicing, giving a list of indices, or a logical expression.

The hstack and vstack functions can join arrays horizontally or vertically. The T attribute can be used to find the transpose of a numpy array. A 2D array will correspond to a matrix, however we also have a matrix class in numpy.

# ex38.py
from __future__ import print_function, division
import numpy as np
A = np.arange(5)
print('A =',A)
print('A.dtype =',A.dtype)
B = A.astype('float')
print('B =',B)
print('B.dtype =',B.dtype)
print('A[2:4] =',A[2:4])
print('A[[2,3]] =',A[[2,3]])
print('A[A==2 | A==3] =',A[(A==2) | (A==3)])
M1 = np.hstack((A,A,A))
print('M1 =',M1)
M2 = np.vstack((A,A,A))
print('M2 =',M2)
M3 = M2.T
print('M3 (transpose M2) =',M3)

#A = [0 1 2 3 4]
#A.dtype = int32
#B = [ 0.  1.  2.  3.  4.]
#B.dtype = float64
#A[2:4] = [2 3]
#A[[2,3]] = [2 3]
#A[A==2 | A==3] = [2 3]
#M1 = [0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]
#M2 = [[0 1 2 3 4]
# [0 1 2 3 4]
# [0 1 2 3 4]]
#M3 (transpose M2) = [[0 0 0]
# [1 1 1]
# [2 2 2]
# [3 3 3]
# [4 4 4]]

aud1. Karplus–Strong string synthesis in Python

String and other sounds can be creating by repeating a random signal while filtering out the high frequency components.

In the note function, we take average of two nearby samples, which is a simple low-pass filter.

Below we have 2 programs. The only function of the first program is to define a dictionary Hz. You can see Dictionary Example. Here I go over the frequency dictionary in more detail.

In the second program, which imports the dictionary, we iterate over the elements of notes list. The duration is the fraction of a second for that particular note. Thus the first two notes are C4 played for 0.25 seconds each.

# Hz.py
Notes = ['C','C#','D','D#','E','F','F#','G','G#','A','A#','B']
Hz ={}
for i in range(88):
    Octave = (i+9)/12
    Pos = (i+9)%12
    S = Notes[Pos]+str(Octave)
    Hz[S] = 27.5*(2**(i/12.0))

# audio1.py

from __future__ import division
import numpy as np
from scipy.io import wavfile
from Hz import Hz
 
SR = 44100
notes = [('C4', 4), ('C4', 4), ('C4', 3), ('D4', 6),
         ('E4', 4), ('E4', 3), ('D4', 6), ('E4', 3),
         ('F4', 6), ('G4', 2), ('C5', 6), ('C5', 6),
         ('C5', 6), ('G4', 6), ('G4', 6), ('G4', 6),
         ('E4', 6), ('E4', 6), ('E4', 6), ('C4', 6),
         ('C4', 6), ('C4', 6), ('G4', 3), ('F4', 6),
         ('E4', 3), ('D4', 6), ('C4', 2)]

max_pos=32767
 
def note(f,num):
    buf=np.random.rand(SR//f)-0.5
    samples=[]
    for i in range(num):
        samples.append(buf[0])
        avg=0.5*(buf[0]+buf[1])
        buf = np.append(buf[1:],avg)
    return np.array([x*max_pos for x in samples])
 
if __name__ == '__main__':
    fname='row.wav'
    out= np.zeros(10)
    for i in range(len(notes)):
        buff = note(Hz[notes[i][0]],SR//notes[i][1])
        out = np.concatenate((out,buff))
    fade_dist = 5000
    x = np.arange(fade_dist-1,-1,-1)
    fade_out = (1-np.exp(-x/fade_dist))/(1-np.exp(-1))
    out[-fade_dist:] *= fade_out 
    wavfile.write(fname,SR,out.astype('int16'))

bpy27. Counting Letters of a biological sequence

We use the createRecords function, that we had used in the last example, by using an appropriate import.

We also define a new function, count, which will count unique letters and return their number of occurrences. This function, in turn, can be used by other modules if they import it.

In the function count, we try to add 1 to number already there, for the appropriate dictionary entry. Should that entry not exist, it is created and set to 1.

# bpy27.py
from __future__ import print_function, division
from bpy26 import createRecords

def count(seq):
    count = {}
    for i in seq:
        try: count[i] = count[i] + 1
        except: count[i] = 1
    return count
    
if __name__ == '__main__':
    records = createRecords('data')
    for record in records:
        print('Entry Name:',record.entry_name)
        print('Sequence counts:')
        print(count(record.sequence))

#Entry Name: IGF1R_HUMAN
#Sequence counts:
#{'A': 72, 'C': 44, 'E': 113, 'D': 62, 'G': 91, 'F': 47,
# 'I': 73, 'H': 23, 'K': 69, 'M': 39, 'L': 116, 'N': 85,
# 'Q': 35, 'P': 81, 'S': 97, 'R': 80, 'T': 69, 'W': 21,
# 'V': 87, 'Y': 63}
#Entry Name: IGF1R_MOUSE
#Sequence counts:
#{'A': 68, 'C': 44, 'E': 114, 'D': 63, 'G': 91, 'F': 48,
# 'I': 73, 'H': 24, 'K': 67, 'M': 40, 'L': 115, 'N': 87,
# 'Q': 36, 'P': 83, 'S': 92, 'R': 82, 'T': 73, 'W': 22,
# 'V': 90, 'Y': 61}

bpy26. Parsing SwissProt files using Biopython

In the last example, we saved the insulin proteins to 'data' subfolder.

We have a function that takes one parameter, the subfolder where *.txt files are stored, each corresponding to a UniProtKB text record.

We could have given our files a different extension during saving, for example .dat or .swiss which would require that term in the list comprehension. We do not have to have a filter term if only UniProtKB records are in the subfolder; in this case we have fils = os.listdir(fol)

Several attributes are printed of the records.

# bpy26.py
from __future__ import print_function, division
import os 
from Bio import SwissProt

def createRecords(fol):
    records = []
    fils = [fil for fil in os.listdir(fol) if fil.endswith('.txt')]
    for fil in fils:
        handle = open(fol + '/' + fil)
        record = SwissProt.read(handle)
        records.append(record)
        handle.close()
    return records

if __name__ == '__main__':
    records = createRecords('data')
    for record in records:
        print('Entry Name:',record.entry_name)
        print('Organism:',record.organism)
        print('Length:',record.sequence_length)
        first_crossref = record.cross_references[0]
        print('First cross ref:')
        for i in first_crossref:
            print('\t',i)
        print()
    
#Entry Name: IGF1R_HUMAN
#Organism: Homo sapiens (Human).
#Length: 1367
#First cross ref:
#         EMBL
#         X04434
#         CAA28030.1
#         -
#         mRNA
#
#Entry Name: IGF1R_MOUSE
#Organism: Mus musculus (Mouse).
#Length: 1373
#First cross ref:
#         EMBL
#         AF056187
#         AAC12782.1
#         -
#         mRNA