Index: mds-and-trees/mds_plot.py
===================================================================
--- mds-and-trees/mds_plot.py	(revision 597)
+++ mds-and-trees/mds_plot.py	(revision 598)
@@ -16,5 +16,50 @@
 
 
-
+#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,7 +76,12 @@
 
 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,6 +91,5 @@
 	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,15 +118,12 @@
 	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])
