Solutions to the exercises of part 4 / Part 4 の練習問題解答

Functions of a matrix / 行列の関数

import Numeric; N = Numeric
import LinearAlgebra; LA = LinearAlgebra

# Calculate an arbitrary function of a matrix
# 任意の行列関数を計算
def matrixFunction(function, matrix):
    eigenvalues, eigenvectors = LA.eigenvectors(matrix)
    eigenvalues = function(eigenvalues)
    return N.dot(eigenvalues*N.transpose(eigenvectors), eigenvectors)

# Calculate an exponential by Taylor expansion
# 指数関数をテーラー展開で計算
def matrixExponential(matrix):
    sum = N.identity(matrix.shape[0])
    product = sum
    for i in range(1, 50):
	product = N.dot(product, matrix)/i
	sum = sum + product
    return sum

# Create a symmetric matrix.
# 対称行列の生成
m = N.reshape(N.arrayrange(9.), (3, 3))**2/10.
m = (m + N.transpose(m))/2.

# Compare the results
# 結果の比較
print matrixExponential(m)
print matrixFunction(N.exp, m)

Visualization of differences in conformations / 配置の差異の可視化

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector
import string
from Scientific.Visualization import VRML

generic_format = FortranFormat('A6')
atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

# Read the PDB file and make a list of all C-alpha positions
# PDB ファイルを読み込み、全α炭素の位置からなるリストを作成
def readCAlphaPositions(filename):
    c_alpha_positions = []
    for line in TextFile(filename):
	record_type = FortranLine(line, generic_format)[0]
	if record_type == 'ATOM  ' or record_type == 'HETATM':
	    data = FortranLine(line, atom_format)
	    atom_name = string.strip(data[2])
	    position = Vector(data[8:11])
	    if atom_name == 'CA':
		c_alpha_positions.append(position)
    return c_alpha_positions

# Read the two conformations
# 二つの配置の読み込み
conf1 = readCAlphaPositions('conf1.pdb')
conf2 = readCAlphaPositions('conf2.pdb')

# Create the VRML scene and the material for the arrows
# VRML シーンと矢印素材の生成
scene = VRML.Scene([])
arrow_material = VRML.DiffuseMaterial('black')
arrow_radius = 0.2

# Create the spheres and put them into the scene
# 球の生成とシーンへの配置
for i in range(len(conf1)):
    scene.addObject(VRML.Arrow(conf1[i], conf2[i], arrow_radius,
			       material=arrow_material))

# View the scene
# シーンの描画
scene.view()