#!/usr/bin/env python
# coding: utf-8
#
Table of Contents
#
#
#
#
# 数式処理group work-4(Google page rank)
#
#
#
# file:/~/python/doing_math_with_python/group_works_4.ipynb
#
# cc by Shigeto R. Nishitani 2009-2022
#
#
#
# # 目的:Google page rankと線形代数
#
# Googleの躍進の土台となったPageRankは線形代数の計算です.
# 線形代数のコマンドを実際に使いながら,その中身を理解して行くことを
# 目的としています.
# # 解説
# ## Google page rank
#
# Googleのpage rankは非常に単純な仮定
# >「多くの良質なページからリンクされているページはやはり良質なページである」
#
# から成り立っている.博士論文のテーマを探していたスタンフォード大学院コンピュータサイエンス学部時代のラリー・ペイジ(Larry Page)が考案しました.
#
# > 『そのアルゴリズムはペイジの名をとって「ページランク(Page Rank)」と呼ばれたが,
# > 特定のサイトに入るリンクの数と,リンクしたサイトのそれぞれに入るリンクの数の,
# > その両方を考慮に入れる.
# > これは学術論文の引用の度数計算の方法を手本にしており、予想通りに機能した』
# > (『ザ・サーチ グーグルが世界を変えた』ジョン・バッテル著、中谷和男訳,2005年日経BP社)
#
# 詳しい解説はhttp://ja.wikipedia.org/wiki/ページランク
# にある.
# ## ランクの計算手順
#
# つぎのようなリンクが張られたページ群を考える.
#
# ![image.png](attachment:image.png)
#
# まずは,リンクを再現する隣接行列を作る.ページに番号をつけて,その間が結ばれている$i-j$要素を1,そうでない要素を0とする.
# 上の例では,
# ```
# 1 2 3 4 5 6 7
# ---
# 0 1 1 1 1 0 1
# 1 0 0 0 0 0 0
# 1 1 0 0 0 0 0
# 0 1 1 0 1 0 0
# 1 0 1 1 0 1 0
# 1 0 0 0 1 0 0
# 0 0 0 0 1 0 0
# ```
# となる.
#
# 2. この隣接行列を転置する.これはページランクが「どれだけリンクしているか」ではなく,「どれだけリンクされているか」を評価するためである.
# 3. 転置した隣接行列tMのそれぞれの列ベクトルの総和が1となるように規格化して「推移確率行列」をつくる.(ヒント参照)
# 4. 初期ベクトルとして,すべての要素が同じ値で,足して1になるような長さ7の列ベクトルを用意する.
# 5. 初期ベクトルに何度か(例えば10回,あるいは収束するまで)推移確率行列を掛ける.この操作は,行列の最大固有値に属する固有ベクトルを見つけることに相当する.
# 6. 得られたベクトルの各要素が対応するページの得点とみなせ,得点順にランクが高くなる.
# # 演習
# ## Page Rankの手順の確認
#
# うえの手順にしたがって,「推移確率行列(transition probability matrix)」trans_matを作り,初期ベクトル(init_v)に5回ほど作用させて数字の変化を観察せよ.ページランクはどうなるか.
#
# ### ヒント
#
# 手順2,3に於いて,和を取って規格化する代わりに,
# 以下のようにして作成した対角(diagonal)行列VAを,
# 転置した隣接行列tMに右からかければ推移確率行列が得られる.
# ``` python
# > V = [1/5,1,1/2,1/3,1/4,1/2,1]
# > diag(*V)
# ```
#
# ### ヒント
# はじめに
# > init_printing()
#
# をいれておくと,固有値ベクトルとかの表示が綺麗.
# ### 解答例
# In[1]:
from sympy import *
init_printing()
M = Matrix([[0, 1, 1, 1, 1, 0, 1],
[1, 0, 0, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0, 0],
[0, 1, 1, 0, 1, 0, 0],
[1, 0, 1, 1, 0, 1, 0],
[1, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0, 0]])
M
# In[2]:
Mt = M.T
Mt
# In[3]:
#V = [1/5,1,1/2,1/3,1/4,1/2,1] # use float
V = [Rational(1,5),1,Rational(1,2),Rational(1,3),Rational(1,4),Rational(1,2),1]
VA = diag(*V)
VA
# In[4]:
tM = Mt*VA
tM
# In[5]:
v_init = Matrix([1/7,1/7,1/7,1/7,1/7,1/7,1/7])
v_init.shape
# In[6]:
tM* v_init
# In[7]:
tM*tM* tM* tM* tM*v_init
# ページランクは1->5->2->3->4->7->6となる
# ## 初期ベクトルがちがうと...
#
# 初期ベクトルを
# ``` python
# > v_init=Matrix([100,0,0,0,0,0,0])
# ```
#
# として,推移確率行列に右から掛け,それに伴うベクトルの各要素の変化を観察し,前問と比較せよ.その結果から,推移確率行列を掛けることによってどのように状態が推移していくか,漸近していく様子を記述せよ.
#
# ### 解答例
# In[8]:
v_init=Matrix([100,0,0,0,0,0,0])
tM*tM* tM* tM* tM* v_init
# ### コメント
# 要素のそれぞれの値は初期値が違うので,異なっている.しかし,何度かかけた後のページランク(順位)は先ほどの等確率の初期値を使った場合と同じである.初期値がどうであろうと,つまりどこから入ったとしても,いくつかリンクをたどった後は同じサイトに行き着く.ステップを繰り返すと定常状態になることが期待される.これがPage Rankの値の安定性を保証し,ランクの信頼性につながっている.
# ## 固有ベクトルから確認
#
# 推移確率行列の固有値・固有ベクトルをもとめ,最大固有値に対応する固有ベクトルを取り出せ.
# 前問までに得られた結果と比較し,一致していることを確かめよ.
# ただし,固有ベクトルの大きさは任意であるため,norm()で規格化しておくと比較しやすい.
# ### 解答例
# In[9]:
tM.eigenvects(error_when_incomplete=False)
#trans_mat.eigenvects()
# In[10]:
#v = Matrix(trans_mat.eigenvects(error_when_incomplete=False)[0][2]) # use float
v = Matrix(tM.eigenvects(error_when_incomplete=False)[1][2])
v
# In[11]:
v.norm().evalf()
# In[12]:
v_n = v.norm().evalf()
for i in range(0,7):
print(v[i]/v_n)
# ### コメント
# Page Rankが最大固有値の固有ベクトルと一致することが,どのような初期値であろうと最終状態のベクトルが安定することを保証している.この導出は,[la_eigen_vectors.ipynb](https://nbviewer.jupyter.org/github/daddygongon/jupyter_num_calc/blob/master/numerical_calc/la_eigen_vectors.ipynb)の3節に示した.
# ### eigenvectsの注意点
# my mac上ではうまく行くが,大学PC上ではエラーが出る場合がある.versionの違いによるらしい.eigenvectsをerror_when_incomplete=Falseのオプションをつけてやると通る.その場合,初めの推移確率行列の値を単に分数で与えるのでなく,Rationalとして与える必要がある.(挙動の細かいところは不明なので,sympyのversionによるかもしれない.要確認19/5/24,つけないと私のでも動かなくなった20/10/29.)
#
# ### eigenvectsに関するコメント
#
# 最後の固有値を求めるeigenvectsで失敗するかもしれません.
# これはsympyのversionによるようです.1.6.1ではダメです.
# それまでのは行けそうです.
#
# 上の方で,対角ベクトルを与えるところで,
# Rationalを使うとちゃんと求めてくれます.
# ```python
# #V = [1/5,1,1/2,1/3,1/4,1/2,1]
# V = [Rational(1,5),1,Rational(1,2),Rational(1,3),Rational(1,4),Rational(1,2),1]
# ```
#
# これに関する問い合わせと開発者の反応は[Githubのissue](https://github.com/sympy/sympy/issues/20351) の通りです.
# この辺りがOSS開発ソフトで時々おこる問題です.
# そういうものだということを覚悟して,
# いくつかのライブラリや他のソフトで結果をクロス検証することを心がけてください.
#
# <2020-11-02 月>
#
#
#