From 97c3105c14cce0327e5e1eef6acff251b6175caf Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sun, 3 Jul 2016 14:18:19 -0700 Subject: [PATCH] Removed anisotropic kernel creation. Figured out that I was unneccesarily exponentiating a diagonalizable matrix. Removed timescaling init, removed tests for anisotropic kernel and timescaling. --- package/MDAnalysis/analysis/diffusionmap.py | 52 ++----------------- .../analysis/test_diffusionmap.py | 27 +--------- 2 files changed, 5 insertions(+), 74 deletions(-) diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index 3ed5c82b106..aea88a728eb 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/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]] diff --git a/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py b/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py index 73149cbad20..2769e88607f 100644 --- a/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py +++ b/testsuite/MDAnalysisTests/analysis/test_diffusionmap.py @@ -37,7 +37,6 @@ 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') @@ -45,35 +44,11 @@ def test_dist_weights(self): 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)