diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index 5a923a5b001..b2203f2afa2 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -100,6 +100,7 @@ logger = logging.getLogger("MDAnalysis.analysis.diffusionmap") + class DiffusionMap(AnalysisBase): """Non-linear dimension reduction method @@ -137,8 +138,8 @@ def __init__(self, u, select='all', epsilon='average', k=10, weights=None, different frames. Water should be excluded. epsilon : float, optional Specifies the epsilon used for the diffusion map. More information - in [1] and [2]. With 'average' the average of the RMSD to the - k-nearest-neighbor will be used. + in [1] and [2]. With 'average' the average of the RMSD to the k-th + greatest value is used. (#TODO, check if actually agrees with theory) k : int, optional specifies the k for the epsilon choice if average epsilon is used. @@ -159,22 +160,25 @@ def __init__(self, u, select='all', epsilon='average', k=10, weights=None, self._u = u self.atoms = u.select_atoms(select) self._natoms = self.atoms.n_atoms - self._k = k + # modulus to prevent index out of bounds exception traj = u.trajectory self._epsilon = epsilon if metric is not None: self._metric = metric else: self._metric = qcp.CalcRMSDRotationalMatrix - + # remember that this must be called before referencing self.nframes self._setup_frames(traj, start, stop, step) + # modulus to prevent array index out of bounds exception + self._k = k % self.nframes if weights is None: # weights do not apply to metric self._weights = np.ones((self.nframes,)) else: if weights.shape[0] != self.nframes: - raise ValueError("The weight should have the same length as the trajectroy") + raise ValueError("The weight should have the same length as the" + 'trajectory') else: # weights are constructed as relative to the mean self._weights = (np.asarray(weights, dtype=np.float64) / @@ -182,7 +186,6 @@ def __init__(self, u, select='all', epsilon='average', k=10, weights=None, def _prepare(self): self.diffusion_matrix = np.zeros((self.nframes, self.nframes)) - if self._epsilon == 'average': self._epsilon = np.zeros((self.nframes, ), ) self._type_epsilon = 'average' @@ -202,11 +205,13 @@ def _single_frame(self): # for significant decimal places to sparsify matrix for j in range(self.nframes-1, self._ts.frame-1, -1): j_ref = self._u.trajectory[j].positions-self.atoms.center_of_mass() - self.diffusion_matrix[traj_index, j] = self._metric(i_ref.T.astype(np.float64), - j_ref.T.astype(np.float64), - self._natoms, self._rot, - weights=None) - + ij_result = self._metric(i_ref.T.astype(np.float64), + j_ref.T.astype(np.float64), self._natoms, + self._rot, weights=None) + if(ij_result < 1E-05): + self.diffusion_matrix[traj_index, j] = 0 + else: + self.diffusion_matrix[traj_index, j] = ij_result def _conclude(self): self.diffusion_matrix = (self.diffusion_matrix + self.diffusion_matrix.T - @@ -243,10 +248,12 @@ def _conclude(self): for i in range(self.nframes): d_vector[i] = np.dot(self._kernel2[i, :], self._weights) + logger.info('printing d_vector matrix: {0}'.format(d_vector)) for i in range(self.nframes): self._kernel2[i, :] = self._kernel2[i, :] * self._weights self._kernel2 /= np.sqrt(d_vector[:, np.newaxis].dot(d_vector[np.newaxis])) + logger.info('printing final kernel: {0}'.format(self._kernel2)) # eigenvalues and eigenvector are the collective coordinates eigenvals, eigenvectors = np.linalg.eig(self._kernel2) diff --git a/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py b/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py index 7bb98cb1a40..821405eee0f 100644 --- a/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py +++ b/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py @@ -14,6 +14,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # from __future__ import print_function +import numpy as np import MDAnalysis import MDAnalysis.analysis.diffusionmap as diffusionmap from numpy.testing import (assert_almost_equal, assert_equal) @@ -33,22 +34,26 @@ def __init__(self): #assert_equal(self.ev.shape, (98,98)) #assert_almost_equal(self.ev[0,0], .095836037343022831) #faster - u = MDAnalysis.Universe(PDB, XTC) - dmap = diffusionmap.DiffusionMap(u, select='backbone', k = 5) - dmap.run() - self.eigvals = dmap.eigenvalues - self.eigvects = dmap.eigenvectors - + self.u = MDAnalysis.Universe(PDB, XTC) + self.dmap = diffusionmap.DiffusionMap(self.u, select='backbone', k=5) + self.dmap.run() + self.eigvals = self.dmap.eigenvalues + self.eigvects = self.dmap.eigenvectors + self.weights = np.ones((self.dmap.nframes, )) def test_eg(self): - assert_equal(self.eigvals.shape, (10, )) + #number of frames is trajectory is now 10? + assert_equal(self.eigvals.shape, (self.dmap.nframes, )) assert_almost_equal(self.eigvals[0], 1.0, decimal=5) - assert_almost_equal(self.eigvals[-1], 0.0142, decimal = 3) - + assert_almost_equal(self.eigvals[-1], 0.0142, decimal=3) def test_ev(self): - assert_equal(self.eigvects.shape, (10, 10)) - assert_almost_equal(self.eigvects[0, 0], -0.3019, decimal = 2) - - - + assert_equal(self.eigvects.shape, (self.dmap.nframes, self.dmap.nframes)) + assert_almost_equal(self.eigvects[0, 0], -0.3019, decimal=2) + + def test_weights(self): + dmap2 = diffusionmap.DiffusionMap(self.u, select='backbone', + weights=self.weights, k=5) + dmap2.run() + assert_almost_equal(self.eigvals, dmap2.eigenvalues, decimal=5) + assert_almost_equal(self.eigvects, dmap2.eigenvectors, decimal=6)