sel.principleAxes() is not correct

[GoogleCodeExporter]

I've used linalg.svd to produce the principal axes of an alpha helix. These are 
different to that produced by principleAxes().

output from principleAxes() [in red in attached figure]

[-0.67825811 -0.20538766 -0.70553656] [-0.51792553 -0.54748404  0.657278  ] 
[-0.5212668   0.81121954  0.26495998]

output from linalg.avd [in blue in figure]

[ 0.67825812  0.5179255   0.52126682] [-0.20538768 -0.54748404  0.81121951] 
[-0.70553654  0.657278    0.26495996]

I have mapped (10x) these vectors back onto the centre of Geometry (only C 
atoms used) in VMD. A snapshot is attached. As you can clearly see svd does a 
good job of finding the helical axis. 

If you look closely at the coordinates they seem to be transposed i.e. the 
first eigenvector from svd is just the first coordinate from each of the 
principleAxes() and the second eigenvector is the second etc.

I'm not sure where this has happened. I've had a look at the src but am not yet 
good enough in python to be sure and this error seemed like one that a 
developer could fix in a trice (prob whilst slapping a forehead).

my python code:

# define the helix
sel = u.selectAtoms(" ( name CA or name C ) and protein and resid 169:190 ")
selCoords = sel.coordinates()
selCentre = sel.centerOfGeometry()
e1, e2, e3 = sel.principleAxes()
uu, dd, vv = linalg.svd(selCoords - selCentre)
v0h = vv[0]/linalg.norm(vv[0])
v2h = vv[2]/linalg.norm(vv[2])
v1h = vv[1]/linalg.norm(vv[1])
print e1,e2,e3
print v0h,v1h, v2h

excerpt from my MDAnalysis source:

    def momentOfInertia(self):
        # Convert to local coordinates
        recenteredpos = self.coordinates() - self.centerOfMass()
        masses = self.masses()
        values = zip(masses, recenteredpos)
        # Create the inertia tensor
        # m_i = mass of atom i
        # (x_i, y_i, z_i) = pos of atom i
        # Ixx = sum(m_i*(y_i^2+z_i^2)); Iyy = sum(m_i*(x_i^2+z_i^2)); Izz = sum(m_i*(x_i^2+y_i^2))
        # Ixy = Iyx = -1*sum(m_i*x_i*y_i)
        # Ixz = Izx = -1*sum(m_i*x_i*z_i)
        # Iyz = Izy = -1*sum(m_i*y_i*z_i)
        Ixx = reduce(lambda t,a: t+a[0]*(a[1][1]*a[1][1]+a[1][2]*a[1][2]), values, 0.)
        Iyy = reduce(lambda t,a: t+a[0]*(a[1][0]*a[1][0]+a[1][2]*a[1][2]), values, 0.)
        Izz = reduce(lambda t,a: t+a[0]*(a[1][0]*a[1][0]+a[1][1]*a[1][1]), values, 0.)
        Ixy = Iyx = -1*reduce(lambda t,a: t+a[0]*a[1][0]*a[1][1], values, 0.)
        Ixz = Izx = -1*reduce(lambda t,a: t+a[0]*a[1][0]*a[1][2], values, 0.)
        Iyz = Izy = -1*reduce(lambda t,a: t+a[0]*a[1][1]*a[1][2], values, 0.)
        return numpy.array([[Ixx, Ixy, Ixz],[Iyx, Iyy, Iyz],[Izx, Izy, Izz]])
    def principleAxes(self):
        from numpy.linalg import eig
        eigenval, eigenvec = eig(self.momentOfInertia())
        # Sort
        indices = numpy.argsort(eigenval)
        return eigenvec[:,indices] 

Original issue reported on code.google.com by philipwf...@gmail.com on 29 Jun 2010 at 9:01

Attachments:

• test-principleAxes.png

[GoogleCodeExporter]

Indeed, the values are exactly transposed (up to spurious numerical differences 
and the sign of one eigenvector (which is irrelevant)):

  >>> pa = numpy.array([[-0.67825811, -0.20538766, -0.70553656] ,[-0.51792553, -0.54748404,  0.657278  ], [-0.5212668,   0.81121954,  0.26495998]])
  >>> sv = numpy.array([[ 0.67825812  0.5179255   0.52126682] [-0.20538768 -0.54748404  0.81121951] [-0.70553654  0.657278    0.26495996]])
  >>> sv[0] *= - 1    # flip sign of first ev 
  >>> pa - sv.T       # difference
  array([[  9.99999994e-09,  -2.99999999e-08,   2.00000000e-08],
         [  2.00000000e-08,   0.00000000e+00,   2.99999999e-08],
         [ -2.00000000e-08,   0.00000000e+00,   2.00000000e-08]])  

The first eigenvector from principleAxes is pa[:,0]; this implementation 
followed the numpy philosophy of having the coordinates in individual columns 
so that one can pass them separately around.

Is this too confusing? 

Right now I am changing the docs to make this clear at least.

Many thanks for checking!

Original comment by orbeckst on 29 Jun 2010 at 9:55

• Changed state: Started

[GoogleCodeExporter]

I agree that this is confusing.

In r335 I added principalAxes() (correct spelling!) which does what you would 
expect:

  e1,e2,e3 = selection.principalAxes()

The old function principleAxes() with the old behaviour is kept but deprecated.

Original comment by orbeckst on 29 Jun 2010 at 10:20

• Changed state: Fixed

[GoogleCodeExporter]

I've compared the behaviour of the new principalAxes() and my results using svd 
using the python code in my original post and they are identical (to within a 
factor of -1)

Original comment by philipwf...@gmail.com on 29 Jun 2010 at 10:37

[GoogleCodeExporter]

That's good – many thanks for checking.

Original comment by orbeckst on 29 Jun 2010 at 10:54

• Changed state: Verified