静かなる名辞

pythonとプログラミングのこと



【python】sklearnのSparsePCAを使ってみる

はじめに

 SparsePCAというものがあることを知ったので、使ってみようと思います。

SparsePCAとは?

 その名の通り、スパースな主成分分析です。スパースな主成分ベクトルを推定します。

Sparse PCA - Wikipedia

 原理などは理解しないで、カジュアルに使えるかどうか試してみるだけという趣旨です。なので「どうやって動いているの?」という質問には答えられません。許してください。

sklearnの実装

 きっちり存在しています(存在しなかったらこんな記事は書きませんが)。

sklearn.decomposition.SparsePCA — scikit-learn 0.20.0 documentation

 主要なパラメータとしては、以下のものがあります。

  • n_components

 PCAのと同じです。

  • alpha

 スパースPCAのキモで、L1正則化の強さを調整できます。

  • ridge_alpha

 こちらはtransformの際に使われるリッジ回帰(L2正則化)の正則化パラメータです。なんでリッジを使うのかは、実のところよくわかりません。

 このパラメータがあるということは、最適化とか勾配法的なもので推定するのだな、というくらいに思っておきます。

  • normalize_components

 主成分ベクトルのノルムを1にするかどうか。Trueにしておくと良いと思います。

 結果に大きな影響を及ぼすのは上くらいだと思います。他のパラメータについてはドキュメントを参照してください。

実験

 今回はwineデータセットでやってみました。素のPCAでやった場合、alphaを0.5と5にした場合の結果をバイプロットで示します。

import matplotlib.pyplot as plt
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA, SparsePCA

def biplot(X_2d, components, target, ax):
    r1 = 5
    r2 = 1.01
    for i, coef in enumerate(components.T):
        ax.arrow(0, 0, coef[0]*r1, coef[1]*r1, color='r')    
        ax.text(coef[0]*r1*r2, coef[1]*r1*r2, i, color='b', fontsize=8)

    ax.scatter(X_2d[:,0], X_2d[:,1], c=target, cmap="rainbow")

def main():
    wine = load_wine()
    ss = StandardScaler()
    X = ss.fit_transform(wine.data)

    pca = PCA(n_components=2)
    spca = SparsePCA(n_components=2,
                     max_iter=3000,
                     n_jobs=-1,
                     normalize_components=True)
    
    fig, axes = plt.subplots(figsize=(12, 6), nrows=1, ncols=3)

    X_2d = pca.fit_transform(X)
    biplot(X_2d, pca.components_, wine.target, axes[0])
    axes[0].set_title("PCA")

    for i,alpha in zip([1, 2], [0.5, 5]):
        spca.set_params(alpha=alpha)
        X_2d = spca.fit_transform(X)
        biplot(X_2d, spca.components_, wine.target, axes[i])
        axes[i].set_title("SPCA alpha={:.2f}".format(alpha))
    plt.savefig("result.png")

    # 図と突き合わせて確認するために特徴量の名前を出力しておく
    for i, name in enumerate(wine.feature_names):
        print(i, name)

if __name__ == "__main__":
    main()

 max_iterをきもち高めにしましたが、結果は数秒程度で出ました。

result.png
result.png

0 alcohol
1 malic_acid
2 ash
3 alcalinity_of_ash
4 magnesium
5 total_phenols
6 flavanoids
7 nonflavanoid_phenols
8 proanthocyanins
9 color_intensity
10 hue
11 od280/od315_of_diluted_wines
12 proline

 とりあえず、PCAの結果とSparsePCAの結果で左右が反転しているのに注意。あとは見ての通りで、alpha=0.5で一部の係数が主成分にべたっと張り付くようになり、alpha=5では大半の係数が主成分に張り付いています。これがSparsePCAの効果で、結果の解釈が容易になるということらしいです(この次元数だとあまり威力はありませんが、高次元では活躍しそうです)。

 ワインにはあまり詳しくないので、今回は結果を細かく解釈することはしませんが……。

まとめ

 使えることがわかりました。

【python】sklearnのRFE(Recursive Feature Elimination)を使ってみる

はじめに

 RFE(Recursive Feature Elimination)というものがあることを知ったので試してみたいと思いました。

 RFEは特徴選択の手法で、その名の通り再帰的にモデルを再構築しながら特徴を選択するという特色があります。

sklearn.feature_selection.RFE — scikit-learn 0.20.0 documentation

RFEとはなにか

 RandomForestのような特徴重要度の出るestimatorがあったとして、たとえば64次元から32次元に特徴選択したいと思ったとします(元の次元数、選択後の次元数は何次元でも良いのですが、あとで実際にこの数字でやるので例として先に示します)。

 真っ先に思いつくのは特徴重要度を見て重要度の低い下位32次元を捨てることですが、ゴミが大量に入った状態では特徴重要度がうまく計算できていないかもしれません。本来有用な特徴を捨ててしまったり、逆に本来ゴミな特徴を残してしまったりする可能性があります。

 一方、RFEは「モデルを作って特徴重要度を計算する」→「下位n次元を捨てる」→「下位n次元を捨てた特徴量でまたモデルを作る」→「下位n次元を捨てる」→……と繰り返していき、望む次元数になるまでこれを続けます。直感的には、うまくいきそうだけど計算コストは大きそう、という気がします。

sklearnのRFEの使い方

 記事の冒頭にも貼りましたが、ドキュメントはここです。

sklearn.feature_selection.RFE — scikit-learn 0.20.0 documentation

 モデルに渡せる引数としては、以下のものがあります。

  • estimator

 好きなestimatorを渡せます……と言いたいところですが、coef_かfeature_importances_を持っている必要があります。coef_のときはデータをスケーリングしないとダメじゃないのかとか、SVMのcoef_はカーネルが絡むから特殊なんじゃないのかとか、色々と微妙な懸念があります(深く追求はしません)。一番安心して使えるのはRandomForestといった決定木のアンサンブル系など、feature_importances_があるようなestimatorです。

  • n_features_to_select

 最終的に選択したい特徴量の次元数です。デフォルトはNoneで、そうすると元の次元数の半分にされます。

  • step

 一度モデルを再構築するたびにどれだけ次元数を減らすか。これを大きくすると速く計算できますが、そのぶん雑な感じの選択になりそうな気がします。デフォルトは1です。

  • verbose

 設定すると学習過程をprintしてくれます。デフォルトは0で何も表示されません。1にすると学習過程が出てきます。1にしたときと2以上にしたときとの違いは私にはわかりませんでした。

 使い方は、適当にモデル作って、fit&predictで動かすだけです。

実際にやってみる

 digitsデータセットの分類をやります。digitsデータセットは64次元なので、

  • そのまま64次元で学習&予測
  • 64次元で学習させた際の特徴重要度を用いて下位32次元を切り捨て、学習&予測
  • RFEを使って32次元まで落とし、学習&予測

 の3パターン試してみました。なお、step=4としました。

 コードを以下に示します。

from sklearn.datasets import load_digits
from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report as clf_report

def main():
    digits = load_digits()
    
    X_train, X_test, y_train, y_test = train_test_split(
        digits.data, digits.target, random_state=0)

    # とりあえず普通にRFCにそのまま入れてみる
    rfc = RFC(n_estimators=100)
    rfc.fit(X_train, y_train)
    preds = rfc.predict(X_test)

    print("通常のRandomForest")
    print(clf_report(y_test, preds, digits=3))

    # 重要度の低い下位32次元を落としてみる
    idx = rfc.feature_importances_.argsort()[32:]
    X_train_32 = X_train[:,idx]
    X_test_32 = X_test[:,idx]
    
    rfc.fit(X_train_32, y_train)
    preds = rfc.predict(X_test_32)

    print("特徴重要度が低いものを32個削除")
    print(clf_report(y_test, preds, digits=3))

    # RFEを使ってみる
    rfec = RFE(rfc, step=4, verbose=1)
    rfec.fit(X_train, y_train)
    preds = rfec.predict(X_test)
    print("RFE + RFC", rfec.n_features_)
    print(clf_report(y_test, preds, digits=3))

if __name__ == "__main__":
    main()

 結果は、

通常のRandomForest
              precision    recall  f1-score   support

           0      0.974     1.000     0.987        37
           1      0.956     1.000     0.977        43
           2      1.000     0.955     0.977        44
           3      0.938     1.000     0.968        45
           4      1.000     0.974     0.987        38
           5      1.000     0.958     0.979        48
           6      1.000     1.000     1.000        52
           7      0.980     1.000     0.990        48
           8      1.000     0.958     0.979        48
           9      0.957     0.957     0.957        47

   micro avg      0.980     0.980     0.980       450
   macro avg      0.980     0.980     0.980       450
weighted avg      0.981     0.980     0.980       450

特徴重要度が低いものを32個削除
              precision    recall  f1-score   support

           0      0.974     1.000     0.987        37
           1      0.956     1.000     0.977        43
           2      1.000     0.909     0.952        44
           3      0.917     0.978     0.946        45
           4      0.974     0.974     0.974        38
           5      0.958     0.958     0.958        48
           6      1.000     0.981     0.990        52
           7      0.959     0.979     0.969        48
           8      0.957     0.938     0.947        48
           9      0.957     0.936     0.946        47

   micro avg      0.964     0.964     0.964       450
   macro avg      0.965     0.965     0.965       450
weighted avg      0.965     0.964     0.964       450

Fitting estimator with 64 features.
Fitting estimator with 60 features.
Fitting estimator with 56 features.
Fitting estimator with 52 features.
Fitting estimator with 48 features.
Fitting estimator with 44 features.
Fitting estimator with 40 features.
Fitting estimator with 36 features.
RFE + RFC 32
              precision    recall  f1-score   support

           0      0.974     1.000     0.987        37
           1      0.977     1.000     0.989        43
           2      1.000     0.909     0.952        44
           3      0.936     0.978     0.957        45
           4      1.000     0.974     0.987        38
           5      0.979     0.958     0.968        48
           6      1.000     1.000     1.000        52
           7      0.960     1.000     0.980        48
           8      0.979     0.958     0.968        48
           9      0.938     0.957     0.947        47

   micro avg      0.973     0.973     0.973       450
   macro avg      0.974     0.973     0.973       450
weighted avg      0.974     0.973     0.973       450

 微妙に単純に下位32次元切り捨てよりRFEの方が良さげにみえますが、はっきり言ってこれはほぼ誤差です。試行によって逆転するときもあります。何回か回した感じ、微妙にRFEの方が優秀そうな気はするので、たぶん100回くらい回して検定すれば有意差は出るんじゃないかなぁ、と思います。面倒なのでやりませんけど。

 何しろ元の分類精度が良いせいで、あまり際立った結果にならない感がありますが、元論文のアブストではベースライン86%に対して98%の精度を達成した、といった数字が示されています。実際のところ、どうなんでしょうね・・・。

 なお、16次元まで削ったところまあまあ顕著な差(それでも1%ない程度ですが、逆転する頻度は4回に1回程度まで減少した)になったことを言い添えておきます。あと、stepsに対しては性能の変化はあまりない気がしました。

 データが差が出づらいものな可能性はありますが、ちょっと微妙すぎるかな・・・と思い、step=1にして10次元まで削ったところ、やっと納得できるような差(見てわかるレベル、数%)が生じ、RFE有利になりました。大きく次元を落とす分には良さそうです。ただ、スコアそのものが0.9前後になってしまうので、絶対性能の観点からするとどうだろうという気はします。

可視化してみる

 digitsはMNISTと同様の、8*8の数字の画像のデータなので、どの特徴が有効になったのかを画像として可視化してみました。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.model_selection import train_test_split

def main():
    digits = load_digits()
    
    X_train, X_test, y_train, y_test = train_test_split(
        digits.data, digits.target, random_state=0)

    rfc = RFC(n_estimators=100)
    rfc.fit(X_train, y_train)

    idx_a = np.array([True]*64)
    idx_a[rfc.feature_importances_.argsort()[:32]] = False
    idx_a = idx_a.reshape(8, 8)

    rfec = RFE(rfc, step=4, verbose=1)
    rfec.fit(X_train, y_train)
    idx_b = rfec.support_.reshape(8, 8)

    fig, axes = plt.subplots(nrows=1, ncols=2)

    axes[0].imshow(idx_a)
    axes[1].imshow(idx_b)
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

可視化した結果
可視化した結果

 左が一気に落とした場合、右がRFEですが、……気のせいレベルの差異しかないですね。ただの乱数のいたずらじゃないの? というレベル。

 step=1として10次元まで落とすと、こうなりました。

step=1で10次元まで落とした場合
step=1で10次元まで落とした場合

 今度は顕著な差がついた・・・ように見えますが、2箇所入れ替わっているだけです。この内容について考察することは控えたいと思いますが(よくわからないので)、とにかくこれくらいの差になります。10次元だと選択した特徴の数が少ないので、ある程度スコアの差にもなってくると思います。

まとめ

 ぶっちゃけ、本当に良いの? という程度ですが、とにかく特徴選択に使えます。

 stepsを大きめにすればたとえば元の10倍とかその程度の処理時間で済むようにすることもできる訳で、気楽に使ってみる分には良いんじゃないでしょうか。

 すごく良い、というものではなさそうです。

ループで辞書の要素を削除しようと思ったらRuntimeError: dictionary changed size during iteration

前提

  • ループで条件に従って辞書の全要素を舐め、条件が真になる要素を削除したい
  • あくまでもin-placeで処理したい(今回はdel文で書いていた)

 要するにこんなコード。

d = {v:"hoge!"*v for v in range(5)}
# => {0: '', 1: 'hoge!', 2: 'hoge!hoge!', 3: 'hoge!hoge!hoge!', 4: 'hoge!hoge!hoge!hoge!'}
for k in d.keys():
    if len(d[k]) > 10: 
        del d[k]

 そしてこうなる。

RuntimeError: dictionary changed size during iteration

原因

 dict.keys()でループ中に辞書を変更してはいけない。これは辞書ビューオブジェクトと言い、他の言語でいうところのジェネレータである(pythonはジェネレータの定義が特殊なので単にジェネレータと呼ぶのは憚られるけど)。

辞書の項目の追加や削除中にビューをイテレートすると、 RuntimeError を送出したり、すべての項目に渡ってイテレートできなかったりします。

https://docs.python.jp/3/library/stdtypes.html#dict-views

 他に、dict.items()やdict.values()等でループしても同じエラーを拝めるはずである。

対策

 他の言語でいうところのジェネレータになっているのが悪いので、ループ前にlistに変換してやる。

d = {v:"hoge!"*v for v in range(5)}
for k in list(d.keys()):
    if len(d[k]) > 10: 
        del d[k]
print(d)  # => {0: '', 1: 'hoge!', 2: 'hoge!hoge!'}

 dict.values()やdict.items()でループする場合も同じ対処でいけます。

 めでたしめでたし。

【python】numpy配列の結合方法まとめ

はじめに

 複数のnumpy配列を一つにまとめたい、連結したい、結合したいというシチュエーションはよくあると思います。numpyでは配列を結合してまとめる様々な方法が存在します。

 色々あるのは嬉しいことですが、「多すぎて覚えきれんわ」と思ったので備忘録としてまとめておきます。

 なお、公式ドキュメントにおける配列の結合関係の記述は、

Array manipulation routines — NumPy v1.15 Manual

 に一覧があり、七種類の関数が存在しています。

 また、この記事の内容を理解するにはnumpyのaxisについて知っている必要があります。ご存知でない方は下のリンクなどを参考にしてください。

NumPyの軸(axis)と次元数(ndim)とは何を意味するのか - DeepAge

NumPyでのaxis指定 - Qiita

 目次

方法一覧

listなどに入れてnumpy配列に変換する

 listなど(tupleとかでも可)に入れてnumpy配列に変換することができます。普通に想像する通りの動作になります。

>>> import numpy as np  # この記事の以下の記述では同様にnumpyをimportしていることを前提とします
>>> a = np.array([1,2])
>>> b = np.array([3,4])
>>> np.array([a,b])
array([[1, 2],
       [3, 4]])
>>> np.array([[a],[b]])
array([[[1, 2]],

       [[3, 4]]])
>>> c = np.array([[1,2],[3,4]])
>>> d = np.array([[5,6],[7,8]])
>>> np.array([c,d])
array([[[1, 2],
        [3, 4]],

       [[5, 6],
        [7, 8]]])

 たとえば画像のlistをnumpy配列として扱いたい……というときなどはこの方法で構わないと思います。

 なお、ndimが揃わないとエラーになります。

>>> np.array([a,c])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: could not broadcast input array from shape (2,2) into shape (2)

 ndimが同じでshapeが揃わない場合、エラーにこそならないもののobject型の配列という扱いになり、うまく変換できません。

>>> np.array([np.array([1,2]), np.array([1,2,3])])
array([array([1, 2]), array([1, 2, 3])], dtype=object)

numpy.concatenate

 numpy.concatenateはすでに存在するaxis上で配列を結合します。画像を横や縦に並べたりするイメージです。

 第一引数に配列を格納したシーケンス(listやtupleなど)を渡します。第二引数axisはオプションで、デフォルトは0です。第三引数outはdestinationを指定できますが、通常の用途では必要ないでしょう。

numpy.stack — NumPy v1.15 Manual

>>> a = np.array([[1,2],[3,4]])
>>> b = np.array([[5,6],[7,8]])
>>> np.concatenate([a,b])
array([[1, 2],
       [3, 4],
       [5, 6],
       [7, 8]])
>>> np.concatenate([a,b], axis=1)
array([[1, 2, 5, 6],
       [3, 4, 7, 8]])

numpy.stack

 numpy.stackは新しいaxisを作成し、そのaxisで配列を連結します。画像を積み重ねて立体にする感じ……でしょうか。

 引数はnumpy.concatenateと同じです。axisの扱いが異なっており、新しい次元をどの方向に作るかを指定できます。うまく言葉で説明できないので、結果を見てもらった方がわかりやすいと思います。

numpy.stack — NumPy v1.15 Manual

>>> np.stack([a,b], axis=0)
array([[[1, 2],
        [3, 4]],

       [[5, 6],
        [7, 8]]])
>>> np.stack([a,b], axis=1)
array([[[1, 2],
        [5, 6]],

       [[3, 4],
        [7, 8]]])
>>> np.stack([a,b], axis=2)
array([[[1, 5],
        [2, 6]],

       [[3, 7],
        [4, 8]]])

 ……なんとなくわかるような、わからないようなって感じかな。なお、axis=0の場合は、listに入れてarrayに変換するのと同じです。

numpy.column_stack

 numpy.column_stackは配列を列方向に積み重ねます。引数は一次元ないし二次元の配列のシーケンスです。concatenate, stackと比べると汎用性の低い関数です。

numpy.column_stack — NumPy v1.15 Manual

 一次元配列の場合はこのような動作になります。

>>> a = np.array([1,2,3])
>>> b = np.array([4,5,6])
>>> np.column_stack([a,b])
array([[1, 4],
       [2, 5],
       [3, 6]])

 二次元配列では、また挙動が異なります。

>>> a = np.array([[1,2,3]])
>>> b = np.array([[4,5,6]])
>>> np.column_stack([a,b])
array([[1, 2, 3, 4, 5, 6]])

 このようにあくまでも二次元配列の列方向に向かってスタックします。また、三次元以上の配列に対して適用することはできません。

numpy.vstack

 numpy.vstackはnumpy配列を垂直方向に結合します。使いやすいので、使った記憶がある人も多いと思います。

numpy.vstack — NumPy v1.15 Manual

 一次元配列は(1, N)にreshapeしてから結合するような挙動になります。

>>> a = np.array([1,2,3])
>>> b = np.array([4,5,6])
>>> np.vstack([a,b])
array([[1, 2, 3],
       [4, 5, 6]])

 二次元以上であればaxis=0での結合です。

>>> a = np.array([[1,2,3]])
>>> b = np.array([[4,5,6]])
>>> np.vstack([a,b])
array([[1, 2, 3],
       [4, 5, 6]])

numpy.hstack

 numpy.hstackはnumpy配列を水平方向に結合します。これもvstackと並んでよく使われる関数だと思います。

numpy.hstack — NumPy v1.15 Manual

 一次元配列の場合はaxis=0で結合されます。二次元以上であれば、axis=1で結合されます。何次元の配列でも使えるようですが、あまりndimが深い配列をこれで操作することもないでしょう。

>>> a = np.array([1,2,3])
>>> b = np.array([4,5,6])
>>> np.hstack([a,b])
array([1, 2, 3, 4, 5, 6])
>>> a = np.array([[1],[2],[3]])
>>> b = np.array([[4],[5],[6]])
>>> np.hstack([a,b])
array([[1, 4],
       [2, 5],
       [3, 6]])

numpy.dstack

 numpy.dstackはnumpy配列を深さ方向に結合します。要はaxis=2の結合であり、axis=0で結合するvstack、axis=1で結合するhstackに続いています。

numpy.dstack — NumPy v1.15 Manual

 直感的にわかりづらいと思うので、実例を示します。

 一次元配列のlistを渡した場合。

>>> a = np.array([1,2,3])
>>> b = np.array([4,5,6])
>>> np.dstack([a,b])
array([[[1, 4],
        [2, 5],
        [3, 6]]])

 なんだかよくわからないのでそのまま二次元配列にしてみます。

>>> a = np.array([[1,2],[3,4]])
>>> b = np.array([[5,6],[7,8]])
>>> np.dstack([a,b])
array([[[1, 5],
        [2, 6]],

       [[3, 7],
        [4, 8]]])

 まだわからないので3次元に。

>>> a = np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
>>> b = np.array([[[9,10],[11,12]],[[13,14],[15,16]]])
>>> np.dstack([a,b])
array([[[ 1,  2,  9, 10],
        [ 3,  4, 11, 12]],

       [[ 5,  6, 13, 14],
        [ 7,  8, 15, 16]]])

 あー、なるほど、なんとなくはわかった。

 dstackはvstack, hstackと比較して、使う機会は少ないと思います。

numpy.block

 これは今までの関数とは少し毛色が違います。これを使うことで、垂直に結合してから水平に結合……といった複雑な操作を一発で簡単に書くことができます(たとえば四枚の画像を縦横に並べるときなど)。挙動としては、numpy配列を入れた多重リストを結合します。

numpy.block — NumPy v1.15 Manual

 最も内側のlistはaxis=-1で結合され、その外側はaxis=-2で……と再帰的に結合される仕様になっています。

>>> a = np.array([[1,2],[3,4]])
>>> b = np.array([[5,6],[7,8]])
>>> c = np.array([[9,10],[11,12]])
>>> d = np.array([[13,14],[15,16]])
>>> np.block([[a,b],[c,d]])
array([[ 1,  2,  5,  6],
       [ 3,  4,  7,  8],
       [ 9, 10, 13, 14],
       [11, 12, 15, 16]])

どれを覚えれば良いのか

 特に覚えておく必要があるのは、

  • listに入れて配列に変換する方法
  • すでに存在するaxisで結合するnumpy.concatenate
  • 新しくaxisを作って結合するnumpy.stack
  • ネストされた配列のlistを結合するnumpy.block

 あたりで、他は簡略記法としてvstack, hstackを覚えておけば大抵の状況には対応できると思います。

 column_stackとdstackは使っているコードをほとんど見たことがないので、忘れても構いませんが、そういうものがあるということは頭の片隅に入れておくと良いと思います。

まとめ

 numpy配列の結合方法についてまとめました。これでnumpy配列を自由自在に結合できるようになると思います。

【python】ctypesのcreate_string_buffer()を使ってみる

はじめに

 以前の記事で、ctypesでバイト列や文字列を受け渡しする方法について述べました。

【python】ctypesでバイト列や文字列を受け渡しする - 静かなる名辞

 しかし、ctypesに存在しているcreate_string_buffer()
create_unicode_buffer()には触れませんでした。

 使ってみたら便利だったので紹介します。

これらは何なのか

 変更可能なバッファをpython側で作成します。

 これの良さそうなところは、python側のインタプリタでリソースが管理されるであろうことです。そのため、処理が簡単になります。

 先に断っておきますが、具体的な仕様はドキュメントを参照してください。この記事はあくまでも簡単に使ってみるだけです。

create_string_buffer

 まずC言語で次のような関数を書きます。

lib1.c

void reverse_bytes(int n, char *s1, char *s2) {
  int i;
  
  for (i=0; i<n; i++) {
    s2[i] = s1[n-i-1];
  }
}
// ファイル名がlib1.cなら「gcc -shared -fPIC lib1.c -o lib1.so」のようにコンパイルしてください

 アホみたいに単純なコードですが、これですべてです。s1の内容を反転した内容をs2に書き込みます。

 python側では、以下のような呼び出しを用意します。

import ctypes

lib = ctypes.cdll.LoadLibrary('./lib1.so')
lib.reverse_bytes.argtypes = (ctypes.c_int, ctypes.c_char_p, ctypes.c_char_p)

def reverse_bytes(buf):
    n = len(buf)
    str_buf = ctypes.create_string_buffer(buf)
    lib.reverse_bytes(n, buf, str_buf)
    return str_buf.value

if __name__ == "__main__":
    print(reverse_bytes(b"hoge"))  # => b'egoh'

 これでちゃんと実行されます。

create_unicode_buffer

 これも基本的に同じです。

lib2.c

#include <wchar.h>

void reverse_str(int n, wchar_t *s1, wchar_t *s2) {
  int i;
  
  for (i=0; i<n; i++) {
    s2[i] = s1[n-i-1];
  }
}
// コンパイル手順:$ gcc -shared -fPIC lib2.c -o lib2.so

 python側のコード。

import ctypes

lib = ctypes.cdll.LoadLibrary('./lib2.so')
lib.reverse_str.argtypes = (ctypes.c_int, ctypes.c_wchar_p, ctypes.c_wchar_p)

def reverse_str(buf):
    n = len(buf)
    str_buf = ctypes.create_unicode_buffer(buf)
    lib.reverse_str(n, buf, str_buf)
    return str_buf.value

if __name__ == "__main__":
    print(reverse_str("hoge"))  # => egoh

 わーい直感的に使える!

速度比較

 前回は組み込みに負けていましたが、今度はどうなるでしょうか。

import ctypes
import timeit

lib = ctypes.cdll.LoadLibrary('./lib2.so')
lib.reverse_str.argtypes = (ctypes.c_int, ctypes.c_wchar_p, ctypes.c_wchar_p)

def reverse_str(buf):
    n = len(buf)
    str_buf = ctypes.create_unicode_buffer(buf)
    lib.reverse_str(n, buf, str_buf)
    return str_buf.value

if __name__ == "__main__":
    s = "hogefugaほげふが"*(10**3)
    print(s[::-1] == reverse_str(s))
    print(timeit.timeit(lambda : reverse_str(s), number=10**3))
    print(timeit.timeit(lambda : s[::-1], number=10**3))
 
""" => 結果
True
0.060752346995286644
0.00984983000671491

# 参考:mallocで実装した前回の結果
True
0.046704520999810484
0.009301142999902368
"""

 かえって遅くなる。悲しい。

まとめ

 何はともあれ簡単にpython側のGCに管理を委ねることができるとわかりました。

 NULL終端の文字列等を渡すのであれば、これらを使えば良さそうです。ただし、本当に生のバイナリデータを渡すのであれば、前回の記事の方法でやった方が良いのかな? という気も少しするので、けっこう微妙なところだと思います。

【python】sklearnのLDA(LatentDirichletAllocation)を試してみる

 注意:線形判別分析(LinearDiscriminantAnalysis)ではありません。トピックモデルのLDAです。

はじめに

 LDAといえば、トピックモデルの代表的な手法であり、一昔前の自然言語処理では頻繁に使われていました(最近は分散表現や深層学習に押されて廃れ気味な気もしますが)。

 普通、pythonでLDAといえばgensimの実装を使うことが多いと思います。が、gensimは独自のフレームワークを持っており、少しとっつきづらい感じがするのも事実です。

gensim: models.ldamodel – Latent Dirichlet Allocation

 このLDA、実はsklearnにもモデルがあるので、そっちを試しに使ってみようと思います。

 ライブラリのリンク
sklearn.decomposition.LatentDirichletAllocation — scikit-learn 0.20.0 documentation

何はともあれ分類タスク

 当ブログでは定番になっている20newsgroupsデータセットで文書分類をやってみましょう。

from sklearn.datasets import fetch_20newsgroups
from sklearn.model_selection import train_test_split
from sklearn.feature_extraction.text import CountVectorizer as CV
from sklearn.decomposition import LatentDirichletAllocation as LDA
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

def main():
    news20 = fetch_20newsgroups()

    X_train, X_test, y_train, y_test = train_test_split(
        news20.data[:5000], news20.target[:5000],
        stratify=news20.target[:5000])

    cv = CV(min_df=0.04, stop_words="english")
    lda = LDA(n_components=100, max_iter=30, n_jobs=-1)
    rfc = RFC(n_estimators=500, n_jobs=-1)
    estimators = [("cv", cv), ("lda", lda), ("rfc", rfc)]
    pl = Pipeline(estimators)

    pl.fit(X_train, y_train)
    y_pred = pl.predict(X_test)
    print(classification_report(
        y_test, y_pred, target_names=news20.target_names))

if __name__ == "__main__":
    main()

 見ての通りのモデルです。CountVectorizerでBoWに変換し、LDAで次元削減したあとランダムフォレストにかけて分類します。パラメータは適当にいじって計算時間と分類性能の兼ね合いで決めました。

 で、肝心の分類性能ですが、classification_reportで出した結果を以下に示します。

                          precision    recall  f1-score   support

             alt.atheism       0.30      0.25      0.27        55
           comp.graphics       0.29      0.27      0.28        63
 comp.os.ms-windows.misc       0.58      0.74      0.65        66
comp.sys.ibm.pc.hardware       0.35      0.31      0.33        64
   comp.sys.mac.hardware       0.16      0.09      0.12        64
          comp.windows.x       0.31      0.25      0.28        67
            misc.forsale       0.46      0.62      0.53        66
               rec.autos       0.48      0.57      0.52        65
         rec.motorcycles       0.22      0.25      0.23        69
      rec.sport.baseball       0.25      0.25      0.25        65
        rec.sport.hockey       0.44      0.60      0.51        65
               sci.crypt       0.47      0.66      0.55        67
         sci.electronics       0.37      0.29      0.33        69
                 sci.med       0.25      0.28      0.27        68
               sci.space       0.65      0.63      0.64        62
  soc.religion.christian       0.49      0.56      0.52        66
      talk.politics.guns       0.31      0.28      0.29        58
   talk.politics.mideast       0.53      0.38      0.45        60
      talk.politics.misc       0.43      0.35      0.38        52
      talk.religion.misc       0.36      0.21      0.26        39

               micro avg       0.40      0.40      0.40      1250
               macro avg       0.39      0.39      0.38      1250
            weighted avg       0.39      0.40      0.39      1250

 ちょっとしょぼい・・・かな。何も考えず、BoWのままランダムフォレストに入れたときでも0.7くらいの評価指標になっていたので、しょぼいと言わざるを得ない感じです。

【python】sklearnのfetch_20newsgroupsで文書分類を試す(2) - 静かなる名辞

 トピック数やmax_iterなどのパラメータを上げると改善する傾向ですが、計算時間は増加します。気合い入れてパラメータを適正に選んでぶん回せばたぶんちゃんとした結果になるのだと思いますが、そうするコストは高そうです。

トピックの中身を見てみる

 LDAは分類の前処理に使うより、むしろ出来上がったトピックを見て「はーんこんな風に分かれるのね」と楽しむ方が面白いような気もするので、軽く見てみます。ただし、膨大な文書をLDAで回すのも、その結果を人手で解釈するのも大変なので、

  • 1000文書だけ入れる
  • 50トピックでLDAモデルを構築し、うち先頭10トピックだけ見ることにする
  • 各トピックの大きい係数5個だけ見る

 という手抜き解釈です。

import numpy as np
from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer as CV
from sklearn.decomposition import LatentDirichletAllocation as LDA

def main():
    news20 = fetch_20newsgroups()
    X = news20.data[:1000]
    y = news20.target[:1000]

    cv = CV(min_df=0.04, stop_words="english")
    lda = LDA(n_components=50, max_iter=50, n_jobs=-1)
    lda.fit(cv.fit_transform(X, y))

    feature_names = np.array(cv.get_feature_names())
    for i, component in enumerate(lda.components_[:10]):
        print("component:", i)
        idx = component.argsort()[::-1][:5]
        for j in idx:
            print(feature_names[j], component[j])

if __name__ == "__main__":
    main()

 結果は以下の通り。コメントで解釈を付けてみました。

component: 0  # carとspeedなので車関係
car 105.0307361645338
speed 66.0512803728286
years 34.28189251426528
year 28.634803273219863
ago 22.20999274754642
component: 1  # よくわからん
time 88.93501118238727
david 40.73976076823079
group 37.700769220589336
information 25.679148778591387
services 24.1215050863642
component: 2  # IBMなので恐らくコンピュータ関連
com 178.86471720750257
ibm 103.60651295996851
writes 53.83219596330523
article 38.98071797104704
lines 31.65901015205818
component: 3  # 教育機関とか?
edu 293.2163826198764
cc 106.47895192641978
writes 96.85399760416148
organization 67.46869749413264
article 67.44141504814338
component: 4  # 上のトピックと似ていてeduがトップに来る。なんでしょうねぇ(メールアドレスとかURLだったりして)
edu 303.2418519434582
article 97.83240562135947
writes 89.61827003849298
technology 78.82229006007034
organization 65.80274084351083
component: 5  # windowsとかソフト関連
windows 150.91253925843284
00 91.50205933157577
15 24.611812739517074
software 19.04945545853369
20 14.902153206775973
component: 6  # よくわからない
thought 28.89413050306475
kind 21.98192709084761
david 13.591746503570576
code 8.357254795587956
little 7.6968699914925
component: 7  # メーリングリスト関連?
list 62.55156212975074
address 50.37113047912307
data 45.349827179735904
information 39.6113013625701
mail 33.98010842887469
component: 8  # カネ、政府、ソフトウェア……
money 36.13150161231052
gov 28.618673196459778
software 21.193931015050016
source 20.603010804981693
cost 20.35807225674739
component: 9  # 大学とかのネットワーク関連っぽい
edu 439.0633403498364
university 165.04391386786645
posting 162.13921320408093
nntp 159.2048859320109
host 158.37699187575328

 恐ろしく雑な解釈ですが、一応トピックごとにまとまっているっぽいことはなんとなくわかりました。はっきり内容がわかるのは半分以下ですが……。

sklearnのLDAの欠点

 たぶんperplexityとかCoherenceとかは見れないと思う。

 体感ではgensimのより遅い気がします。ちゃんと測っていないので、あくまでもなんとなくですが。

まとめ

 sklearnでも一応できることはわかりました。

 機能がしょぼいので本格的に使いたければgensimの方がおすすめです。gensimが使えない環境でLDAしないといけなくなったときには代用品として使えなくはありません。

C言語でshellの多段パイプを実装

はじめに

 学校の課題でCでshellもどきを書きました。

 今後、同じ目にあう人のために、「shellの多段パイプをどうやって実装したら良いのか」を記事としてまとめておきます。

 目次

パイプの概要

 shellのパイプとは……という話はさすがに要らないと思いますが、以下のような機能があります。なお、カレントディレクトリに今回の記事で紹介するCのコードtest.cを置いてあるとします(内容については後述します)。

$ cat test.c | head | grep char
char *cmd1[] = {"cat", "test.c", NULL};
char *cmd2[] = {"head", NULL};
char *cmd3[] = {"grep", "char", NULL};
char **cmds[] = {cmd1, cmd2, cmd3};

 「cat test.c」はtest.cの内容を標準出力に吐き出します。が、今回はその出力はパイプによって「head」に繋がれます。「head」は入力の先頭10行を出力するコマンドです。その出力も「grep char」の入力に繋がれて、先頭10行の中でcharにマッチする行だけが出てきます。

 C言語のコードでこれと同じものが動くことが今回の目標です。

使用する関数

 shellのパイプ機能を実装するために最低限必要な関数を示します。なお、上の例のtest.cで察した方もおられるかと思いますが、コマンドの入力やパースは今回省略し、あくまで「実行するとパイプでコマンドをつないで一回動作するプログラム」を作ります。

 更に、エラー処理等も省略し、コードを極力シンプルな形になるまで削ぎ落としてあります。ヘッダファイルは「unistd.h」だけincludeすればコンパイルできます。使う関数はたった6種類です。

 以下に使用する関数の簡単な説明を記述します。あくまでも簡単な説明なので、ちゃんとした説明が必要ならmanなどを読んでください。

int pipe(int pipefd[2])

 名前無しパイプを生成します。pipefdは配列のアドレスを渡し、ファイルディスクリプタを受け取ります。pipefd[1]にデータを書き込むとpipefd[0]から読み出せます。

int close(int fd)

 引数に渡されたファイルディスクリプタを閉じます。

 パイプをうまく機能させようと思うと、必要のないファイルディスクリプタは片っ端から閉じておく必要があります。無駄に開いている読み出し口や書き込み口があると、入力が終わってもEOFが返されません。するとパイプで繋がれたプログラムが終了しないので、いつまでも待ち続ける羽目になります。必要なものだけ開いた状態にするのが鉄則です。

int dup2(int oldfd, int newfd)

 newfdをoldfdのコピーとして作成します。

 と急に言われても何をするのかよくわかりませんが、上で説明したpipe()でパイプを作っておいたとして、

dup2(pipefd[1], 1);

 で標準出力をパイプの書き込み口に繋ぐことができ、同様に

dup2(pipefd[0], 0);

 とすれば標準入力をパイプの読み出し口に繋ぐことができる、ということだけ覚えておけば、今回は十分です。

pid_t fork(void)

 プロセスをforkします。返り値はpid_tという型ですが、これはただの整数型です。forkが成功した場合、pid_tが0なら子プロセス、0以外(実際には子プロセスのPID)なら親プロセスです。失敗すると-1が返ります。

pid_t wait(int *status)

 子プロセスが返るのを待ちます。成功した場合、返り値は子プロセスのPIDで、statusに終了情報が格納されます。

int execvp(const char *file, char *const argv[])

 コマンドを実行します。第一引数はコマンドの文字列、第二引数は引数の配列でNULL終端とする必要があります。

 要するに、

char *cmd1[] = {"ls", NULL};

 と定義しておけば、

execvp(cmd1[0], cmd1);

 でlsが実行できます。

パイプの実装方針

 さて、シェルのパイプを作ることを考えます。ここで考え込んでもあまり良いアイデアは浮かんでこないので、とりあえず実際のコマンドを見ます。

$ cat test.c | head | grep char

 では3つのコマンドを2つのパイプで繋いでいます。要するにコマンド数-1のパイプを作れば良い訳です。

 ということは、パイプを配列で管理するのかな? と一瞬思いますが、それでも確かにできるのですが、ちょっと煩雑そうです。

 もう少し簡単にする方法はないでしょうか? あります。再帰を使います。

 まずメインのプロセスからforkして、パイプを作り、更にforkします。親はstdinをパイプに繋いで右端(右から0番目)のコマンドの実行、子はstdoutをパイプに繋いで更にパイプを作ってforkして、今度は親になった先程の子が右端から1番目のコマンドを実行、子はまたforkしてパイプを作り……と繰り返していって、左端のコマンドに達したら単に実行して終わりです。この手続きは再帰的に行えます。

 「なぜ素直に左端からforkしないの?」と疑問を持つ方もいると思いますが、実は左から始めてもパイプそのものはできます。ただし、execしてしまうとプロセスの制御は呼び出し元に戻ってきません。execの中でexitされると思ってください。

 なので、左端からforkすると左端のコマンドが終了した段階でメインプロセス側のwaitが返り、他のコマンドがまだ実行途中であっても制御が戻ってしまいます。実際にやるとわかりますが、出力の途中でプロンプトが出てきたりして、ちょっと不格好な結果になります。右端からforkすれば、右端のコマンドは最後に終了するので、確実にメインプロセス側でコマンドの終了を検知できます。左端からforkする方法でこの問題を回避しようとすると、何らかの制御手段を追加する必要があります。

 左端からやろうと右端からやろうと、execしてしまう以上、途中では親が子をwaitできないことに違いはありません。ゾンビにならないの? と思うかもしれませんが、この場合は最終的に親が死ぬので、initが引き取ってくれてゾンビになりません。大本の親だけ回収しておけば、それほど気にする必要はありません。

 説明だけ読んでいてもどうなっているのかよくわからないと思うので、図にしてみました。

パイプ実行時のforkの流れ
パイプ実行時のforkの流れ

 この図は上から下に実行されていると思ってください。分岐はfork、合流はwaitでプロセスを看取っていることを示します。分岐の左側がforkの親で、右側が子です。

 cmd3はメインプロセスが看取ります。省略していますが、cmd1とcmd2はcmd3が看取られた後にinitが看取ります。

コード

 実際に書いたコードを以下に示します。上の説明はこのコードを書いてから起こしたものなので、ここまでの内容を読んだ方であれば簡単に理解できると思います。50行ちょっとなので読みやすいはずです。

test.c

#include <unistd.h>

char *cmd1[] = {"cat", "test.c", NULL};
char *cmd2[] = {"head", NULL};
char *cmd3[] = {"grep", "char", NULL};
char **cmds[] = {cmd1, cmd2, cmd3};
int cmd_n = 3;

void dopipes(i) {
  pid_t ret;
  int pp[2] = {};
  if (i == cmd_n - 1) {
    // 左端なら単にexecvp
    execvp(cmds[0][0], cmds[0]);
  }
  else {
    // 左端以外ならpipeしてforkして親が実行、子が再帰
    pipe(pp);
    ret = fork();

    if (ret == 0) {
      // 子プロセスならパイプをstdoutにdup2してdopipes(i+1)で再帰し、
      // 次のforkで親になった側が右からi+1番目のコマンドを実行
      close(pp[0]);
      dup2(pp[1], 1);
      close(pp[1]);
      
      dopipes(i+1);
    }
    else {
      // 親プロセスならパイプをstdinにdup2して、
      // 右からi番目のコマンドを実行
      close(pp[1]);
      dup2(pp[0], 0);
      close(pp[0]);
      
      execvp(cmds[cmd_n-i-1][0], cmds[cmd_n-i-1]);
    }
  }  
}

int main(void) {
  pid_t ret;
  
  ret = fork();
  if (ret == 0)
    dopipes(0);
  else
    wait(NULL);

  return 0;
}

 コンパイルして実行すると(実行ファイルはソースコードと同一ディレクトリに置いてください)、

char *cmd1[] = {"cat", "test.c", NULL};
char *cmd2[] = {"head", NULL};
char *cmd3[] = {"grep", "char", NULL};
char **cmds[] = {cmd1, cmd2, cmd3};

 と最初のshellから打ち込んだパイプコマンドと同じ結果が出力されます。

改良した方が良い点

 とりあえず、システムコールは失敗することもあり得るので、ちゃんとエラー処理しましょう。この記事ではわかりやすさを重視してすべて端折っていますが、

 あとはdup2に関してですが、

// stdoutにdup2
close(pp[0]);
dup2(pp[1], 1);
close(pp[1]);

// ------

// stdinにdup2
close(pp[1]);
dup2(pp[0], 0);
close(pp[0]);

 dup2の際にnewfdが開いていれば勝手に閉じられます。これに失敗する可能性があり、その場合エラー情報は握りつぶされます。なので、stdinとstdoutも明示的に閉じてエラー処理をした方が良いとされます。

まとめ

 これでパイプを実装しないといけなくなっても大丈夫!

参考にしたサイト

シェルの多段パイプを自作してみる | 慶應義塾大学ロボット技術研究会
 こちらは配列で管理する方法でパイプを実装しています。

linux上で動くシェルを自作しています。多段階のパイプを実装方法を教... - Yahoo!知恵袋
 結構ヒントになりました。ここの回答をコードに起こしたようなものでした。