d = np.diag(D/ np.diag(D).sum())
d

"""
Explanation: Markov chain stationary distribution
$
\pi(y) = \frac{d(y)}{\sum_{z \in X} d(z)}
$
End of explanation
"""


print d[0] * p[0,1]
print d[1] * p[1,0]

"""
Explanation: The chain is reversible following the detailed balanced condition
$
\pi(x) p(x,y) = \pi(y) p(y,x)
$
End of explanation
"""


e = p0
for i in range(1):
    e = e.dot(M)
print e, sum(e)

w, v = eig(Ms)
w = w.real

print 'phi\n %s' % (v * (d ** .5))
print 'psi\n %s' % (v * (d ** -.5))

w, v = eig(Ms)
w = w.real
print w
print v

c = inv(v).dot(np.eye(3))
#c = inv(v).dot(np.array([1,0,0]))
c

c.dot(p0)

e

c[0]*v[:,0] + c[2]*w[2]*v[:,2] + c[1]*w[1]*v[:,1]

"""
Explanation: 
End of explanation
"""


p = []
for x in range(3):
    _ = []
    for y in range(3):
        _.append(L[x,y]/d[x])
    p.append(_)
p = np.array(p)
print p

"""
Explanation: $
p(x,y) = \frac{k(x,y)}{d(x)}
$
End of explanation
"""


a = []
for i in range(3):
    b = []
    for j in range(3):
        b.append(L[i,j]/(np.sqrt(d[i])*np.sqrt(d[j])))
    a.append(b)
a = np.array(a)
print a

"""
Explanation: $
a(x,y) = \frac{k(x,y)}{\sqrt{\pi(x)} \sqrt{\pi(y)}}
$
End of explanation
"""


a = []
for x in range(3):
    _ = []
    for y in range(3):
        _.append(np.sqrt(e[x])/np.sqrt(e[y])*p[x,y])
    a.append(_)
a = np.array(a)
print a

"""
Explanation: $
a(x,y) = \frac{\sqrt{\pi(x)}}{\sqrt{\pi(y)}} p(x, y)
$
End of explanation
"""


phi = v * (d ** .5)
print phi

"""
Explanation: $
\phi_j = \mathbf{v}_j D^{1/2}
$
End of explanation
"""


print v * v[:,0]

"""
Explanation: $
\phi_l(x) = \Phi_l(x) \times \pi(x)
$
End of explanation
"""


psi = v * (d**-.5)
print psi

"""
Explanation: $
\psi_j = \mathbf{v}_j D^{-1/2} 
$
End of explanation
"""


psi = v / v[:,0]
print psi

"""
Explanation: $
\psi_l(x) = \Phi_l(x) / \pi(x)
$
End of explanation
"""


a = []
for _ in range(3):
    a.append(w[_] * phi[:,0] * phi[:,0])
a = np.array(a)
a

"""
Explanation: $
a(x,y) = \sum_{l \geq 0} \lambda_l \phi_l(x) \phi_l(y)
$
End of explanation
"""


e = p0
for i in range(1):
    e = e.dot(Ms)
print e

p = []
for _ in range(0,3):
    p.append(inv(v).dot(p0)[_] * w[_] * v[:,_])
p = np.array(p)
print p
print p.sum(axis=0)

"""
Explanation: $
p(t,y|x) = \phi_0(y) + \sum_{j \geq 1} a_j(x) \lambda_j^t \phi_j(y)
$
End of explanation
"""


diffmap = w * psi.T
diffmap

"""
Explanation: Diffusion map
$
\Psi_t(x) = (\lambda_1^t \psi_1(x), \lambda_2^t \psi_2(x), ..., \lambda_k^t \psi_k(x), )
$
End of explanation
"""


pairwise_distances(diffmap[:,1:])

"""
Explanation: The diffusion distance is equal to Euclidean distance in the diffusion map space with all (n-1) eigenvectors.
End of explanation
"""


