Diffusion Map Refactor/Implementation

package/MDAnalysis/analysis/diffusionmap.py

@@ -260,8 +260,7 @@ class DiffusionMap(object):
         the collective coordinates.
     """
 
-    def __init__(self, u, epsilon=1, manifold_density=None,
-                 timescale=1, **kwargs):
+    def __init__(self, u, epsilon=1, **kwargs):
         """
         Parameters
         -------------
@@ -273,10 +272,6 @@ def __init__(self, u, epsilon=1, manifold_density=None,
         epsilon : Float
             Specifies the method used for the choice of scale parameter in the
             diffusion map. More information in [1], [2] and [3], Default: 1.
-        manifold_density: list, optional
-            The list has to have the same length as the trajectory.
-            With 'None' the weight of each frame of the trajectory will be the
-            same, Default : None
         timescale: int, optional
             The number of steps in the random walk, large t reflects global
             structure whereas small t indicates local structure, Default: 1
@@ -296,21 +291,6 @@ def __init__(self, u, epsilon=1, manifold_density=None,
                           "be very slow to compute. Consider picking a larger "
                           "step size in distance matrix initialization.")
 
-        # determines length of diffusion process
-        self._timescale = timescale
-
-        if manifold_density is None:
-            # weights do not apply to metric but density of data
-            self._weights_ker = np.ones((self._nframes,))
-        else:
-            if manifold_density.shape[0] != self._nframes:
-                raise ValueError("The weight should have the same length as "
-                                 'the trajectory')
-            else:
-                # density weights are constructed as relative to the mean
-                self._weights_ker = (np.asarray(manifold_density,
-                                     dtype=np.float64) /
-                                     np.mean(manifold_density))
 
     def run(self):
         # run only if distance matrix not already calculated
@@ -320,33 +300,9 @@ def run(self):
                                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
-        q_vector = np.zeros((self._nframes, ))
-        # weights should reflect the density of the points on the manifold
-        for i in range(self._nframes):
-            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])))
-
-        # Weighted Graph Laplacian normalization on this new graph
-        d_vector = np.zeros((self._nframes, ))
-        for i in range(self._nframes):
-            d_vector[i] = np.dot(self._kernel[i, :], self._weights_ker)
-
-        for i in range(self._nframes):
-            self._kernel[i, :] = self._kernel[i, :] * self._weights_ker
-
-        # Define anisotropic transition by dividing kernel by this term
-        self._kernel /= np.sqrt(d_vector[:, np.newaxis].dot(d_vector[np.newaxis]))
-
-        # Apply timescaling
-        if self._timescale > 1:
-            self._kernel = np.linalg.matrix_power(self._kernel,
-                                                  self._timescale)
-
-        eigenvals, eigenvectors = np.linalg.eig(self._kernel)
+        D_inv = np.diag(1 / self._kernel.sum(1))
+        self._diff = np.dot(D_inv, self._kernel)
+        eigenvals, eigenvectors = np.linalg.eig(self._diff)
         eg_arg = np.argsort(eigenvals)
         self.eigenvalues = eigenvals[eg_arg[::-1]]
         self.eigenvectors = eigenvectors[eg_arg[::-1]]

testsuite/MDAnalysisTests/analysis/test_diffusionmap.py

@@ -37,43 +37,18 @@ def test_eg(self):
         # number of frames is trajectory is now 10 vs. 98
         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))
         self.dist = diffusionmap.DistanceMatrix(self.u, select='backbone',
                                            weights=weights_atoms)
         self.dist.run()
-        del self.dist
-
-    def test_kernel_weights(self):
-        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)
-        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):
-        self.dmap = diffusionmap.DiffusionMap(self.dist,
-                                         manifold_density=np.ones((2,)))
-        del self.dmap
-
-    def test_timescaling(self):
-        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):
         self.dmap = diffusionmap.DiffusionMap(self.u, select='backbone', step=3)
         self.dmap.run()
-        del self.dmap, self.dist
-
+
     def test_transform(self):
         self.n_eigenvectors = 4
         self.dmap = diffusionmap.DiffusionMap(self.u)