In [7]:

from numpy.linalg import inv
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from scipy.linalg import eig
from mpl_toolkits.mplot3d import Axes3D

from diffmaps_util import k, diag




In [8]:

X = np.array([.9,1.1,1.2,1.3]).reshape(2,2)
X = np.array([.9,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8]).reshape(3,3)




Out[8]:

array([[ 0.9,  1.1,  1.2],
[ 1.3,  1.4,  1.5],
[ 1.6,  1.7,  1.8]])




In [9]:

%matplotlib inline




In [10]:

fig = plt.figure()
ax.scatter(X[:,0], X[:,1], X[:,2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()






$M = D^{-1}L$



In [11]:

L = k(X, .7)
D = diag(L)
M = inv(D).dot(L)
# Mi,j denotes the transition probability
# from the point xi to the point xj in one time step
print M




[[ 0.55778813  0.34318337  0.09902851]
[ 0.26806012  0.43568764  0.29625225]
[ 0.09557879  0.36606396  0.53835725]]



Equivalente a

$M_{i,j}=\frac{k_\epsilon(x_i, x_j)}{p_\epsilon(x_j)}$

onde

$p_\epsilon(x_j) = \sum_i k_\epsilon(x_i, x_j)$



In [12]:

L/L.sum(axis=1).reshape(-1,1)




Out[12]:

array([[ 0.55778813,  0.34318337,  0.09902851],
[ 0.26806012,  0.43568764,  0.29625225],
[ 0.09557879,  0.36606396,  0.53835725]])



$Ms = D^{1/2}LD^{-1/2}$



In [17]:

Ms = (diag(D,.5)).dot(M).dot(diag(D,-.5))
Ms




Out[17]:

array([[ 0.55778813,  0.30330475,  0.09728836],
[ 0.30330475,  0.43568764,  0.32931333],
[ 0.09728836,  0.32931333,  0.53835725]])



Equivalente a $Ms = \frac{L_{i,j}}{(d(x_i) \times d(x_j))^{1/2}}$



In [16]:

p = L.sum(axis=1)
for i in range(0,3):
a = []
for j in range(0,3):
a.append(L[i,j]/(p[i]*p[j])**.5)
print a




[0.55778812593133109, 0.30330475353984732, 0.097288359660921531]
[0.30330475353984748, 0.43568763730560173, 0.32931333276734809]
[0.097288359660921614, 0.32931333276734809, 0.53835724641066274]





In [15]:

w, v0, v1 = eig(Ms, left=True)
w = w.real
print '%s\n%s' % (w, v0)




[ 1.          0.45188532  0.07994769]
[[-0.54912401 -0.73404914  0.39954309]
[-0.62132302  0.03883155 -0.78259173]
[-0.55894591  0.67798523  0.47740497]]



$P \times \psi_l = \lambda_l \times \psi_l$



In [12]:

Ms.dot(v0)




Out[12]:

array([[-0.54912401, -0.33170603,  0.03194255],
[-0.62132302,  0.01754741, -0.0625664 ],
[-0.55894591,  0.30637157,  0.03816742]])




In [213]:

w * v0




Out[213]:

array([[-0.54912401, -0.33170603,  0.03194255],
[-0.62132302,  0.01754741, -0.0625664 ],
[-0.55894591,  0.30637157,  0.03816742]])




In [40]:

w = w[::-1]




In [71]:

phi = v0.T
phi = phi[::-1]
#print w, '\n', phi
print w, '\n', phi




[ 1.          0.45188532  0.07994769]
[[ 0.39954309 -0.78259173  0.47740497]
[-0.73404914  0.03883155  0.67798523]
[-0.54912401 -0.62132302 -0.55894591]]




In [202]:

psi = []
for i in range(3):
psi.append([])
for j in range(2):
psi[i].append(phi[j+1,i]/M[i])
psi




Out[202]:

[[array([-1.31599994, -2.13894149, -7.41250336]),
array([-0.98446701, -1.60008922, -5.54511042])],
[array([ 0.14486135,  0.08912704,  0.13107597]),
array([-2.31784955, -1.42607448, -2.097277  ])],
[array([ 7.09346943,  1.85209498,  1.25935934]),
array([-5.84801195, -1.52690777, -1.03824349])]]




In [204]:

fig = plt.figure()

for p in psi:
ax.scatter(p[0][0],p[0][1],p[0][2])
ax.scatter(p[1][0],p[1][1],p[1][2])
#print p[0][0], p[0][1], p[0][2]



$$D_t(x,y) = (\sum_{l \geq 1} \lambda_t^{2t} (\psi_l(x) - \psi_l(y))^2)^{1/2}$$


In [75]:

l = w.real[::-1]
print l
psi = v0.T[::-1]
print psi




[ 0.07994769  0.45188532  1.        ]
[[ 0.39954309 -0.78259173  0.47740497]
[-0.73404914  0.03883155  0.67798523]
[-0.54912401 -0.62132302 -0.55894591]]




In [79]:

phi = []
for i in range(3):
phi.append(l[i] * psi[i])
phi




Out[79]:

[array([ 0.03194255, -0.0625664 ,  0.03816742]),
array([-0.33170603,  0.01754741,  0.30637157]),
array([-0.54912401, -0.62132302, -0.55894591])]




In [ ]:

pairwise_distances(phi[1:])**2