# Source code for Vector

```#===============================================================================
#     This file is part of TEMPy.
#
#     TEMPy is a software designed to help the user in the manipulation
#     and analyses of macromolecular assemblies using 3D electron microscopy maps.
#
#	  Copyright  2015 Birkbeck College University of London.
#
#				Authors: Maya Topf, Daven Vasishtan, Arun Prasad Pandurangan,
#						Irene Farabella, Agnel-Praveen Joseph, Harpal Sahota
#
#
#
#
#     Farabella, I., Vasishtan, D., Joseph, A.P., Pandurangan, A.P., Sahota, H. & Topf, M. (2015). J. Appl. Cryst. 48.
#
#===============================================================================

from numpy import sqrt,matrix,random
from math import cos,sin,pi,acos,asin, atan2, radians, sqrt as msqrt, pow as mpow

[docs]class Vector:
"""A class representing Cartesian 3-dimensonal vectors."""

def __init__(self, x,y,z):
"""x, y, z = Cartesian co-ordinates of vector."""
self.x = x
self.y = y
self.z = z

def __repr__(self):
return "(%.3f,%.3f,%.3f)" %(self.x, self.y, self.z)

def __getitem__(self, index):
l = [self.x, self.y, self.z]
return l[index]

def __iter__(self):
l = [self.x, self.y, self.z]
return l.__iter__()

[docs]    def copy(self):
"""
Return:
A copy of Vector instance
"""
return Vector(self.x, self.y, self.z)

[docs]    def mod(self):
"""
Return:
The modulus (length) of the vector.
"""
return sqrt(self.x**2 + self.y**2 + self.z**2)

[docs]    def reverse(self):
"""
Flip the direction of a Vector instance.

Return:
A Vector instance
"""
return Vector(-self.x,-self.y,-self.z)

[docs]    def arg(self, vector):
"""Return the argument (angle) between this and another vector.RAD"""
top  = self.dot(vector)
bottom = self.mod()*vector.mod()
#print top, bottom
#print 'top/bottom, top, bottom ', top/bottom, top, bottom
if abs(top-bottom) < 0.00001:
return 0.0
else:
#print 'top/bottom, top, bottom ', top/bottom, top, bottom
#return acos(top/bottom)
return acos(min(max(top/bottom,-1.0),1.0))

[docs]    def times(self, factor):
"""
Multiplies a Vector instance by a scalar factor.

Return:
A Vector instance
"""
return Vector(factor*self.x, factor*self.y, factor*self.z)

[docs]    def dot(self, vector):
"""
Return:
The dot product of this and another vector specified as input parameter.
"""
return vector.x * self.x + vector.y * self.y + vector.z * self.z

[docs]    def cross(self, vector):
"""
Return:
A Vector instance of the cross product of this and another vector specified as input parameter
"""
newX = self.y*vector.z - self.z*vector.y
newY = self.z*vector.x - self.x*vector.z
newZ = self.x*vector.y - self.y*vector.x
return Vector(newX, newY, newZ)

def __sub__(self, vector):
"""Return a Vector instance of the subtraction of a vector from this one."""
newX = self.x - vector.x
newY = self.y - vector.y
newZ = self.z - vector.z
return Vector(newX, newY, newZ)

[docs]    def dist(self, vector):
"""
Return:
The distance between this and another vector specified as input parameter.
"""
return (self-vector).mod()

"""Return a Vector instance of the addition of a vector from this one."""
newX = self.x + vector.x
newY = self.y + vector.y
newZ = self.z + vector.z
return Vector(newX, newY, newZ)

def __mul__(self, prod):
newX = self.x*prod
newY = self.y*prod
newZ = self.z*prod
return Vector(newX, newY, newZ)

def __div__(self, divisor):
newX = self.x/float(divisor)
newY = self.y/float(divisor)
newZ = self.z/float(divisor)
return Vector(newX, newY, newZ)

def __neg__(self):
return Vector(-self.x, -self.y, -self.z)

[docs]    def translate(self, x, y, z):
"""
Translate a Vector instance.

Arguments:
*x, y, z*
distance in Angstroms in respective Cartesian directions to translate vector.

Return:
Vector instance.
"""
newX = self.x + x
newY = self.y + y
newZ = self.z + z
return Vector(newX, newY, newZ)

[docs]    def matrix_transform(self, rot_mat):
"""
Transform the vector using a transformation matrix.

Arguments:
*rot_mat*
a 3x3 Python matrix instance.
Return:
A vector instance
"""
vec_mat = matrix([[self.x],[self.y],[self.z]])
new_pos = rot_mat*vec_mat
x = float(new_pos[0])
y = float(new_pos[1])
z = float(new_pos[2])
return Vector(x,y,z)

[docs]    def to_atom(self):
"""
Create an Atom instance based on Vector instance.

Return:
Atom instance
"""
from ProtRep_Biopy import BioPyAtom
#template = 'ATOM      1  C   NOR A   1      23.161  39.732 -25.038  1.00 10.00             C'
a = BioPyAtom([])
a.x = self.x
a.y = self.y
a.z = self.z
return a

[docs]    def unit(self):
"""
Return:
Vector instance of a unit vector.
"""
mod = self.mod()
if mod==0:
return Vector(0,0,0)
return Vector(self.x/mod, self.y/mod, self.z/mod)

###########################################################################
###########################################################################
###########################################################################
#### def out of the class .
#### better have them separate as these definition are an adaptation of
#### Transformations Python Module from Christoph Gohlke
#### http://www.lfd.uci.edu/~gohlke/
#### http://www.lfd.uci.edu/~gohlke/code/transformations.py.html
###########################################################################
###########################################################################
###########################################################################

[docs]def random_vector(min_v, max_v):
"""
Generate a random vector.
The values for the vector component x, y, and z are randomly sampled between minimum and maximum values specified.

Argument:
*min_v, max_v*
minimum and maximum value
Return:
A Vector instance.
"""

x = random.uniform(min_v, max_v)
y = random.uniform(min_v, max_v)
z = random.uniform(min_v, max_v)
return Vector(x,y,z)

[docs]def axis_angle_to_matrix(x, y, z, turn, rad=False):
"""
Converts the axis angle rotation to a matrix form.

Arguments:
*x, y, z*
axis of rotation (does not need to be normalised).
*turn*
Return:
A 3X3 transformation matrix.
"""
turn = turn*pi/180
c_a = cos(turn)
s_a = sin(turn)
v = Vector(x,y,z).unit()
x = v.x
y = v.y
z = v.z

rot_mat = matrix([[x**2+(1-x**2)*c_a, x*y*(1-c_a)-z*s_a, x*z*(1-c_a)+y*s_a],
[ x*y*(1-c_a)+z*s_a, y**2+(1-y**2)*c_a, y*z*(1-c_a)-x*s_a],
[x*z*(1-c_a)-y*s_a, y*z*(1-c_a)+x*s_a, z**2+(1-z**2)*c_a]])
return rot_mat

"""
Converts an euler rotation to a matrix form.

Arguments:
*x_turn, y_turn, z_turn*
rotation angles around respective axis, in radians if rad=True, else in degrees.
Return:
A 3X3 transformation matrix.
"""
x_turn = x_turn*pi/180
y_turn = y_turn*pi/180
z_turn = z_turn*pi/180
return x_mat*y_mat*z_mat

"""
Converts the axis angle rotation to an Euler form.

Arguments:
*x, y, z*
axis of rotation (does not need to be normalised).
*turn*
Returns:
A 3-tuple (x,y,z) containing the Euler angles. .
"""
turn = turn*pi/180
z_rot = atan2(y*sin(turn)-x*z*(1-cos(turn)), 1-(y**2+z**2)*(1-cos(turn)))
x_rot = asin(x*y*(1-cos(turn))+z*sin(turn))
y_rot = atan2(x*sin(turn)-y*z*(1-cos(turn)), 1-(x**2+z**2)*(1-cos(turn)))
return (x_rot, y_rot, z_rot)

# -- Vector methods for torsion angle geometry -- #

[docs]def torsion(a, b, c):
"""
Find the torsion angle between planes ab and bc.

Arguments:

*a,b,c*
Vector instances.

Returns:
"""
n1 = a.cross(b)
n2 = b.cross(c)
return n1.arg(n2)

[docs]def calcMtrx(arr):
"""
Calculate 3 x 4 transformation matrix from Euler angles and offset.
Arguments:
*arr*
[psi,theta,phi,offsetx,offsety,offsetz].
Returns:
3 x 4 transformation matrix
"""
# Taken from http://mathworld.wolfram.com/EulerAngles.html
# Order psi/theta/phi often reversed

res = [[0 for row in range(4)] for col in range(4)]

res[0][0] = Cpsi * Cphi - Ctheta * Sphi * Spsi
res[0][1] = Cpsi * Sphi + Ctheta * Cphi * Spsi
res[0][2] = Spsi * Stheta
res[1][0] = -Spsi * Cphi - Ctheta * Sphi * Cpsi
res[1][1] = -Spsi * Sphi + Ctheta * Cphi * Cpsi
res[1][2] = Cpsi * Stheta
res[2][0] = Stheta * Sphi
res[2][1] = -Stheta * Cphi
res[2][2] = Ctheta
'''
# y convention
res[0][0] = -Spsi * Sphi + Ctheta * Cphi * Cpsi
res[0][1] = Spsi * Cphi + Ctheta * Sphi * Cpsi
res[0][2] = -Cpsi * Stheta
res[1][0] = -Cpsi * Sphi - Ctheta * Cphi * Spsi
res[1][1] = Cpsi * Cphi - Ctheta * Sphi * Spsi
res[1][2] = Spsi * Stheta
res[2][0] = Stheta * Cphi
res[2][1] = Stheta * Sphi
res[2][2] = Ctheta
'''

res[0][3] = float(arr[3])
res[1][3] = float(arr[4])
res[2][3] = float(arr[5])
return res

def _rotmat_to_axisangle(mat):
"""
Convert rotation matrix to axisangle.
Arguments:
*mat*
Rotation matrix.
Returns:
axis-angle
"""
trace_mat=float(mat[0,0])+float(mat[1,1])+float(mat[2,2])

if trace_mat > 3.0:
print '_rotmat_to_axisangle trace_mat to large:', trace_mat
trace_mat =  3.0
if trace_mat < -1.0:
print '_rotmat_to_axisangle trace_mat to small:', trace_mat
trace_mat = -1.0

theta1=acos((trace_mat-1)/2.0)
if theta1 == 0.0: return (0.0, 1.0,1.0,1.0)
k=1/(2.0*sin(theta1))
# vector axis:
n1=k*(mat[2,1]-mat[1,2])
n2=k*(mat[0,2]-mat[2,0])
n3=k*(mat[1,0]-mat[0,1])
return (theta1,n1,n2,n3)

[docs]def cps(mat_1,mat_2):
"""
Find rotation and translation difference between two transformations.
Arguments:
*mat_1,mat_2*
Transformation matrices.
Returns:
The translation and rotation differences
"""
mat1 = matrix(mat_1)
mat2 = matrix(mat_2)

# mat to euler (angular distance not accurate!)
#t1,p1,s1 = _rotmat_to_euler(mat1)
#t2,p2,s2 = _rotmat_to_euler(mat2)
#ang_magnitude = msqrt(mpow(t2-t1,2)+mpow(p2-p1,2)+mpow(s2-s1,2))
matR1 = matrix([mat_1[0][:-1],mat_1[1][:-1],mat_1[2][:-1]])
matR2 = matrix([mat_2[0][:-1],mat_2[1][:-1],mat_2[2][:-1]])
matR = matR2*matR1.transpose()
ang_magnitude,xV,yV,zV = _rotmat_to_axisangle(matR)
ang_magnitude = ang_magnitude*(180.0/pi)
shift = msqrt(mpow(mat2[0,3]-mat1[0,3],2)+mpow(mat2[1,3]-mat1[1,3],2)+mpow(mat2[2,3]-mat1[2,3],2))
#print (mpow(mat2[0,3]-mat1[0,3],2)+mpow(mat2[1,3]-mat1[1,3],2)+mpow(mat2[2,3]-mat1[2,3],2)), ang_magnitude
#acps_score = (pi/360)*(mpow(mat2[0,3]-mat1[0,3],2)+mpow(mat2[1,3]-mat1[1,3],2)+mpow(mat2[2,3]-mat1[2,3],2))*abs(ang_magnitude)
return shift, ang_magnitude

[docs]def altTorsion(a,b,c):
"""
An alternate and better way to find the torsion angle between planes ab and bc.

Arguments:
*a,b,c*
Vector instances.
Return:
"""
A = a.dot(b.cross(c))*b.mod()
B = (a.cross(b)).dot(b.cross(c))
return atan2(A,B)
[docs]def random_vector2(ul_list):
x = random.uniform(ul_list[0],ul_list[1])
y = random.uniform(ul_list[0],ul_list[1])
z = random.uniform(ul_list[0],ul_list[1])
return Vector(x,y,z)

[docs]def align_2seqs(seq1,seq2):
try:
#NW from biopython
import Bio.pairwise2
aln1 = Bio.pairwise2.align.globalms(''.join(seq1),''.join(seq2),1,-1,-0.5,-0.1)[0]
align_seq1,align_seq2 = list(aln1[0]),list(aln1[1])
except:
print 'ImportError: Bio.pairwise2'
return
return align_seq1,align_seq2
```