UMAP 0.4の新機能で遊ぶ（プロット、非ユークリッド空間への埋め込み、逆変換）

UMAPがバージョンアップしてv0.4が公開された。

2020/02/10現在では、```pip install --pre umap-learn```でバージョンを上げることができる。

疎行列をそのまま入力できたりいろんな機能が追加されているらしいけど、ここではプロット機能、非ユークリッド空間への埋め込み、逆変換を試してみる。

データだけ変えてほぼドキュメントに書いてあるコード例そのままやってるだけなので、それぞれについて詳しくは[UMAPドキュメント](https://umap-learn.readthedocs.io)へ。

## データ

[PARC](https://github.com/ShobiStassen/PARC)のレポジトリに置いてあったscRNA-seqのデータセットとアノテーション（Zheng et al., 2017, 10X PBMC）を使って実験する。68,579細胞、事前にPCAで50次元に圧縮済み。気軽にやるにはちょっと大きすぎるデータなので適当に1万細胞くらいに落として使う。


```python
import numpy as np
import pandas as pd
import umap
```

データ読み込み。1万細胞、50主成分。


```python
dat = np.loadtxt('./data/pca50_pbmc10k.txt', delimiter=',')
dat.shape
```




    (10000, 50)



アノテーションは以下のような感じ。


```python
labels = []
for line in open('./data/zheng17_annotations_10k.txt'):
    labels.append(line.rstrip())

import collections
import pprint
pprint.PrettyPrinter(indent=4).pprint(collections.Counter(labels))
```

    Counter({   'CD4+/CD45RA+/CD25- Naive T': 1853,
                'CD8+ Cytotoxic T': 1694,
                'CD8+/CD45RA+ Naive Cytotoxic': 1599,
                'CD4+/CD45RO+ Memory': 1425,
                'CD56+ NK': 1211,
                'CD4+/CD25 T Reg': 828,
                'CD14+ Monocyte': 595,
                'CD19+ B': 583,
                'Dendritic': 127,
                'CD4+ T Helper2': 55,
                'CD34+': 30})


とりあえずは普通にUMAPしておく。


```python
model = umap.UMAP(verbose=True)
model.fit(dat)
```

## プロット機能

umap.plotで新たに、umap単体でプロットができるようになった。

matplotlibで描くのとほとんど変わらないけど、手軽ではある。

umap.plotを使う場合はumapと別にdatashader, bokeh, holoviewsが必要なので注意。それぞれインストールしておく。


```python
import umap.plot
```

学習させたモデルのインスタンスを与えると、scatter plotを描いてくれる。


```python
umap.plot.points(model)
```
![output_14_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/bdc66175-a804-3ee2-e093-9f4c0dbf8aeb.png)

また、```labels```のパラメータにそれぞれの点のラベルデータを与えると色分けして描いてくれる。

```python
umap.plot.points(model, labels=labels)
```
![output_16_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/c88e6060-d645-4021-e981-48b2aeccecd1.png)


与えるデータは連続値でもOK。ここでは適当なデータがないので、元データの平均値を与える。

背景やカラーマップなどの組み合わせを用意してくれてるので、```theme```パラメータで好きなものを選んで使える。

選べるthemeは次の9種類。

```python
umap.plot._themes.keys()
```

    dict_keys(['fire', 'viridis', 'inferno', 'blue', 'red', 'green', 'darkblue', 'darkred', 'darkgreen'])


```python
umap.plot.points(model, 
                 values=dat.mean(axis=1),
                 theme='viridis')
```

![output_19_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/3c7353ee-5bfe-aca5-1de4-24c6afc259f6.png)

また、Bokehを使ってインタラクティブなプロットを描くことも可能。

```python
umap.plot.output_notebook()
```

マウスオーバーしたときに表示する情報は事前に```pandas.dataframe```のかたちで用意しておく。

```python
df_labels = pd.DataFrame(labels, columns=['celltype'])
p = umap.plot.interactive(model, labels=labels, 
                          hover_data=df_labels, point_size=2)
umap.plot.show(p)
```

<img width="796" alt="output_19_2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/e26dce34-7780-5521-5639-e5aae7513998.png">


さらに、UMAP埋め込みのときに使われたneighborhood graphを可視化してくれる機能もある。エッジの重みもグラデーションで表示してくれる。どういったconnectivityが学習されたのかを検証するときに使えるかもしれない。


```python
umap.plot.connectivity(model, show_points=True,
                      edge_cmap='viridis')
```

![output_26_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/b3d1ff5b-8dcb-4991-aae5-e2f889631e68.png)


他にも、様々な指標で埋め込みを診断するためのdiagnostic plot機能があるらしい。

## 非ユークリッド空間への埋め込み

デフォルトではUMAPはユークリッド空間に埋め込む（低次元空間のユークリッド距離をターゲットに最適化する）が、球面など他のタイプの空間に埋め込むことも可能らしい。

これは、```output_metric```パラメータで低次元側の距離計算手法を指定することで可能。

まず、球面に埋め込んでみる。球面埋め込みの場合は、Haversine式を指定する。


```python
sphere_mapper = umap.UMAP(output_metric='haversine')
sphere_mapper.fit(dat)
```

結果の座標は[球面座標系](https://ja.wikipedia.org/wiki/球面座標系)で出てくるので、そのまま描いてもよくわからない。なので教科書通りに直交座標系に変換してからプロットしてみる。

```python
x = np.sin(sphere_mapper.embedding_[:, 0]) * np.cos(sphere_mapper.embedding_[:, 1])
y = np.sin(sphere_mapper.embedding_[:, 0]) * np.sin(sphere_mapper.embedding_[:, 1])
z = np.cos(sphere_mapper.embedding_[:, 0])
```

このプロットは勝手にやってくれないみたいなので、自分でmatplotlibで描画。

```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
sns.set(style='white')

categories = sorted(list(set(labels)))
label_ids = [categories.index(l) for l in labels]

fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c=label_ids, cmap='Spectral')
```
![output_36_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/0b485b52-0d99-c8c8-8b61-0b15133a0f9f.png)

セルタイプごとにわりときれいに、球面上にまとまってくれた。

3次元でちょっとわかりづらいので、このまま2次元に展開してみる。

```python
x = np.arctan2(x, y)
y = np.arccos(z)
fig = plt.figure(figsize=(12, 12))
plt.scatter(x, y, c=label_ids, cmap='Spectral')
```

![output_39_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/f5daaa69-be3c-a333-a443-d6630ab196a2.png)


左右は本来つながってる。

どういったときに使えばいいのかよくわからないけど、データが本質的に周期的な性質を持つときなどに有効？

また、双曲空間への埋め込みも紹介されている。これについてはほんとによくわからないのでドキュメントの流れそのままやってみる。[ポアンカレ円板モデル](https://ja.wikipedia.org/wiki/ポワンカレの円板モデル)そのものは最適化が難しいらしいので、[Hyperboloid model](https://en.wikipedia.org/wiki/Hyperboloid_model)をターゲットに学習しているらしい。

```python
hyperbolic_mapper = umap.UMAP(output_metric='hyperboloid')
hyperbolic_mapper.fit(dat)

x = hyperbolic_mapper.embedding_[:, 0]
y = hyperbolic_mapper.embedding_[:, 1]
z = np.sqrt(1 + np.sum(hyperbolic_mapper.embedding_**2, axis=1))

disk_x = x / (1 + z)
disk_y = y / (1 + z)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111)
ax.scatter(disk_x, disk_y, c=label_ids, cmap='Spectral')
boundary = plt.Circle((0,0), 1, fc='none', ec='k')
ax.add_artist(boundary)
ax.axis('off');
```
![output_45_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/13136fdf-83e0-252f-5d40-38dd34894ac7.png)

## 逆変換

埋め込みが学習された低次元側の座標から、対応する高次元サンプルを生成する手法。

VAEみたいな生成モデルと同じように使うのは難しいかもしれないけど、埋め込まれた低次元がどんな空間なのかざっと確認するのには便利かも。

高次元サンプル生成についてはシングルセルのデータはあんまりおもしろくないので、ここでは[Kuzushiji-MNIST](http://codh.rois.ac.jp/kmnist/)のデータ（を適当に1万個ランダムに選んだもの）でやってみる。

```python
dat = np.load('./data/kmnist-train-imgs_10k.npy')
dat.shape
```
    (10000, 784)

```python
labels = np.load('./data/kmnist-train-labels_10k.npy')
labels
```
    array([4, 5, 0, ..., 6, 9, 0], dtype=uint8)

まずは普通にUMAPを計算してみる。

```python
model = umap.UMAP(n_epochs=500, verbose=True).fit(dat)

umap.plot.points(model, labels=labels)
```
![output_52_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/5db19110-7bd9-f635-d4ca-e10846cc792a.png)

この空間から、補間しておもしろそうな領域について逆変換してみたいので、

クラスタ8の左上、クラスタ9の左下、クラスタ0の右上、クラスタ3の右下あたりの点を選んで、テスト用の点を100個作ってみる。

```python
x = model.embedding_

top_left = x[labels == 8, :][x[labels == 8, 0].argmin()]
btm_left = x[labels == 9, :][x[labels == 9, 1].argmin()]
top_right = x[labels == 0, :][x[labels == 0, 0].argmax()]
btm_right = x[labels == 3, :][x[labels == 3, 1].argmin()]

test_pts = np.array([
    (top_left*(1-x) + top_right*x)*(1-y) +
    (btm_left*(1-x) + btm_right*x)*y
    for y in np.linspace(0, 1, 10)
    for x in np.linspace(0, 1, 10)
])

print(top_left)
print(btm_left)
print(top_right)
print(btm_right)
```

    [-3.0056033  9.982167 ]
    [2.0912035 2.2021403]
    [14.147088  10.8581085]
    [10.669375   2.6710818]

逆変換を実行する。```inverse_transform```関数に調べたい点の座標を与える。

```python
inv_transformed_points = model.inverse_transform(test_pts)

from matplotlib.gridspec import GridSpec

fig = plt.figure(figsize=(12,6))
gs = GridSpec(10, 20, fig)
scatter_ax = fig.add_subplot(gs[:, :10])
kuzushiji_axes = np.zeros((10, 10), dtype=object)
for i in range(10):
    for j in range(10):
        kuzushiji_axes[i, j] = fig.add_subplot(gs[i, 10 + j])

scatter_ax.scatter(model.embedding_[:, 0], model.embedding_[:, 1],
                   c=labels.astype(np.int32), cmap='Spectral', s=0.1)
scatter_ax.set(xticks=[], yticks=[])

scatter_ax.scatter(test_pts[:, 0], test_pts[:, 1], marker='x', c='k', s=15)

for i in range(10):
    for j in range(10):
        kuzushiji_axes[i, j].imshow(inv_transformed_points[i*10 + j].reshape(28, 28), cmap='viridis')
        kuzushiji_axes[i, j].set(xticks=[], yticks=[])
```
![output_60_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/176102/a67e74c5-ae23-5e70-7cf9-151f628c0b97.png)


ちょっと粗いが、軸に沿ってどのように遷移するのか、また、クラスタどうしがなぜ近い距離に配置されているのかがなんとなく類推できる。