r28103 - /trunk/lib/structure/pca.py -- November 25, 2015 - 18:38

Author: bugman
Date: Wed Nov 25 18:38:04 2015
New Revision: 28103

URL: http://svn.gna.org/viewcvs/relax?rev=28103&view=rev
Log:
Implemented the SVD algorithm for the PCA analysis in the relax library.

This simply calls numpy.linalg.svd().

Modified:
    trunk/lib/structure/pca.py

Modified: trunk/lib/structure/pca.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/lib/structure/pca.py?rev=28103&r1=28102&r2=28103&view=diff
==============================================================================
--- trunk/lib/structure/pca.py  (original)
+++ trunk/lib/structure/pca.py  Wed Nov 25 18:38:04 2015
@@ -24,9 +24,10 @@
 
 # Python module imports.
 from numpy import float64, outer, zeros
-from numpy.linalg import eigh
+from numpy.linalg import eigh, svd
 
 # relax library module imports.
+from lib.errors import RelaxError
 from lib.structure.statistics import calc_mean_structure
 
 
@@ -77,19 +78,28 @@
     # Calculate the covariance matrix for the structures.
     covariance_matrix = calc_covariance_matrix(coord)
 
-    # Perform an eigenvalue decomposition on the covariance matrix.
-    values, vectors = eigh(covariance_matrix)
+    # Perform an eigenvalue decomposition of the covariance matrix.
+    if algorithm == 'eigen':
+        values, vectors = eigh(covariance_matrix)
 
-    # Sort the values and vectors.
-    indices = values.argsort()[::-1]
-    values = values[indices]
-    vectors = vectors[:, indices]
+        # Sort the values and vectors.
+        indices = values.argsort()[::-1]
+        values = values[indices]
+        vectors = vectors[:, indices]
+
+    # Perform a singular value decomposition of the covariance matrix.
+    elif algorithm == 'svd':
+        vectors, values, V = svd(covariance_matrix)
+
+    # Invalid algorithm.
+    else:
+        raise RelaxError("The '%s' algorithm is unknown.  It should be 
either 'eigen' or 'svd'." % algorithm)
 
     # Truncation to the desired number of modes.
     values = values[:num_modes]
     vectors = vectors[:,:num_modes]
 
     # Printout.
-    print("\nThe eigenvalues are:")
+    print("\nThe eigenvalues/singular values are:")
     for i in range(num_modes):
         print("Mode %i:  %10.5f" % (i+1, values[i]))