Author: bugman Date: Wed Nov 25 18:38:06 2015 New Revision: 28104 URL: http://svn.gna.org/viewcvs/relax?rev=28104&view=rev Log: The PCA analysis in the relax library now calculates the per structure projections along the PCs.
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=28104&r1=28103&r2=28104&view=diff ============================================================================== --- trunk/lib/structure/pca.py (original) +++ trunk/lib/structure/pca.py Wed Nov 25 18:38:06 2015 @@ -23,7 +23,7 @@ """Module for performing a principle component analysis (PCA).""" # Python module imports. -from numpy import float64, outer, zeros +from numpy import dot, float64, outer, sqrt, zeros from numpy.linalg import eigh, svd # relax library module imports. @@ -36,14 +36,15 @@ @keyword coord: The list of coordinates of all models to superimpose. The first index is the models, the second is the atomic positions, and the third is the xyz coordinates. @type coord: list of numpy rank-2, Nx3 arrays - @return: The covariance matrix. - @rtype: numpy rank-2, 3N array + @return: The covariance matrix and the deviation matrix. + @rtype: numpy rank-2 3Nx3N array, numpy rank-2 MxNx3 array """ # Init. M = len(coord) N = len(coord[0]) covariance_matrix = zeros((N*3, N*3), float64) + deviations = zeros((M, N, 3), float64) mean_struct = zeros((N, 3), float64) # Calculate the mean structure. @@ -52,16 +53,24 @@ # Loop over the models. for i in range(M): # The deviations from the mean. - deviations = coord[i] - mean_struct + deviations[i] = coord[i] - mean_struct # Sum the covariance element. - covariance_matrix += outer(deviations, deviations) + covariance_matrix += outer(deviations[i], deviations[i]) - # Normalise. + # Average. covariance_matrix /= M # Return the matrix. - return covariance_matrix + return covariance_matrix, deviations + + +def calc_projections(coord=None, num_modes=4): + """Calculate the PCA projections. + + @keyword num_modes: The number of PCA modes to calculate. + @type num_modes: int + """ def pca_analysis(coord=None, algorithm='eigen', num_modes=4): @@ -75,8 +84,12 @@ @type num_modes: int """ + # Init. + M = len(coord) + N = len(coord[0]) + # Calculate the covariance matrix for the structures. - covariance_matrix = calc_covariance_matrix(coord) + covariance_matrix, deviations = calc_covariance_matrix(coord) # Perform an eigenvalue decomposition of the covariance matrix. if algorithm == 'eigen': @@ -95,11 +108,16 @@ 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/singular values are:") for i in range(num_modes): print("Mode %i: %10.5f" % (i+1, values[i])) + + # Calculate the projection for each structure. + proj = zeros((num_modes, M), float64) + for s in range(M): + for mode in range(num_modes): + proj[mode, s] = dot(vectors[:, mode], deviations[s].reshape(N*3)) + + # Truncation to the desired number of modes. + values = values[:num_modes] _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-commits mailing list relax-commits@gna.org To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-commits