Fixes to docs, Transform function.

package/MDAnalysis/analysis/diffusionmap.py

@@ -115,13 +115,17 @@
 
 If you use this Dimension Reduction method in a publication, please
 reference:
-..[Ferguson1] Ferguson, A. L.; Panagiotopoulos, A. Z.; Kevrekidis, I. G.;
-Debenedetti, P. G. Chem. Phys. Lett. 2011, 509, 1−11
-..[Lafon1] Coifman, Ronald R., Lafon, Stephane (2006) Diffusion maps.
-Appl. Comput. Harmon. Anal. 21, 5–30.
-..[Clementi1] Rohrdanz, M. A, Zheng, W, Maggioni, M, & Clementi, C. (2013)
+..[Ferguson1] Ferguson, A. L.; Panagiotopoulos, A. Z.; Kevrekidis, I. G.
+Debenedetti,  P. G. Nonlinear dimensionality reduction in molecular
+simulation: The diffusion map approach  Chem. Phys. Lett. 509, 1−11 (2011)
+..[Lafon1] Coifman, Ronald R., Lafon, Stephane Diffusion maps.
+Appl. Comput. Harmon. Anal. 21, 5–30  (2006).
+..[Lafon2] Boaz Nadler, Stéphane Lafon, Ronald R. Coifman, Ioannis G. Kevrekidis.
+Diffusion maps, spectral clustering and reaction coordinates
+of dynamical systems. Appl. Comput. Harmon. Anal. 21 (2006) 113–127
+..[Clementi1] Rohrdanz, M. A, Zheng, W, Maggioni, M, & Clementi, C.
 Determination of reaction coordinates via locally scaled
-diffusion map. Journal of Chemical Physics.
+diffusion map. J. Chem. Phys. 134, 124116 (2011).
 
 .. If you choose the default metric, this module uses the fast QCP algorithm
 [Theobald2005]_ to calculate the root mean square distance (RMSD) between
@@ -243,15 +247,15 @@ class DiffusionMap(object):
     eigenvectors: array
         Eigenvectors of the diffusion map
     diffusion_space : array
-        After calling `transform(num_eigenvectors)` the diffusion map embedding
+        After calling `transform(n_eigenvectors)` the diffusion map embedding
         into the lower dimensional diffusion space will exist here.
 
     Methods
     -------
     run()
         Constructs an anisotropic diffusion kernel and performs eigenvalue
         decomposition on it.
-    transform(num_eigenvectors)
+    transform(n_eigenvectors, time)
         Perform an embedding of a frame into the eigenvectors representing
         the collective coordinates.
     """
@@ -293,7 +297,7 @@ def __init__(self, u, epsilon=1, manifold_density=None,
                           "step size in distance matrix initialization.")
 
         # determines length of diffusion process
-        self._t = timescale
+        self._timescale = timescale
 
         if manifold_density is None:
             # weights do not apply to metric but density of data
@@ -310,10 +314,10 @@ def __init__(self, u, epsilon=1, manifold_density=None,
 
     def run(self):
         # run only if distance matrix not already calculated
-        if not self.dist_matrix._calculated
+        if not self._dist_matrix._calculated:
             self._dist_matrix.run()
-        self._scaled_matrix = self._dist_matrix.dist_matrix / self._epsilon
-
+        self._scaled_matrix = (self._dist_matrix.dist_matrix ** 2 /
+                               self._epsilon)
         # take negative exponent of scaled matrix to create Isotropic kernel
         self._kernel = np.exp(-self._scaled_matrix)
         # define an anistropic diffusion term q
@@ -323,7 +327,8 @@ def run(self):
             q_vector[i] = np.dot(self._kernel[i, :], self._weights_ker)
 
         # Form a new kernel from the anisotropic diffusion term q
-        self._kernel /= np.sqrt(q_vector[:, np.newaxis].dot(q_vector[np.newaxis]))
+        self._kernel /= (np.sqrt(q_vector[:, np.newaxis
+                         ].dot(q_vector[np.newaxis])))
 
         # Weighted Graph Laplacian normalization on this new graph
         d_vector = np.zeros((self._nframes, ))
@@ -337,26 +342,27 @@ def run(self):
         self._kernel /= np.sqrt(d_vector[:, np.newaxis].dot(d_vector[np.newaxis]))
 
         # Apply timescaling
-        if self._t > 1:
-            self._kernel = np.linalg.matrix_power(self._kernel, self._t)
+        if self._timescale > 1:
+            self._kernel = np.linalg.matrix_power(self._kernel,
+                                                  self._timescale)
 
         eigenvals, eigenvectors = np.linalg.eig(self._kernel)
         eg_arg = np.argsort(eigenvals)
         self.eigenvalues = eigenvals[eg_arg[::-1]]
-        self.eigenvalues = self.eigenvalues[1:]
-        self.eigenvectors = eigenvectors[eg_arg[::-1], :]
-        self.eigenvectors = self.eigenvectors[1:]
+        self.eigenvectors = eigenvectors[eg_arg[::-1]]
         self._calculated = True
 
-    def transform(self, n_eigenvectors):
+    def transform(self, n_eigenvectors, time):
         """ Embeds a trajectory via the diffusion map
 
         Parameter
         ---------
         n_eigenvectors : int
             The number of dominant eigenvectors to be used for
             diffusion mapping
-
+        time : float
+            Exponent that eigenvalues are raised to for embedding, for large
+            values, more dominant eigenvectors determine diffusion distance.
         Return
         ------
         diffusion_space : array
@@ -365,5 +371,5 @@ def transform(self, n_eigenvectors):
             between the higher dimensional space and the space spanned by
             the eigenvectors.
         """
-        return (self.eigenvectors.T[:,:n_eigenvectors] *
-                self.eigenvalues[:n_eigenvectors])
+        return (self.eigenvectors[1:n_eigenvectors+1,].T *
+                (self.eigenvalues[1:n_eigenvectors+1]**time))

testsuite/MDAnalysisTests/analysis/test_diffusionmap.py

@@ -17,64 +17,68 @@
 import numpy as np
 import MDAnalysis
 import MDAnalysis.analysis.diffusionmap as diffusionmap
-from numpy.testing import (assert_almost_equal, assert_equal, raises)
+from numpy.testing import (assert_almost_equal, assert_equal,
+                           assert_array_almost_equal,raises)
 
 
 from MDAnalysisTests.datafiles import PDB, XTC
 
 
 class TestDiffusionmap(object):
-    def __init__(self):
+    def setUp(self):
         self.u = MDAnalysis.Universe(PDB, XTC)
         self.dist = diffusionmap.DistanceMatrix(self.u, select='backbone')
         self.dmap = diffusionmap.DiffusionMap(self.dist)
         self.dmap.run()
         self.eigvals = self.dmap.eigenvalues
         self.eigvects = self.dmap.eigenvectors
-        self.weights_ker = np.ones((self.dist.nframes, ))
 
     def test_eg(self):
         # number of frames is trajectory is now 10 vs. 98
-        assert_equal(self.eigvals.shape, (self.dist.nframes-1, ))
+        assert_equal(self.eigvals.shape, (self.dist.nframes, ))
         # makes no sense to test values here, no physical meaning
+        del self.dmap
 
     def test_dist_weights(self):
         backbone = self.u.select_atoms('backbone')
         weights_atoms = np.ones(len(backbone.atoms))
-        dist = diffusionmap.DistanceMatrix(self.u, select='backbone',
+        self.dist = diffusionmap.DistanceMatrix(self.u, select='backbone',
                                            weights=weights_atoms)
-        dist.run()
+        self.dist.run()
+        del self.dist
 
     def test_kernel_weights(self):
-        dist = diffusionmap.DistanceMatrix(self.u, select='backbone')
-        dmap = diffusionmap.DiffusionMap(dist,
+        self.dist = diffusionmap.DistanceMatrix(self.u, select='backbone')
+        self.weights_ker = np.ones((self.dist.nframes, ))
+        self.dmap = diffusionmap.DiffusionMap(self.dist,
                                          manifold_density=self.weights_ker)
-        dmap.run()
-        assert_almost_equal(self.eigvals, dmap.eigenvalues, decimal=5)
-        assert_almost_equal(self.eigvects, dmap.eigenvectors, decimal=6)
+        self.dmap.run()
+        assert_array_almost_equal(self.eigvals, self.dmap.eigenvalues, decimal=5)
+        assert_array_almost_equal(self.eigvects, self.dmap.eigenvectors, decimal=6)
+        del self.dist, self.dmap
 
     @raises(ValueError)
     def test_wrong_kernel_weights(self):
-        dmap = diffusionmap.DiffusionMap(self.dist,
+        self.dmap = diffusionmap.DiffusionMap(self.dist,
                                          manifold_density=np.ones((2,)))
+        del self.dmap
 
     def test_timescaling(self):
-        dmap = diffusionmap.DiffusionMap(self.dist, timescale=2)
-        dmap.run()
-        assert_equal(dmap.eigenvalues.shape, (self.dist.nframes-1, ))
+        self.dmap = diffusionmap.DiffusionMap(self.u, timescale=2)
+        self.dmap.run()
+        assert_equal(self.dmap.eigenvalues.shape, (self.dmap._nframes, ))
+        del self.dmap
 
     def test_different_steps(self):
-        dist = diffusionmap.DistanceMatrix(self.u, select='backbone', step=3)
-        dmap = diffusionmap.DiffusionMap(dist)
-        dmap.run()
+        self.dmap = diffusionmap.DiffusionMap(self.u, select='backbone', step=3)
+        self.dmap.run()
+        del self.dmap, self.dist
 
     def test_transform(self):
-        self.num_eigenvectors = 4
-        self.dmap.transform(self.num_eigenvectors)
-        assert_equal(self.dmap.diffusion_space.shape,
+        self.n_eigenvectors = 4
+        self.dmap = diffusionmap.DiffusionMap(self.u)
+        self.dmap.run()
+        diffusion_space = self.dmap.transform(self.n_eigenvectors,1)
+        assert_equal(diffusion_space.shape,
                      (self.eigvects.shape[0],
-                      self.num_eigenvectors))
-
-    def test_universe_init(self):
-        dmap = diffusionmap.DiffusionMap(self.u)
-        dmap.run()
+                      self.n_eigenvectors))