ループで辞書の要素を削除しようと思ったら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、合流は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!知恵袋
　結構ヒントになりました。ここの回答をコードに起こしたようなものでした。

【python】numpyで多次元配列のargsortと値の取り出し

はじめに

　numpy配列のargsort()メソッドは値をソートした結果のインデックスの配列を返します。

>>> import numpy as np
>>> a = np.array([2,0,1,8,1,1,0,7])  # 適当な配列を定義
>>> idx = a.argsort()  # argsort
>>> idx  # こんな配列になる
array([1, 6, 2, 4, 5, 0, 7, 3])
>>> a[idx]  # 配列をインデックスとして使う
array([0, 0, 1, 1, 1, 2, 7, 8])

numpy.ndarray.argsort — NumPy v1.15 Manual

　同様のことを多次元配列でも行うことができます。その方法について説明します。

多次元配列でargsort()

　argsort()メソッド自体は多次元配列でも問題なく呼べます。が、ソートなので、axisを指定する必要があります。

>>> a = np.array([[2,0,1,8],[1,1,0,7]])
>>> a
array([[2, 0, 1, 8],
       [1, 1, 0, 7]])
>>> a.argsort()
array([[1, 2, 0, 3],
       [2, 0, 1, 3]])
>>> a.argsort(axis=0)  # axis=0では縦にソート
array([[1, 0, 1, 1],
       [0, 1, 0, 0]])
>>> a.argsort(axis=1)  # axis=1で横にソート。デフォルトと同じ動作
array([[1, 2, 0, 3],
       [2, 0, 1, 3]])

　ドキュメントによると、デフォルトのaxisは-1のようです。ちょっとめずらしい気がしますね。

値の取り出し

　残念ながら２次元配列なので、値を取り出そうとするとなかなか思い通りにはいきません。

>>> a[a.argsort()]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: index 2 is out of bounds for axis 0 with size 2

　こんなコードを書いた人も多いのでは。

sorted_array = a.sort()
sorted_index = a.argsort()

　どう考えても二回ソートするのは二度手間です。

　そこで、numpy.take_along_axis()関数を使います。ちなみに、これはnumpy 1.15から入った新しい関数です。

numpy.take_along_axis — NumPy v1.15 Manual

　以下のように使えます。

>>> np.take_along_axis(a, a.argsort(axis=1), axis=1)
array([[0, 1, 2, 8],
       [0, 1, 1, 7]])

　便利ですね。

【python】引数のデフォルト値は定義時評価なので注意

はじめに

　pythonでは関数の引数にデフォルト値を設定することができます。この機能を使うと、引数が与えられなかったときの挙動を定義することができ、とても便利です。

>>> def f(x="hoge"):
...     print(x)
... 
>>> f("aiu")
aiu
>>> f(x="aiu")
aiu
>>> f()
hoge

　ただし、デフォルト値にlistやdictなどのmutableな型を使うと容易にハマります。それについて説明します。

ハマる例

　デフォルト値に空のリストを追加したとします。

>>> def g(x, l=[]):
...     l.append(x)
...     return l
... 
>>> g(10)
[10]
>>> g(20)
[10, 20]

　さっそくおかしな感じになっています。lが呼び出しのたびに更新されていないようです。そのせいで、前の値が残ってしまうようです。

検証

　listを継承したクラスを作り、コンストラクタが呼ばれるとメッセージを表示するようにします。

>>> class mylist(list):
...     def __init__(self, *args, **kargs):
...         print("__init__ of mylist")
...         list.__init__(self, *args, **kargs)
... 
>>> l = mylist()
__init__ of mylist
>>> l
[]
>>> l.append(1)
>>> l
[1]

　__init__にメッセージを挟んだ以外は何もいじっていないので、元のlistと同様に扱えます。

　このインスタンスをデフォルト値にしてみます。

>>> def h(x, l=mylist()):
...     l.append(x)
...     return l
... 
__init__ of mylist
>>> h(10)
[10]
>>> h(20)
[10, 20]

　なんということでしょう、defが実行されたタイミングで評価されて、それっきりです。

解説と対策

　これは、実はドキュメントにも説明がある事項です。

重要な警告: デフォルト値は 1 度だけしか評価されません。デフォルト値がリストや辞書のような変更可能なオブジェクトの時にはその影響がでます。

4. その他の制御フローツール — Python 3.6.5 ドキュメント

　要するに「仕様」ということですね。

　対策は、これも上のページに書いてあることですが、デフォルト値はNoneにしておき、都度関数の中でインスタンス化することです。

>>> def i(x, l=None):  # ここは定義時評価
...     if l is None:  # ここは実行時評価
...         l = []
...     l.append(x)
...     return l
... 
>>> i(10)
[10]
>>> i(20)
[20]
>>> i(10, l=[-10, 0])
[-10, 0, 10]

　わかってしまえば簡単な話ですが、気づかないと大変なミスをやらかすので注意してください。

　逆に、これを利用してクロージャもどきとして使ったり、難読コードを書くことも可能です（やらないけど）。

まとめ

　再帰的にリストを処理するような関数を書くとき、def f(l=[]):のように書きたくなってしまいますが、バグを産むので注意してください。