Author: bugman Date: Mon Feb 18 18:41:19 2013 New Revision: 18487 URL: http://svn.gna.org/viewcvs/relax?rev=18487&view=rev Log: Spun out the atomic_rmsd() and calc_mean_structure() functions into their own module. They were previously in the generic_fns.structure.superimpose module but are now in the new generic_fns.structure.statistics module. Added: trunk/generic_fns/structure/statistics.py - copied, changed from r18475, trunk/generic_fns/structure/superimpose.py Modified: trunk/generic_fns/structure/superimpose.py Copied: trunk/generic_fns/structure/statistics.py (from r18475, trunk/generic_fns/structure/superimpose.py) URL: http://svn.gna.org/viewcvs/relax/trunk/generic_fns/structure/statistics.py?p2=trunk/generic_fns/structure/statistics.py&p1=trunk/generic_fns/structure/superimpose.py&r1=18475&r2=18487&rev=18487&view=diff ============================================================================== --- trunk/generic_fns/structure/superimpose.py (original) +++ trunk/generic_fns/structure/statistics.py Mon Feb 18 18:41:19 2013 @@ -1,6 +1,6 @@ ############################################################################### # # -# Copyright (C) 2011 Edward d'Auvergne # +# Copyright (C) 2011-2013 Edward d'Auvergne # # # # This file is part of the program relax (http://www.nmr-relax.com). # # # @@ -20,59 +20,11 @@ ############################################################################### # Module docstring. -"""Module for handling all types of structural superimpositions.""" +"""Module for handling all types of structural statistics.""" # Python module imports. -from copy import deepcopy -from math import pi -from numpy import diag, dot, eye, float64, outer, sign, sqrt, transpose, zeros -from numpy.linalg import det, norm, svd - -# relax module import. -from maths_fns.rotation_matrix import R_to_axis_angle - - -class Pivot_finder: - """Class for finding the optimal pivot point for motions between the given models.""" - - def __init__(self, models, coord): - """Set up the class for pivot point optimisation. - - @keyword models: The list of models to use. If set to None, then all models will be used. - @type models: list of int or None - @keyword coord: The array of molecular coordinates. The first dimension corresponds to the model, the second the atom, the third the coordinate. - @type coord: rank-3 numpy array - """ - - # Store the args. - self.models = models - self.coord = coord - - # Store a copy of the coordinates for restoration. - self.orig_coord = deepcopy(coord) - - - def func(self, params): - """Target function for the optimisation of the motional pivot point. - - @param params: The parameter vector from the optimisation algorithm. - @type params: list - @return: The target function value defined as the combined RMSD value. - @rtype: float - """ - - # The fit to mean algorithm. - T, R, pivot = fit_to_mean(models=self.models, coord=self.coord, centroid=params, verbosity=0) - - # The RMSD. - val = atomic_rmsd(self.coord) - - # Restore the coordinates. - self.coord = deepcopy(self.orig_coord) - - # Return the RMSD. - return val - +from numpy import float64, sqrt, zeros +from numpy.linalg import norm def atomic_rmsd(coord): @@ -136,237 +88,3 @@ # Average. mean[i] = mean[i] / M - - -def find_centroid(coords): - """Calculate the centroid of the structural coordinates. - - @keyword coord: The atomic coordinates. - @type coord: numpy rank-2, Nx3 array - @return: The centroid. - @rtype: numpy rank-1, 3D array - """ - - # The sum. - centroid = coords.sum(0) / coords.shape[0] - - # Return. - return centroid - - -def fit_to_first(models=None, coord=None, centroid=None): - """Superimpose a set of structural models using the fit to first algorithm. - - @keyword models: The list of models to superimpose. - @type models: list of int - @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 - @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. - @type centroid: list of float or numpy rank-1, 3D array - @return: The lists of translation vectors, rotation matrices, and rotation pivots. - @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays - """ - - # Print out. - print("\nSuperimposition of structural models %s using the 'fit to first' algorithm." % models) - - # Init (there is no transformation for the first model). - T_list = [zeros(3, float64)] - R_list = [eye(3, dtype=float64)] - pivot_list = [zeros(3, float64)] - - # Loop over the ending models. - for i in range(1, len(models)): - # Calculate the displacements (Kabsch algorithm). - trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='model %s'%models[i], coord_from=coord[i], coord_to=coord[0], centroid=centroid) - - # Store the transforms. - T_list.append(trans_vect) - R_list.append(R) - pivot_list.append(pivot) - - # Return the transform data. - return T_list, R_list, pivot_list - - -def fit_to_mean(models=None, coord=None, centroid=None, verbosity=1): - """Superimpose a set of structural models using the fit to first algorithm. - - @keyword models: The list of models to superimpose. - @type models: list of int - @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 - @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. - @type centroid: list of float or numpy rank-1, 3D array - @keyword verbosity: The amount of information to print out. If 0, nothing will be printed. - @type verbosity: int - @return: The lists of translation vectors, rotation matrices, and rotation pivots. - @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays - """ - - # Print out. - if verbosity: - print("\nSuperimposition of structural models %s using the 'fit to mean' algorithm." % models) - - # Duplicate the coordinates. - orig_coord = deepcopy(coord) - - # Initialise the displacement lists. - T_list = [] - R_list = [] - pivot_list = [] - - # Initialise the mean structure. - N = len(coord[0]) - mean = zeros((N, 3), float64) - - # Iterative fitting to mean. - converged = False - iter = 0 - while not converged: - # Print out. - if verbosity: - print("\nIteration %i of the algorithm." % iter) - print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)", "Rotation (deg)")) - - # Calculate the mean structure. - calc_mean_structure(coord, mean) - - # Fit each model to the mean. - converged = True - for i in range(len(models)): - # Calculate the displacements (Kabsch algorithm). - trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='mean', coord_from=coord[i], coord_to=mean, centroid=centroid, verbosity=0) - - # Table printout. - if verbosity: - print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 / pi * 360.0))) - - # Shift the coordinates. - for j in range(N): - # Translate. - coord[i][j] += trans_vect - - # The pivot to atom vector. - coord[i][j] -= pivot - - # Rotation. - coord[i][j] = dot(R, coord[i][j]) - - # The new position. - coord[i][j] += pivot - - # Convergence test. - if trans_dist > 1e-10 or angle > 1e-10: - converged = False - - # Increment the iteration number. - iter += 1 - - # Perform the fit once from the original coordinates to obtain the full transforms. - for i in range(len(models)): - # Calculate the displacements (Kabsch algorithm). - trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[i], name_to='the mean structure', coord_from=orig_coord[i], coord_to=mean, centroid=centroid, verbosity=0) - - # Store the transforms. - T_list.append(trans_vect) - R_list.append(R) - pivot_list.append(pivot) - - # Return the transform data. - return T_list, R_list, pivot_list - - -def kabsch(name_from=None, name_to=None, coord_from=None, coord_to=None, centroid=None, verbosity=1): - """Calculate the rotational and translational displacements between the two given coordinate sets. - - This uses the Kabsch algorithm (http://en.wikipedia.org/wiki/Kabsch_algorithm). - - - @keyword name_from: The name of the starting structure, used for the printouts. - @type name_from: str - @keyword name_to: The name of the ending structure, used for the printouts. - @type name_to: str - @keyword coord_from: The list of atomic coordinates for the starting structure. - @type coord_from: numpy rank-2, Nx3 array - @keyword coord_to: The list of atomic coordinates for the ending structure. - @type coord_to: numpy rank-2, Nx3 array - @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. - @type centroid: list of float or numpy rank-1, 3D array - @return: The translation vector T, translation distance d, rotation matrix R, rotation axis r, rotation angle theta, and the rotational pivot defined as the centroid of the ending structure. - @rtype: numpy rank-1 3D array, float, numpy rank-2 3D array, numpy rank-1 3D array, float, numpy rank-1 3D array - """ - - # Calculate the centroids. - if centroid != None: - centroid_from = centroid - centroid_to = centroid - else: - centroid_from = find_centroid(coord_from) - centroid_to = find_centroid(coord_to) - - # The translation. - trans_vect = centroid_to - centroid_from - trans_dist = norm(trans_vect) - - # Calculate the rotation. - R = kabsch_rotation(coord_from=coord_from, coord_to=coord_to, centroid_from=centroid_from, centroid_to=centroid_to) - axis, angle = R_to_axis_angle(R) - - # Print out. - if verbosity >= 1: - print("\n\nCalculating the rotational and translational displacements from %s to %s using the Kabsch algorithm.\n" % (name_from, name_to)) - print("Start centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_from[0], centroid_from[1], centroid_from[2])) - print("End centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_to[0], centroid_to[1], centroid_to[2])) - print("Translation vector: [%20.15f, %20.15f, %20.15f]" % (trans_vect[0], trans_vect[1], trans_vect[2])) - print("Translation distance: %.15f" % trans_dist) - print("Rotation matrix:") - print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0, 2])) - print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1, 2])) - print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2, 2])) - print("Rotation axis: [%20.15f, %20.15f, %20.15f]" % (axis[0], axis[1], axis[2])) - print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0)) - - # Return the data. - return trans_vect, trans_dist, R, axis, angle, centroid_to - - -def kabsch_rotation(coord_from=None, coord_to=None, centroid_from=None, centroid_to=None): - """Calculate the rotation via SVD. - - @keyword coord_from: The list of atomic coordinates for the starting structure. - @type coord_from: numpy rank-2, Nx3 array - @keyword coord_to: The list of atomic coordinates for the ending structure. - @type coord_to: numpy rank-2, Nx3 array - @keyword centroid_from: The starting centroid. - @type centroid_from: numpy rank-1, 3D array - @keyword centroid_to: The ending centroid. - @type centroid_to: numpy rank-1, 3D array - @return: The rotation matrix. - @rtype: numpy rank-2, 3D array - """ - - # Initialise the covariance matrix A. - A = zeros((3, 3), float64) - - # Loop over the atoms. - for i in range(coord_from.shape[0]): - # The positions shifted to the origin. - orig_from = coord_from[i] - centroid_from - orig_to = coord_to[i] - centroid_to - - # The outer product. - A += outer(orig_from, orig_to) - - # SVD. - U, S, V = svd(A) - - # The handedness of the covariance matrix. - d = sign(det(A)) - D = diag([1, 1, d]) - - # The rotation. - R = dot(transpose(V), dot(D, transpose(U))) - - # Return the rotation. - return R Modified: trunk/generic_fns/structure/superimpose.py URL: http://svn.gna.org/viewcvs/relax/trunk/generic_fns/structure/superimpose.py?rev=18487&r1=18486&r2=18487&view=diff ============================================================================== --- trunk/generic_fns/structure/superimpose.py (original) +++ trunk/generic_fns/structure/superimpose.py Mon Feb 18 18:41:19 2013 @@ -25,10 +25,11 @@ # Python module imports. from copy import deepcopy from math import pi -from numpy import diag, dot, eye, float64, outer, sign, sqrt, transpose, zeros +from numpy import diag, dot, eye, float64, outer, sign, transpose, zeros from numpy.linalg import det, norm, svd # relax module import. +from generic_fns.structure.statistics import atomic_rmsd, calc_mean_structure from maths_fns.rotation_matrix import R_to_axis_angle @@ -73,69 +74,6 @@ # Return the RMSD. return val - - -def atomic_rmsd(coord): - """Determine the RMSD for the given atomic coordinates. - - This is the per atom RMSD to the mean structure. - - - @keyword coord: The array of molecular coordinates. The first dimension corresponds to the model, the second the atom, the third the coordinate. - @type coord: rank-3 numpy array - """ - - # Init. - M = len(coord) - N = len(coord[0]) - model_rmsd = zeros(M, float64) - mean = zeros((N, 3), float64) - - # Calculate the mean structure. - calc_mean_structure(coord, mean) - - # Loop over the models. - for i in range(M): - # Loop over the atoms. - for j in range(N): - # The vector connecting the mean to model atom. - vect = mean[j] - coord[i][j] - - # The atomic RMSD. - model_rmsd[i] += norm(vect)**2 - - # Normalise, and sqrt. - model_rmsd[i] = sqrt(model_rmsd[i] / N) - - # Return the average RMSD. - return sum(model_rmsd) / M - - -def calc_mean_structure(coord=None, mean=None): - """Average the coordinates. - - @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 - @keyword mean: The data storage for the mean structure. - @type mean: numpy rank-2, Nx3 array - """ - - # The number of atoms. - N = len(coord[0]) - M = len(coord) - - # Clear the mean data structure. - for i in range(N): - mean[i] = [0.0, 0.0, 0.0] - - # Loop over the atoms. - for i in range(N): - # Loop over the models. - for j in range(M): - mean[i] += coord[j][i] - - # Average. - mean[i] = mean[i] / M def find_centroid(coords):