Changeset 598

Timestamp:

08/25/16 17:44:00 (8 years ago)

Author:

oriona

Message:

MDS algorithm changed to cmdscale. Coordinates export to file added.

File:

• mds-and-trees/mds_plot.py (modified) (4 diffs)

mds-and-trees/mds_plot.py

@@ -16 +16 @@
 
 
-
+#http://www.nervouscomputer.com/hfs/cmdscale-in-python/
+def cmdscale(D):
+    """                                                                                       
+    Classical multidimensional scaling (MDS)                                                  
+                                                                                               
+    Parameters                                                                                
+    ----------                                                                                
+    D : (n, n) array                                                                          
+        Symmetric distance matrix.                                                            
+                                                                                               
+    Returns                                                                                   
+    -------                                                                                   
+    Y : (n, p) array                                                                          
+        Configuration matrix. Each column represents a dimension. Only the                    
+        p dimensions corresponding to positive eigenvalues of B are returned.                 
+        Note that each dimension is only determined up to an overall sign,                    
+        corresponding to a reflection.                                                        
+                                                                                               
+    e : (n,) array                                                                            
+        Eigenvalues of B.                                                                     
+                                                                                               
+    """
+    # Number of points                                                                        
+    n = len(D)
+ 
+    # Centering matrix                                                                        
+    H = np.eye(n) - np.ones((n, n))/n
+ 
+    # YY^T                                                                                    
+    B = -H.dot(D**2).dot(H)/2
+ 
+    # Diagonalize                                                                             
+    evals, evecs = np.linalg.eigh(B)
+ 
+    # Sort by eigenvalue in descending order                                                  
+    idx   = np.argsort(evals)[::-1]
+    evals = evals[idx]
+    evecs = evecs[:,idx]
+ 
+    # Compute the coordinates using positive-eigenvalued components only                      
+    w, = np.where(evals > 0)
+    L  = np.diag(np.sqrt(evals[w]))
+    V  = evecs[:,w]
+    Y  = V.dot(L)
+ 
+    return Y, evals
 
 def rand_jitter(arr):
@@ -31 +76 @@
 
 def compute_mds(distance_matrix, dim):
-        seed = np.random.RandomState(seed=3)
-        mds = manifold.MDS(n_components=int(dim), metric=True, max_iter=3000, eps=1e-9, random_state=seed, dissimilarity="precomputed")
-        embed = mds.fit(distance_matrix).embedding_
+        embed, evals = cmdscale(distance_matrix)
+        variances = compute_variances(embed)
+        embed = np.asarray([embed[:,i] for i in range(dim)]).T
+
+        percent_variances = [sum(variances[:i+1])/sum(variances) for i in range(dim)]
+        for i,pv in enumerate(percent_variances):
+                print(i+1,"dimension:",pv)
+
         return embed
 
@@ -41 +91 @@
         for i in range(len(embed[0])):
                 variances.append(np.var(embed[:,i]))
-        percent_variances = [sum(variances[:i+1])/sum(variances) for i in range(len(variances))]
-        return percent_variances
+        return variances
 
 
@@ -69 +118 @@
         else:
                 plt.savefig(outname+".pdf")
+                np.savetxt(outname+".csv", coordinates, delimiter=";")
 
 
 def main(filename, dimensions=3, outname="", jitter=0, separator='\t'):
+        dimensions = int(dimensions)
         distances = read_file(filename, separator)
         embed = compute_mds(distances, dimensions)
 
-        variances_perc = compute_variances(embed)
-        for i,vc in enumerate(variances_perc):
-                print(i+1,"dimension:",vc)
-
-        dimensions = int(dimensions)
         if dimensions == 1:
                 embed = np.array([np.insert(e, 0, 0, axis=0) for e in embed])