#%%
import numpy as np
W= np.array([
  [0.6, 0.1, 0.1, 0.1, 0.1],
  [0.1, 0.6, 0.1, 0.1, 0.1],
  [0.1, 0.1, 0.6, 0.1, 0.1],
  [0.1, 0.1, 0.1, 0.6, 0.1],
  [0.1, 0.1, 0.1, 0.1, 0.6]])

u = np.array([.6, .5, .6, .2, .1])

# var1
M = np.array([
  [-.25, 0, .25, .25, 0],
  [0, -.25, 0, .25, .25],
  [.25, 0, -.25, 0, .25],
  [.25, .25, 0, -.25, 0],
  [0, .25, .25, 0, -.25]])

# other variants (var4)
# M = M * 2

# v_ss = sigmal ((h . e_i) / 1-lambda_i) * e_i

h = np.dot(W,u)

evals, e = np.linalg.eigh(M)

print(e[:,1])

v_ss = np.zeros_like(u).reshape(1,-1)
for i in range(M.shape[0]):
  # print(evecs[i].reshape(1,-1).T)
  coeff = np.dot(h.T, e[:,i]) / (1-evals[i])
  v_ss += coeff * e[:,i]

print(v_ss.T)