より高速なKModes, KPrototypesを実装

# より高速なKModes, KPrototypesを実装
　こちらの[コード](https://github.com/nicodv/kmodes )を参考に、より高速化したKModes、KPrototypesを実装しました。距離尺度については、カテゴリ特徴量はハミング距離、数値特徴量はユークリッド距離のみを実装しています。

https://github.com/NaoMatch/FasterKModes

## 高速化
 　処理時間がかかる部分をＣで置き換えました。KMeansと同じく処理の重い箇所は、
  
  * 入力データとクラスター中心との間の距離計算（距離行列）
  * クラスター中心の更新

  の二つです。カテゴリ特徴量の場合、前者はSIMDによる置き換えも容易なので、非常に簡易的なものですがそちらも実装しました。数値特徴量では、scikit-learnのKMeansでも用いられているように[GEMMを活用](https://github.com/scikit-learn/scikit-learn/blob/d666202a9349893c1bd106cc9ee0ff0a807c7cf3/sklearn/cluster/_k_means_lloyd.pyx#L200)することで、かなりの高速化が実現できます。おそらくカテゴリ特徴量もGEMM実装のいろいろなテクニックを適用できればもっと速くなるかと思うのですが、難易度が高過ぎたため、着手していません。またSIMD実装は今回実験したデータセットでは意味がありませんでした。
　後者のクラスタ更新は、メモリへの連続アクセス方法が思い浮かばなかったため、SIMD対応はしていません。ただ、カテゴリ数が多い場合には高速になるように（多分）キャッシュを意識した実装になっているようにしてみました。
    
　上記のC言語で実装した関数の一部は、入力データの型と特徴量数に応じて自動的に内部で生成しています。特徴量数が16/32未満ではSIMDを用いる意味がないため、単純な実装を呼び出します。それぞれ以下のファイル内の処理を確認してみてください。
  * 入力データとクラスター中心との間の距離計算（距離行列）
      * `dist_mat_256_hamming.c`(256は列数。入力データの形状により異なる)
          * 入力データ`X`(NxK)とクラスタ中心`C`(MxK)を元に、距離行列`D`(NxM)を計算する
          * 同じく`dist_vec_256_hamming.c`もあるがこちらは、入力データ`X`(NxK)と１個のクラスタ中心`c`(K)を元に、距離ベクトル`d`(N)を計算する
  * クラスター中心の更新
      * `common_funcs.c`

# 使い方
　1. 入力データは`np.ndarray`のみ
　2. KModesの場合は`uint8`, `uint16`のみ。Kprototypesの場合は`float`, `double`も受け付けるが、内部でカテゴリ特徴量部分は`uint8` or `uint16` に変換される。


## KModes
　KModesが受け付ける入力データの型は`uint8`、 `uint16` のみです。`uint16`で６万以上のカテゴリを表現でき、それ以上のカテゴリ数があったとしても、解釈性から言って意味のないものと考えたためです（100以上の時点で十分だと思いますが）。入力データは事前に`sklearn.preprocessing.LabelEncoder`などによる整数値への置換が必要です。
 ```python
import numpy as np
from FasterKModes import FasterKModes

N = 1100 # 行数
K = 31 # 列数
C = 8 # クラスタ数

X = np.random.randint(0, 256, (N, K)).astype(np.uint8)
X_train = X[:1000,:]
X_test = X[1000:,:]

fKModes = FasterKModes(n_clusters=C, init="random", n_init=10)
fKModes.fit(X_train)

labels = fKModes.predict(X_test)
```

 ## KPrototypes
 　KPrototypesは`uint8`、`uin16`に加え、`float`、`double`を受け付けます。KPrototypesは入力`X`とカテゴリ特徴量がカラムのどの位置にあるかを表す`list`を引数とします。
  　`float`か`double`を入力データの型とした際に、カテゴリ特徴量列は最大値に応じて自動的に`uint8`か`uin16`に変換される(小数点以下が無視される)ため、注意が必要です。また、最大値が'uint16'上限を超えている場合にはエラーが出ます。
  
```python
import numpy as np
import FasterKPrototypes as fkp

N = 1100 # 行数
K = 31 # 列数
C = 8 # クラスタ数

X = np.random.randint(0, 256, (N, K)).astype(np.uint8)
X_train = X[:1000,:]
X_test = X[1000:,:]
categorical_features = [1, 3, 5, 7, 9]

fKProto = fkp.FasterKPrototypes(n_clusters=C, init="random", n_init=10)
fKProto.fit(X_train, categorical_feature_indices=categorical_features)

labels = fKProto.predict(X_test)
```




 # 精度・速度比較
 　オリジナルの実装と精度と速度を比較しました。データは、`sklearn.datasets.make_classification`で生成したデータ`X`を、LightGBMで学習し、そのLeafIndexを出力させた出力データ`X'`を用いました。比較の際に変化させた点は、
  * 行数：生成データ`X`の行数
  * 列数：LightGBMの`n_estimators`。生成データ`X`の列数ではない
  * 最大値：`X'`の取り得る最大値。LightGBMの`max_depth`から`pow(2, max_depth)`で計算される＝LeafNodeの取り得る最大の数に相当する
  * クラスタ数：KModesのクラスタ数

 ## 精度
 　精度（cost）は、各データ点とそれらに最も近いクラスター中心の距離までの総和です。完全に同一の値にすることは困難なため、同じ入力データ・同じクラスタ数に対して、同程度の精度が出ればよいものとしています。

  　下記の表がその結果です。初期値はランダムを用いています。おおむね同程度の結果が出ているかと思います。


### KModesのコスト比較
<table border="1" style="border-collapse: collapse; width: 100%;">
  <tr>
    <th rowspan="2" style="text-align: center;">クラスタ数</th>
    <th colspan="2" style="text-align: center;">(1000, 8)</th>
    <th colspan="2" style="text-align: center;">(10000, 8)</th>
    <th colspan="2" style="text-align: center;">(100000, 8)</th>
  </tr>
  <tr>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
  </tr>
  <tr>
    <td style="text-align: right;">2</td>
    <td style="text-align: right;">5,466</td>
    <td style="text-align: right;">5,466</td>
    <td style="text-align: right;">56,118</td>
    <td style="text-align: right;">55,937</td>
    <td style="text-align: right;">554,695</td>
    <td style="text-align: right;">554,695</td>
  </tr>
  <tr>
    <td style="text-align: right;">4</td>
    <td style="text-align: right;">4,353</td>
    <td style="text-align: right;">4,365</td>
    <td style="text-align: right;">43,956</td>
    <td style="text-align: right;">44,325</td>
    <td style="text-align: right;">448,733</td>
    <td style="text-align: right;">445,752</td>
  </tr>
  <tr>
    <td style="text-align: right;">8</td>
    <td style="text-align: right;">3,105</td>
    <td style="text-align: right;">3,086</td>
    <td style="text-align: right;">32,222</td>
    <td style="text-align: right;">33,052</td>
    <td style="text-align: right;">323,439</td>
    <td style="text-align: right;">323,521</td>
  </tr>
  <tr>
    <td style="text-align: right;">16</td>
    <td style="text-align: right;">2,264</td>
    <td style="text-align: right;">2,332</td>
    <td style="text-align: right;">23,035</td>
    <td style="text-align: right;">24,451</td>
    <td style="text-align: right;">241,809</td>
    <td style="text-align: right;">243,484</td>
  </tr>
  <tr>
    <td style="text-align: right;">32</td>
    <td style="text-align: right;">1,632</td>
    <td style="text-align: right;">1,724</td>
    <td style="text-align: right;">17,443</td>
    <td style="text-align: right;">18,388</td>
    <td style="text-align: right;">168,879</td>
    <td style="text-align: right;">182,752</td>
  </tr>
  <tr>
    <td style="text-align: right;">64</td>
    <td style="text-align: right;">1,170</td>
    <td style="text-align: right;">1,147</td>
    <td style="text-align: right;">12,065</td>
    <td style="text-align: right;">12,261</td>
    <td style="text-align: right;">115,768</td>
    <td style="text-align: right;">119,523</td>
  </tr>
  <tr>
    <td style="text-align: right;">128</td>
    <td style="text-align: right;">712</td>
    <td style="text-align: right;">680</td>
    <td style="text-align: right;">7,835</td>
    <td style="text-align: right;">7,723</td>
    <td style="text-align: right;">74,772</td>
    <td style="text-align: right;">78,843</td>
  </tr>
  <tr>
    <td style="text-align: right;">256</td>
    <td style="text-align: right;">402</td>
    <td style="text-align: right;">328</td>
    <td style="text-align: right;">4,992</td>
    <td style="text-align: right;">4,523</td>
    <td style="text-align: right;">48,861</td>
    <td style="text-align: right;">48,304</td>
  </tr>
</table>

### KPrototypesのコスト比較
　クラスタ数16以上はKprototypesでは実行できなかったため、スキップしています。
 https://github.com/nicodv/kmodes#:~:text=Q%3A%20I%27m%20getting%20the%20following%20error%3A%20%22ValueError%3A%20Clustering%20algorithm%20could%20not%20initialize.%20Consider%20assigning%20the%20initial%20clusters%20manually.%22
<table border="1" style="border-collapse: collapse; width: 100%;">
  <tr>
    <th rowspan="2" style="text-align: center;">クラスタ数</th>
    <th colspan="2" style="text-align: center;">(1000, 8)</th>
    <th colspan="2" style="text-align: center;">(10000, 8)</th>
    <th colspan="2" style="text-align: center;">(100000, 8)</th>
  </tr>
  <tr>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
    <th style="text-align: center;">元実装</th>
    <th style="text-align: center;">再実装</th>
  </tr>
  <tr>
    <td style="text-align: right;">2</td>
    <td style="text-align: right;">46,108</td>
    <td style="text-align: right;">46,108</td>
    <td style="text-align: right;">460,290</td>
    <td style="text-align: right;">460,290</td>
    <td style="text-align: right;">6,021,514</td>
    <td style="text-align: right;">6,021,515</td>
  </tr>
  <tr>
    <td style="text-align: right;">4</td>
    <td style="text-align: right;">22,208</td>
    <td style="text-align: right;">21,901</td>
    <td style="text-align: right;">235,031</td>
    <td style="text-align: right;">233,540</td>
    <td style="text-align: right;">1,703,346</td>
    <td style="text-align: right;">1,697,376</td>
  </tr>
  <tr>
    <td style="text-align: right;">8</td>
    <td style="text-align: right;">14,645</td>
    <td style="text-align: right;">14,654</td>
    <td style="text-align: right;">88,984</td>
    <td style="text-align: right;">84,085</td>
    <td style="text-align: right;">1,277,582</td>
    <td style="text-align: right;">1,244,299</td>
  </tr>
</table>



## 速度
　行数、列数、最大値(カテゴリのユニーク数)、クラスタ数をいろいろ変化させてオリジナル実装との速度を比較します。（もっと行数を増やした場合も比較したかったのですが、オリジナル実装の4万行以上の実行は時間がかかりすぎたので、スキップしています。）

### KModes
![可視化結果_max000001.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/323271/cef13f7f-a28f-4c7f-8441-ac5d9d140e28.png)
![可視化結果_max000003.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/323271/6a1284de-a237-441d-bab3-3326a5cb3cdd.png)
![可視化結果_max000015.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/323271/303b4074-bc62-4159-9578-e02999372080.png)
![可視化結果_max000255.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/323271/8015335b-041e-4e01-8eca-b5b586c1c312.png)

 # その他

 ## [KPrototype plus (kpplus)](https://github.com/youbao88/KPrototypes_plus)
 　実験を終わらせてから気づいたのリポジトリです。オリジナル実装にNumbaによる最適化を施しているようです。

## 距離尺度の拡張
　現状ユークリッド距離、ハミング距離しか実装ができていないため、いくつか追加していきたいと考えています。


 以上