Source Code for Module geometry.mds_algos

  1  from . import (best_similarity_transform, contract, np, assert_allclose, 
  2      project_vectors_onto_sphere, eigh) 
  3  from contracts import check_multiple 
  4  import itertools 
  5   
  6   
  7  @contract(S='array[KxN]', returns='array[NxN](>=0)') 
  8 -def euclidean_distances(S): 
  9      ''' Computes the euclidean distance matrix for the given points. ''' 
 10      K, N = S.shape 
 11      D = np.zeros((N, N)) 
 12      for i in range(N): 
 13          p = S[:, i] 
 14          pp = np.tile(p, (N, 1)).T 
 15          assert pp.shape == (K, N) 
 16          d2 = ((S - pp) * (S - pp)).sum(axis=0) 
 17          d2 = np.maximum(d2, 0) 
 18          D[i, :] = np.sqrt(d2) 
 19      return D 
 20   
 21   
 22 -def double_center(P): 
 23      n = P.shape[0] 
 24   
 25      grand_mean = P.mean() 
 26      row_mean = np.zeros(n) 
 27      col_mean = np.zeros(n) 
 28      for i in range(n): 
 29          row_mean[i] = P[i, :].mean() 
 30          col_mean[i] = P[:, i].mean() 
 31   
 32      R = row_mean.reshape(n, 1).repeat(n, axis=1) 
 33      assert R.shape == (n, n) 
 34      C = col_mean.reshape(1, n).repeat(n, axis=0) 
 35      assert C.shape == (n, n) 
 36   
 37      B2 = -0.5 * (P - R - C + grand_mean) 
 38   
 39      if False: 
 40          B = np.zeros(P.shape) 
 41          for i, j in itertools.product(range(n), range(n)): 
 42              B[i, j] = -0.5 * (P[i, j] - col_mean[j] - row_mean[i] + grand_mean) 
 43          assert_allclose(B2, B) 
 44   
 45      return B2 
 46   
 47   
 48  @contract(C='array[NxN]', ndim='int,>0,K', returns='array[KxN]') 
 49 -def inner_product_embedding_slow(C, ndim): 
 50      U, S, V = np.linalg.svd(C, full_matrices=0) 
 51      check_multiple([('array[NxN]', U), 
 52                      ('array[N]', S), 
 53                      ('array[NxN]', V)]) 
 54      coords = V[:ndim, :] 
 55      for i in range(ndim): 
 56          coords[i, :] = coords[i, :] * np.sqrt(S[i]) 
 57      return coords 
 58   
 59   
 60  @contract(C='array[NxN]', ndim='int,>1,K', returns='array[KxN]') 
 61 -def inner_product_embedding(C, ndim): 
 62      n = C.shape[0] 
 63      eigvals = (n - ndim, n - 1) 
 64      S, V = eigh(C, eigvals=eigvals) 
 65   
 66      assert S[0] <= S[1] # eigh returns in ascending order  
 67   
 68      check_multiple([('K', ndim), 
 69                      ('array[NxK]', V), 
 70                      ('array[K]', S)]) 
 71      coords = V.T 
 72      for i in range(ndim): 
 73          coords[i, :] = coords[i, :] * np.sqrt(S[i]) 
 74      return coords 
 75   
 76   
 77 -def truncated_svd_randomized(M, k): 
 78      ''' Truncated SVD based on randomized projections. ''' 
 79      p = k + 5 # TODO: add parameter 
 80      Y = np.dot(M, np.random.normal(size=(M.shape[1], p))) 
 81      Q, r = np.linalg.qr(Y) #@UnusedVariable 
 82      B = np.dot(Q.T, M) 
 83      Uhat, s, v = np.linalg.svd(B, full_matrices=False) 
 84      U = np.dot(Q, Uhat) 
 85      return U.T[:k].T, s[:k], v[:k] 
 86   
 87   
 88  @contract(C='array[NxN]', ndim='int,>0,K', returns='array[KxN]') 
 89 -def inner_product_embedding_randomized(C, ndim): 
 90      '''  
 91          Best embedding of inner product matrix based on  
 92          randomized projections.  
 93      ''' 
 94      U, S, V = truncated_svd_randomized(C, ndim) #@UnusedVariable. 
 95      check_multiple([('K', ndim), 
 96                      ('array[KxN]', V), 
 97                      ('array[K]', S)]) 
 98      coords = V 
 99      for i in range(ndim): 
100          coords[i, :] = coords[i, :] * np.sqrt(S[i]) 
101      return coords 
102   
103   
104  @contract(D='array[MxM](>=0)', ndim='K,int,>=1', returns='array[KxM]') 
105 -def mds(D, ndim, embed=inner_product_embedding): 
106      diag = D.diagonal() 
107      assert_allclose(diag, 0) 
108      # Find centered cosine matrix 
109      P = D * D 
110      B = double_center(P) 
111      return embed(B, ndim) 
112   
113   
114  @contract(D='array[MxM](>=0)', ndim='K,int,>=1', returns='array[KxM]') 
115 -def mds_randomized(D, ndim): 
116      ''' MDS based on randomized projections. ''' 
117      return mds(D, ndim, embed=inner_product_embedding_randomized) 
118   
119   
120  @contract(C='array[NxN]', ndim='int,>0,K', returns='array[KxN],directions') 
121 -def spherical_mds(C, ndim, embed=inner_product_embedding): 
122      # TODO: check cosines 
123      coords = embed(C, ndim) 
124      proj = project_vectors_onto_sphere(coords) 
125      return proj 
126   
127  # TODO: spherical_mds_randomized 
128   
129  best_embedding_on_sphere = spherical_mds 
130   
131   
132  @contract(references='array[KxN]', distances='array[N](>=0)') 
133 -def place(references, distances): 
134      # TODO: docs 
135      K, N = references.shape 
136   
137      # First do MDS on all data 
138      D = np.zeros((N + 1, N + 1)) 
139      D[:N, :N] = euclidean_distances(references) 
140      D[N, N] = 0 
141      D[N, :N] = distances 
142      D[:N, N] = distances 
143      Sm = mds(D, K) 
144   
145      # Only if perfect   
146      # Dm = euclidean_distances(Sm) 
147      # assert_almost_equal(D[:N, :N], Dm[:N, :N]) 
148      # new in other frame  
149      R, t = best_similarity_transform(Sm[:, :N], references) 
150   
151      Sm_aligned = np.dot(R, Sm) + t 
152      result = Sm_aligned[:, N] 
153   
154      return result 
155