ネットワークのクラスタリング - GraphRicciCurvature - 備忘録

#概要
コミュニティのクラスタリングがいずれ仕事にも使えそうなので、お試ししてみた備忘録を残す。
いろいろ調べててGraphRicciCurvatureを見つけた。ユークリッド空間ではなく位相空間で引っ張ったり伸ばしたりするやつらしい。なんか面白そう。
でもGraphRicciCurvatureで検索してもヒットしなかったので、Qitta初か？
それほど新しいアルゴリズムではないようだが、日本語情報はいつにも増して少ない。

- 実施期間: 2021年8月
- 環境：Ubuntu20.04 LTS
- パケージ：GraphRicciCurvature, networkxなど

##1. 準備
パケージをインストールする。[前回のpost](https://qiita.com/Ihmon/items/9da571b83973a513cb6c)が動作するパケージはインストール済みとする。
[オフィシャルの通り](https://github.com/saibalmars/GraphRicciCurvature)にインストールするだけだが、依存パケージのコンパイルの必要があるので諸々必要。

```terminal
# 入っていなければ。
sudo apt install g++
sudo apt install cmake
sudo apt install clang

# 各パケージのオフィシャルの通り
pip3 install networkit
pip install POT
pip3 install GraphRicciCurvature
```

今回使用するパケージをすべてimportする。
今回もお試しのデータセットはkarateを使用する。
基本的には[オフィシャルチュートリアル](https://graphriccicurvature.readthedocs.io/en/v0.5.1/tutorial.html)に沿ってかいつまんで記載する。
チュートリアルはColabでも動作するように書かれているが、何故かcommunityが見つからないとエラーが出る。Colab上でGithubのソースからpython-louvainをインストールすれば良いような気もするがGPUは必要ないので試していない。

```python
import networkx as nx
import numpy as np
import math

# matplotlib setting
%matplotlib inline
import matplotlib.pyplot as plt

# to print logs in jupyter notebook
import logging
logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.ERROR)

# load GraphRicciCuravture package
from GraphRicciCurvature.OllivierRicci import OllivierRicci
from GraphRicciCurvature.FormanRicci import FormanRicci

# load python-louvain for modularity computation
import community as community_louvain

# for ARI computation
from sklearn import preprocessing, metrics

import GraphRicciCurvature
print(GraphRicciCurvature.__version__)
```

*0.5.3*

helper関数は長いので最後に記載する。



##2. GraphRicciCurvatureでクラスタリング
トポロジーの知識が必要なようだが、そんなものは専攻したことがないので[理屈](https://ja.wikipedia.org/wiki/%E3%83%AA%E3%83%83%E3%83%81%E3%83%86%E3%83%B3%E3%82%BD%E3%83%AB)は理解できない。
パケージの解説によると、
「edgeのRicci曲率がグラフの構造に重要な役割を果たしていることが観察されています。正の曲率を持つedgeはクラスタ内のedgeを表し、負の曲率を持つedgeはクラスタ内の架け橋となります。また、負の曲率を持つedgeはグラフの連結性に大きく関係しており、連結されたグラフから負の曲率を持つedge取り除くと、グラフはすぐに切断されてしまいます。
Ricci flowは、グラフのエッジのRicci曲率を均一化するプロセスである。Ricci flowは、与えられたグラフに対して、各edgeに"Ricci flow metric"をedge weightとして与え、このedge weightより小さいければグラフのRicci曲率がどこでもほぼ等しくなるようにするものです。この"Ricci flow metric"がコミュニティを検出できます。
Ricci曲率とRicci flow metricの両方とも、グラフ分類のためのグラフフィンガープリントとして機能します。グラフが異なれば、edgeのRicci曲率分布もRicci flow metricも異なります。」

あわせて[ココのレクチャ](https://www.youtube.com/watch?v=PwRl5W-whTs&t=6s)によると、metric tensor: gは球体のどの部分を伸長・収縮させるかを示す。その曲率をRicci curveture: Rと呼び、一般式は[ココ](https://en.wikipedia.org/wiki/Ricci_curvature)に定義あり。
従い、gがわかればRもわかる。
これをRicci flowに適用する。Ricci flowではgをiterativeに変化させ球体を球に近づける。
このとき球体の凹んだ部分(R<0)は膨らみ、凸った部分(R>0)はへこむ。つまり凹んだ部分では膨らますためgは増加し、凸った部分ではgは減少する。
R = -dg/dt
従いほっておくとR>0なら縮退し消滅する。
縮退の過程で球体の一部が限りなくゼロに近づく（ネック）が放置するとゼロに縮退するのでこのネック部分をマニュアルで削除し、その部分に球体をかぶせる(Surgery Theory)。この後の解説がどうしても納得できなかった。

まぁ、このパケージではこのRがiterationを通じて一定になるようにしているのでグラフは消滅しないし、更新されるedge weightを閾値としてクラスタ間を切断(surgery)することもできる。

![Screenshot from 2021-08-29 11-55-26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/3bc2b9e3-a4c6-c667-6b4f-9ea69ac2de36.png)


アプローチが異なる2つの手法が用意されている。

**Ollivier-Ricci curvature:** 
　最適輸送理論に基づく曲率。各edgeの曲率をミクロに見ることができますが、計算量が多くなります。
**Forman-Ricci curvature:**
　CWコンプレックスに基づく曲率です。各エッジの曲率をマクロ的に見ることができ、計算速度も非常に軽い。

使い方は同じなので前者のコードを残す。
全APIは[ココ](https://graphriccicurvature.readthedocs.io/en/latest/GraphRicciCurvature.html)を参照

###2.1 Ollivier-Ricci curvature

```python
G = nx.karate_club_graph()
draw_graph(G_rf)

orc = OllivierRicci(G, alpha=0.5, verbose="INFO")

orc.compute_ricci_curvature()
G_orc = orc.G.copy()  # save an intermediate result

show_results(G_orc)
```

*Karate Club Graph, first 5 edges: 
Ricci curvature of edge (0,1) is 0.111111
Ricci curvature of edge (0,2) is -0.143750
Ricci curvature of edge (0,3) is 0.041667
Ricci curvature of edge (0,4) is -0.114583
Ricci curvature of edge (0,5) is -0.281250*

![Screenshot from 2021-08-29 12-07-38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/e9e14cf4-d508-a1ba-ace6-3a2d7d10beba.png)
通常、Ricci curvature(以降R)は[-1, 1]の間に入る。
weightは重要なパラメータだが、元のに設定されていなければdefaultで1となる。

###2.2 Ricci flow
Rを一定に、かつ曲率をなめらかにしながらedge weightを更新する。Rが大きくマイナス/プラスになるとそのedgeを伸ばす/縮めることとなる。
経験的には20~50回iterationすると、Rはなめらかになるらしい。
とりあえず一度iterationを回すとどうなるか見てみる。

```python
# Start a Ricci flow with Lin-Yau's probability distribution setting with 4 process.
orf = OllivierRicci(G, alpha=0.5, base=1, exp_power=0, proc=4, verbose="INFO")

# Do Ricci flow for 2 iterations
orf.compute_ricci_flow(iterations=2)

orf.set_verbose("ERROR") # mute logs
orf.compute_ricci_flow(iterations=50)  # 上の2回(iterations=2)に続き、行なうことになる。
G_rf = orf.G.copy()

show_results(G_rf)
```

*Karate Club Graph, first 5 edges: 
Ricci curvature of edge (0,1) is -0.002740
Ricci curvature of edge (0,2) is -0.002740
Ricci curvature of edge (0,3) is -0.002740
Ricci curvature of edge (0,4) is -0.002740
Ricci curvature of edge (0,5) is -0.002740*

![Screenshot from 2021-08-29 14-20-06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/215d545a-a816-c297-70ac-2a7c4de57c9f.png)

Rが一定値に収束し、全edge weightが調整されていることがわかる。


```python
draw_graph(G_rf)
```

![Screenshot from 2021-08-29 14-23-36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/b2a534d2-c05b-bac6-51ba-36f1192198a0.png)

surgery実行していないので全edgeはそのまま。色は元の"club"属性を示す（以降同じ）。

###2.3 Ricci flow for Community detection
クラスタリングにおいては、weight閾値を見つけそれよりも大きなweightを持つedgeを取り除く。

例えば,weight=1.5と1.0でカットした場合を見てみる。

```python
draw_graph(my_surgery(G_rf, cut=1.5))
```
![Screenshot from 2021-08-29 14-38-36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/d678c0a6-71e3-fd1e-3ba3-017b7bdbd764.png)

左上のnode=[4,5,6,10,15]は[前回のLouvain](https://qiita.com/Ihmon/items/9da571b83973a513cb6c#3-louvain%E3%81%A7%E3%82%AF%E3%83%A9%E3%82%B9%E3%82%BF%E3%83%AA%E3%83%B3%E3%82%B0)のデンドロが最初に見つけた小グループと同じである。全く違うアルゴリズムだが見つけるものは見つけてくれる。

```python
draw_graph(my_surgery(G_rf, cut=1.0))
```
![Screenshot from 2021-08-29 14-39-17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/adaa95a0-c08d-2cb9-f226-efef84e1f06e.png)

cut=1.5や1.0はaccuracyを見ながら最適値に調整する。次項で説明する。

###2.4 Optimal transportation and Communities
別のインスタンス(orf2)を作成して、それでやってみる。

```python
draw_graph(my_surgery(G_rf, cut=1.0))
orf2 = OllivierRicci(G, alpha=0.5, base=math.e, exp_power=1, verbose="ERROR")
orf2.compute_ricci_flow(iterations=50)
G_rf2 = orf2.G.copy()

show_results(G_rf2)
```
*Karate Club Graph, first 5 edges: 
Ricci curvature of edge (0,1) is -0.019349
Ricci curvature of edge (0,2) is -0.001516
Ricci curvature of edge (0,3) is -0.010322
Ricci curvature of edge (0,4) is -0.001809
Ricci curvature of edge (0,5) is -0.002310*
![Screenshot from 2021-08-29 14-47-12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/e3d48f9c-adab-7d17-1ecd-38678930a7df.png)

ヒストグラムが前項と異なるのは引数が**base=math.e, exp_power=1**に変わっているから。これと**alpha**を調整する必要がありそう。チュートリアルによると、もしくは**base=math.e, exp_power=2**が理論的に最適（熱拡散に近い）になるとのこと。なお、exp_power=0にすると元論文となる。
**alpha**は、確率分布の元nodeと隣接nodeの分割割合になる。[0, 1.0]

![Screenshot from 2021-08-29 15-23-46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/29df4ed1-b9b4-3fc9-d6f1-033c62cf1ed0.png)


このときのaccuracyを確認する。

```python
check_accuracy(G_rf2, clustering_label="club",plot_cut=True)
```

![Screenshot from 2021-08-29 14-50-50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/4d5a5619-e5d5-88ee-909f-82d15c111805.png)
チャート中のGood cutoffは多分weight=0から増加させながらaccが最初に落ちたときを拾っているぽい。
全体を見るとaccが最高なのは3.0 ~ 3.88の間なので3.88でcutしてみる。

```python
draw_graph(my_surgery(G_rf2, cut=3.88))
```

![Screenshot from 2021-08-29 15-05-08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/6d1059ac-4b02-59da-2cdd-0acccb70644e.png)


うまく分離できている。node=9とnode=8を取り違えているのは前回のLouvainと同じだが、大きく異なるのはRによりedgeを任意に取り除くことができる点で、こっちのほうがよい。
ただ、3.0~3.88を求めるロジックやalphaの調整ロジックは考えないといけない。
一意に決まらないようであれば使うのは難しいか？


ちなみにcutoffを2.0にすると例の小クラスタが分離される。

```python
draw_graph(my_surgery(G_rf2, cut=2.0))
```

![Screenshot from 2021-08-29 15-42-49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205257/c182a053-e386-4c6d-33a2-79c3f9e6b961.png)

答えはないのでクラスタリングを何に使うかによって調整しなければならない。

### 2.5 Helper関数
helper関数をオフィシャルから転記しておく。

```python
def check_accuracy(G_origin, weight="weight", clustering_label="value", plot_cut=True):
    """To check the clustering quality while cut the edges with weight using different threshold

    Parameters
    ----------
    G_origin : NetworkX graph
        A graph with ``weight`` as Ricci flow metric to cut.
    weight: float
        The edge weight used as Ricci flow metric. (Default value = "weight")
    clustering_label : str
        Node attribute name for ground truth.
    plot_cut: bool
        To plot the good guessed cut or not.

    """
    G = G_origin.copy()
    modularity, ari = [], []
    maxw = max(nx.get_edge_attributes(G, weight).values())
    cutoff_range = np.arange(maxw, 1, -0.025)
    print(maxw)

    for cutoff in cutoff_range:
        edge_trim_list = []
        for n1, n2 in G.edges():
            if G[n1][n2][weight] > cutoff:
                edge_trim_list.append((n1, n2))
        G.remove_edges_from(edge_trim_list)

        # Get connected component after cut as clustering
        clustering = {c: idx for idx, comp in enumerate(nx.connected_components(G)) for c in comp}

        # Compute modularity and ari
        modularity.append(community_louvain.modularity(clustering, G, weight))
        ari.append(ARI(G, clustering, clustering_label=clustering_label))

    plt.xlim(maxw, 0)
    plt.xlabel("Edge weight cutoff")
    plt.plot(cutoff_range, modularity, alpha=0.8)
    plt.plot(cutoff_range, ari, alpha=0.8)

    if plot_cut:
        good_cut = -1
        mod_last = modularity[-1]
        drop_threshold = 0.01  # at least drop this much to considered as a drop for good_cut

        # check drop from 1 -> maxw
        for i in range(len(modularity) - 1, 0, -1):
            mod_now = modularity[i]
            if mod_last > mod_now > 1e-4 and abs(mod_last - mod_now) / mod_last > drop_threshold:
                if good_cut != -1:
                    print("Other cut:%f, diff:%f, mod_now:%f, mod_last:%f, ari:%f" % (
                        cutoff_range[i + 1], mod_last - mod_now, mod_now, mod_last, ari[i + 1]))
                else:
                    good_cut = cutoff_range[i + 1]
                    print("*Good Cut:%f, diff:%f, mod_now:%f, mod_last:%f, ari:%f" % (
                        good_cut, mod_last - mod_now, mod_now, mod_last, ari[i + 1]))
            mod_last = mod_now

        plt.axvline(x=good_cut, color="red")
        plt.legend(['Modularity', 'Adjust Rand Index', 'Good cut'])
    else:
        plt.legend(['Modularity', 'Adjust Rand Index'])


def show_results(G, curvature="ricciCurvature"):

    # Print the first five results
    print("Karate Club Graph, first 5 edges: ")
    for n1,n2 in list(G.edges())[:5]:
        print("Ricci curvature of edge (%s,%s) is %f" % (n1 ,n2, G[n1][n2][curvature]))

    # Plot the histogram of Ricci curvatures
    plt.subplot(2, 1, 1)
    ricci_curvtures = nx.get_edge_attributes(G, curvature).values()
    plt.hist(ricci_curvtures,bins=20)
    plt.xlabel('Ricci curvature')
    plt.title("Histogram of Ricci Curvatures (Karate Club)")

    # Plot the histogram of edge weights
    plt.subplot(2, 1, 2)
    weights = nx.get_edge_attributes(G, "weight").values()
    plt.hist(weights,bins=20)
    plt.xlabel('Edge weight')
    plt.title("Histogram of Edge weights (Karate Club)")

    plt.tight_layout()


def ARI(G, clustering, clustering_label="club"):
    """
    Computer the Adjust Rand Index (clustering accuracy) of "clustering" with "clustering_label" as ground truth.

    Parameters
    ----------
    G : NetworkX graph
        A given NetworkX graph with node attribute "clustering_label" as ground truth.
    clustering : dict or list or list of set
        Predicted community clustering.
    clustering_label : str
        Node attribute name for ground truth.

    Returns
    -------
    ari : float
        Adjust Rand Index for predicted community.
    """

    complex_list = nx.get_node_attributes(G, clustering_label)

    le = preprocessing.LabelEncoder()
    y_true = le.fit_transform(list(complex_list.values()))

    if isinstance(clustering, dict):
        # python-louvain partition format
        y_pred = np.array([clustering[v] for v in complex_list.keys()])
    elif isinstance(clustering[0], set):
        # networkx partition format
        predict_dict = {c: idx for idx, comp in enumerate(clustering) for c in comp}
        y_pred = np.array([predict_dict[v] for v in complex_list.keys()])
    elif isinstance(clustering, list):
        # sklearn partition format
        y_pred = clustering
    else:
        return -1

    return metrics.adjusted_rand_score(y_true, y_pred)


def my_surgery(G_origin: nx.Graph(), weight="weight", cut=0):
    """A simple surgery function that remove the edges with weight above a threshold

    Parameters
    ----------
    G_origin : NetworkX graph
        A graph with ``weight`` as Ricci flow metric to cut.
    weight: str
        The edge weight used as Ricci flow metric. (Default value = "weight")
    cut: float
        Manually assigned cutoff point.

    Returns
    -------
    G : NetworkX graph
        A graph after surgery.
    """
    G = G_origin.copy()
    w = nx.get_edge_attributes(G, weight)

    assert cut >= 0, "Cut value should be greater than 0."
    if not cut:
        cut = (max(w.values()) - 1.0) * 0.6 + 1.0  # Guess a cut point as default

    to_cut = []
    for n1, n2 in G.edges():
        if G[n1][n2][weight] > cut:
            to_cut.append((n1, n2))
    print("*************** Surgery time ****************")
    print("* Cut %d edges." % len(to_cut))
    G.remove_edges_from(to_cut)
    print("* Number of nodes now: %d" % G.number_of_nodes())
    print("* Number of edges now: %d" % G.number_of_edges())
    cc = list(nx.connected_components(G))
    print("* Modularity now: %f " % nx.algorithms.community.quality.modularity(G, cc))
    print("* ARI now: %f " % ARI(G, cc))
    print("*********************************************")

    return G

def draw_graph(G, clustering_label="club"):
    """
    A helper function to draw a nx graph with community.
    """
    complex_list = nx.get_node_attributes(G, clustering_label)

    le = preprocessing.LabelEncoder()
    node_color = le.fit_transform(list(complex_list.values()))

    nx.draw_spring(G, nodelist=G.nodes(),
                   node_color=node_color,
                   cmap=plt.cm.rainbow,
                   alpha=0.8, with_labels=True)

@article{ni2019community,
  title={Community detection on networks with ricci flow},
  author={Ni, Chien-Chun and Lin, Yu-Yao and Luo, Feng and Gao, Jie},
  journal={Scientific reports},
  volume={9},
  number={1},
  pages={1--12},
  year={2019},
  publisher={Nature Publishing Group}
}
```

以上