test updates and bug fixes

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,30 +160,32 @@ 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) /
                                  np.mean(weights))
 
     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)
 

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)