Python 40行でbash/zshのヒストリファイルをマージするスクリプト

自分はzshユーザなのですが、色々なマシンで作業していると、各マシンのヒストリファイル.zsh-historyをマージしたくなることってありますよね。「あれ、あのコマンド、前に使ったはずなのに、コマンドヒストリにない...」と思ったら、他のマシンで作業した時に入力したコマンドだった、というような場合です。

Pythonのgeneratorやiterator周りの練習として、zshのコマンドヒストリファイル(.zsh-historyや.bash_history)をマージするmergehist.pyを、短めに40行で書いてみました。多分、書ける人が書けば、もっとずっと短く書けると思います。

コマンドヒストリファイルの構造は単純で、例えば、
: 1443359084:0;ls
といったように、基本、各行は
: UNIXタイムスタンプ:0;打ったコマンド
の羅列です。が、必ずUNIXタイムスタンプの昇順に並んでいなければなりません。また、末尾に、バックスラッシュを入れて複数行のコマンドを入力した場合は、「打ったコマンド」の部分に複数行がそのまま記録されますから、その辺も処理してやる必要があります。

使い方は、
python mergehist.py ヒストリファイル1 ヒストリファイル2 ...
という具合に、マージしたいヒストリファイルを複数(3個以上でもよい)、ちゃんとUNIXタイムスタンプの昇順などを考慮してマージされたものが標準出力に出てくる...はずです。

#!/usr/bin/env python
import sys

def histkey(line): return int(line[2:12])

def concat_multiline(it):
  it = iter(it)
  toyield = ''
  for line in it:
    toyield += line 
    if line.strip().endswith('\\'):
      continue
    yield toyield
    toyield = ''

def merger(*iterables, **kwds):
  iters    = map(iter, iterables)
  tonext   = [1] * len(iters)
  stack    = []
  prevpopi = -1

  mykeyf   = lambda x: kwds.get('key')(x[0]) if kwds.get('key') else None
  while True:
    for i, it in enumerate(iters):
      try:
        if tonext[i] > 0:
          stack.append((it.next(),i))
          tonext[i] = 0
      except StopIteration:
        tonext[i] = 0
    stack.sort(key = mykeyf)
    if stack:
      x, prevpopi = stack.pop(0)
      tonext[prevpopi] = 1
      yield x
    else:
      raise StopIteration()

sys.stdout.writelines(l for l in merger(*[concat_multiline(line for line in \
    open(sys.argv[i], r'r')) for i in range(1,len(sys.argv))], key=histkey))

word2vecリンク集

word2vecに関する資料へのリンク集です.全く網羅的でないですが,とりあえず.

中でやっている計算の資料

中身どうでもよくて使いたいだけなら

博士課程一問一答

博士課程 Advent Calendar 2014 のエントリーです。


Q0. 筆者プロフィール
東大情報系博士課程を13年9月学位取得し修了。
専門は自然言語処理


Q1. 博士課程ってなんですか

今から振り返ってみると、「新人がいきなり一人前の研究者として活動されると危なっかしいから、修了するまでに最低限の経験は積んでおいてね。博士の間は失敗は許してやるからさ。」
という教育としての側面が大きいと思います。
博士課程にいた頃は、「研究者として使い物になる人材を選別するための期間」と思っていましたが、いかに天才を持ってきても、やはり、未経験のうちは失敗はするものだと思います。
(もっとも、優秀な人は、1教えれば10の成果を出すのも事実なんですが)




Q2. なんで博士課程に進学するんですか

「したんですか」ですね。
自分が技術者として成功する確率より、研究者として成功する確率の方が高かそう、と思ったからですね。
卒論が、配属初日に通学電車内で30分ぐらいで思いついた研究テーマを持って行ったら、卒論の指導教員の先生(博士の指導教員とは違います)に「それは、確かにまだないね」と言われて、そのまま研究テーマにOKが出て、
その研究テーマで、卒論、国内学会デビュー、分野の2nd Tier国際会議、初ジャーナル、までトントン拍子にいったのが大きいですね。
この研究は、特に賞とか取ってないので評価は、「まぁまぁ」だとは思いましたが、「自分には、どうやら、競争相手のいない研究テーマを見つけ出す直感みたいなものはあるようだ」、という自信がつきました。
結果だけ見ると、最初に30分で思いついたアイデアがトントン拍子に形になったので(今考えると、偶々類似研究がなかっただけかもしれませんけど)、同じレベルの応用ならいくらでも思いつける、と思ったことがありますね。



Q3. 普段どんな生活してんのさ

これも「してたんですか」ですね。
昼夜逆転したりとか、色々生活が乱れたこともありましたが、そういう生活の乱れはタイムロスも大きいですので、週の総労働時間で見ると、普通に9時−5時ペースで働くのが良いと思います。

Q4. 在学時の金銭問題に関して一言

私個人について言えば、一人っ子+実家という家庭事情により、金銭問題はなかったです。
「金銭問題がない」という家庭事情が、博士課程では大きなアドバンテージになる事も見込んで、進学しました。
(逆に、そういう立場の人間が、扶養者の収入が関係する奨学金をもらうのはアンフェアだと思ったので、親にいくら言われても奨学金は申請せず、親に払わせました。親の収入からすれば、私の学費など大したものではないことは分かっていたからです。)
僕が博士に進学する頃、東大は、博士の学生全員に、学費と同額の給付金を出して、実質学費無料化しようという計画がありました。
結果的に無料にはならなかったんですけど、学費半額分の給付金はもらえたので、学振を取るまではそれを使ってました(この給付金は本人の業績評価で扶養者の収入関係なくもらえるので、アンフェアではないと思います)。
実質、学費は年間30万円になるので、この半額制度も大変ありがたい制度だと思います。

Q5. 学振について
学振DC2を博士3年の後半(10月入学だったので)4月から、2年間もらいました。
邦文ジャーナル1本、ACM Transaction1本(内1本は当時採録決定)、トップカンファレンス1本、2nd Tier1本、国内の若手の会の奨励賞が一つ、IPSJの研究会の学生奨励賞が1つ。
お金のことよりも、ようやく、自分の研究が認められたのだ、という感慨が大きかったです。
誰もいないMTGルームで、ネットから結果をこっそり見て、喜びの叫びをあげたら、廊下を通りがかった人が何事かと見に来たのを思い出します。


Q6. インターンについて

楽天技術研究所NYに行きました。
結果的に、博士の業績には繋がりませんでしたが、大変楽しかったです。

情報系なら、Microsoft Research系の研究所インターンには、ぜひ行くべきだと思います。というか、色々あってMSR系のインターンに行けなかったのが博士課程の心残りです。



Q7. 楽しいこと3つ

1) 誰からも認められる確実な「世界初」の安心感
2) 役に立たなくても新しければ良い
3) 自分の名前が残る

研究以外の世界でよくある「世界初」は、「世界初(自社調べ)」、つまり、「自分の所属組織が世界初だと評価した」というだけの事も多いです。
でも、研究では、第三者から見ても確実に「世界初」なものだけが査読に通ります。従って、良くも悪くも自分に対する評価を安心して信じることが出来ます。だからこそ、認められた時の感慨も大きいのだと思います。

「役に立つ」ということと、「新しい」事は全く違います。役に立つまでは膨大な作りこみが必要ですが、その中で新しいことはほとんど必要ない…という事が、往々にしてあります。
ベンチャーにアルバイト的にちょっと働いていたこともあるのですが、役に立つほうが、新しいことをするよりもずっと大変だと感じました。
それも博士課程に進学した理由の一つですね。

後、自分の名前が残るのも大きいです。これにどれだけ意義を感じられるのかは、人によると思います。


Q8. つらいこと3つ

1) 不安定な未来
2) 業績が修了要件になかなか達さない
3) 博論執筆そのもの

Q9. 現時点で後悔していること

やはり、 色々あってMSR系のインターンに行けなかったのが博士課程の心残りです。

Q10. 博士課程に向いている人

「Q7.楽しいこと3つ」に書いたことに意義があると感じられる人が良いのではないでしょうか。

Cのrand()よりmt19937の方が速いことがあるという話

おはようございます。2年ぶりの記事ですね。
もう1月程前になってしまいましたが、id:sleepy_yoshi:20130720 で id:sleepy_yoshi さんが高速な非復元抽出をやっておられ、その中で、Cのrand関数を使っておられました。僕は、普段、std::mt19937を使っていたので、ちょっと比較してみた、という記事です。

C++11では、大別して、2つの擬似乱数生成の方法があります。1つはC(cstdlib)のrand関数で、高速ですが乱数の質が低く、もう1つはrandomヘッダのmt19937(メルセンヌ・ツイスタ)で、低速ですが乱数の質が高い(科学実験に適する)と、一般には思われていると思います。この高速・低速ですが、mt19937を使うことがボトルネックになるほど遅いことは殆どない、というのが今までの実感でした。なので、僕は、非復元抽出のような処理では、特にボトルネックにならない限り、mt19937の方を用いてきました。擬似乱数の質に起因する数値的なバグが万が一出たりすると、デバッグが困難、というのがその理由です。

が、ちゃんと速度比較をしたことがなかったので、比較をしてみました。すると…定説に反して、Cのrand()よりmt19937の方が、使い方によっては速い、という結果になりました。

「使い方」というのはこういう事です。今、[0,1)の範囲の実数の乱数が得られるとします。これをそのまま使う事もありますが、よくあるのは、ランダムに要素を選ぶために、[0,N)のような範囲の整数の乱数が欲しい、という場合です。このような整数の乱数を、実数の乱数から得る時によく使われているテクニックとしては、N * [0,1)の範囲の実数、を計算した上、整数型に変換して、[0,N)の範囲の乱数を求める、という方法があります。一方、mt19937を使う場合は、C++11ではそのような変換は必要なく、std::uniform_int_distributionという専用の関数があります。
比較してみました。

ソースコード(github)

#include <iostream>
#include <random>
#include <chrono>
#include <cstdlib>

namespace{
  class mtrand_t{
    std::mt19937 mt_;
    public:
    inline std::string name() const{ return "mt19937"; }
    inline double      zeroonerand() const { std::uniform_real_distribution<double> gen(0.0, 1.0); return gen(mt_);} 
    inline size_t      nrand(size_t n_rand) const { std::uniform_int_distribution<size_t> gen(0, n_rand-1); return gen(mt_);} 
    inline void        seed(int seed) { mt_.seed(seed); } 
  };
  class crand_t{
    public:
    inline std::string name() const{ return "std::rand"; }
    inline double      zeroonerand() const { return (((double) std::rand()-1)/((double)RAND_MAX)); } 
    inline size_t      nrand(size_t n_rand) const{ return n_rand * zeroonerand();}
    inline void        seed(int seed) { std::srand(seed); }
  };

  class zeroone_caller_t{
    public:
    template<class P>
    inline void call(const P randgen) const{ randgen->zeroonerand( ); }
  };
  class nrand_caller_t{
    size_t n_rand_;
    public:
    template<class P>
    inline void call(const P randgen) const{ randgen->nrand( n_rand_); }
    inline size_t n_rand()const{ return n_rand_;}
    inline void set_n_rand(size_t n_rand){ n_rand_ = n_rand; }
    nrand_caller_t(size_t n_rand): n_rand_(n_rand) {}
  };
  template<class T, class CALLER>
  void measure_time(T* randgen, const CALLER& caller, size_t n_iter, int seed){
    randgen->seed(seed);
    auto start = std::chrono::system_clock::now();
    for(size_t i = 0 ; i < n_iter ; ++i){  caller.call(randgen);  }
    auto d = std::chrono::system_clock::now() - start;
    std::cout << randgen->name() << ":\t" << std::chrono::duration<double>(d).count() << std::endl;
  }
}// namespace

int main(int argc, char* argv[]){
  size_t n_iter = 1000*1000*1000;
  size_t n_rand = 1000;
  zeroone_caller_t zeroone_caller;
  crand_t crand;
  measure_time(&crand, zeroone_caller, n_iter,  7);
  mtrand_t mtrand;
  measure_time(&mtrand, zeroone_caller, n_iter, 7);

  nrand_caller_t nrand_caller(n_rand);
  measure_time(&crand, nrand_caller, n_iter,  7);
  measure_time(&mtrand, nrand_caller, n_iter, 7);
}

実験結果:単位(秒)

std::rand: 7.17306
mt19937: 8.87919
std::rand: 7.18293
mt19937: 4.70852

ソースコードでは、1セットにつき、1000^3=1G個=10億個の乱数を生成するのに必要な時間を測っています。std::randがCのrand関数、mt19937が、std::19937です。最初の2行が、[0,1)の範囲の実数の乱数生成の実行時間、次の2行が、[0,N)の範囲の整数の乱数生成の実行時間です。[0,N)範囲の整数の乱数生成の実行時間は、mt19937の方が約2.5秒も短い事が分かりました。

もちろん、実装・環境依存だと思います。が、いずれにしても数秒/1G回程度の速度差を理由に、mt19937の使用を回避しなければならない場合は殆どないはずです。従って、殆どの場合で、mt19937の使用が勧められるという結論になりました。気になる人は、ソースコードコンパイルして実行して、自分の環境でどうなのか、試してみて下さい。乱数の質は、mt19937の方が明らかに高いです。また、最近では、SIMD命令を使った高速化版のSFMTも出ています。

なお、rand()は元々[0,RAND_MAX]の整数を返してくるのに、(((double) std::rand()-1)/*1;のように、わざわざ割り算を使ってdoubleに変換いるから遅くなるのだ、整数のままmodulo演算(%)を使って、rand()%Nとすればいいじゃないか、という人がいるかも知れません。が、その場合、modulo biasというのが出ます。仮に、RAND_MAX=4、N=3だとすると、RAND_MAXがNで割り切れないので、rand()%3では、0,1は2/5の確率で、2は1/5の確率で出る事になり、0,1,2が等確率で出なくなってしまいます。一方、同じrand()を使う方法でも、挙げたコードのように、割り算を使って一回doubleに変換すれば、そのような事にはなりません。

2013/08/19 1:14追記

rand()を使って[0,1)の範囲の乱数を作るのには、もっと深い話があるようです id:takeda25 さんが関連記事を書いてくださりました。
いつからその方法で偏りのない乱数が得られると錯覚していた? - アスペ日記

また、id:yatt さん他何人かの方から指摘されましたが、そもそも、mt19937のページに、「Cの標準ライブラリのrand()より高速なこともあります。」と書かれているようです。

*1:double)RAND_MAX

Windows 7でnumpy, scipyインストール

明日は,Tokyo.SciPyですね.Windows 7にnumpy, scipyをインストールしましょう.

  1. http://www.python.org/download/Python 2.7.2 Windows Installerをダウンロードしてインストール.僕の環境は64bitで64bitと書いてあるインストーラもありますが,こちらを使ってみました.numpy, scipyのバイナリはwin32と書かれたものしかないので,何かエラーが出そうな予感がしたからです.
  2. http://sourceforge.net/projects/numpy/files/NumPy/1.6.1/で,numpy-1.6.1-win32-superpack-python2.7.exeをダウンロードしてインストール.
  3. http://sourceforge.net/projects/scipy/files/scipy/0.9.0/で,scipy-0.9.0-win32-superpack-python2.7.exeをダウンロードしてインストール.

これだけです.
基本的には,こちらの記事を参考にしました.

余談:
Windowsでは,cygwinのパッケージとしてpythonとnumpyをインストールすることは簡単に出来ますが,scipyのパッケージはcygwinにはないようです.そこで,最初,cygwinの方のpythonにscipyをインストールしようとしたところ,次のように苦労した挙句失敗しました:
scipyが裏で使っているBLASコンパイルから始めなければならない→BLASはインターフェースの名前なので,ATLASとかGotoBLASとかたくさんある→とりあえず,gcc4とgfortran4を入れて,GotoBLASをコンパイルしようとしてコンパイルにこける.

このように,本気を出すならBLASコンパイルからcygwinmingwのgcc4を使ってやった方が良いのでしょうが,そこからscipyのインストールまでやろうとすると,一日仕事になってしまうので,用意されているバイナリをそのままインストールするだけで手を打つのが楽なようです.

snappy_unittestをやってみたよ

DSIRNLPの発表にあったsnappyは全く知らなかったので,大変参考になりました.かなり速いマシンで,snappy_unittestをやってみました.

コンパイルでちょっとはまったのは,GTEST_CONFIG環境変数に,google testのgtest_configへのパスを設定してからconfigureしないとコンパイルが通りませんでした.
また,以下のgflagsも必要なようです.

google testのインストールは,1.6.0からmake installではなくてちょっと一手間かかるようになっているので(google testのreadmeなどを参照),gtest-1.6.0は入っているものとして,たぶん,これでsnappy_unittestが試せると思います.

tar zvxf gflags-1.5.0.tar.gz
tar zxvf snappy-1.0.3.tar.gz
cd gflags-1.5.0
./configure --prefix=~/usr/local
make
make install
cd ..
cd snappy-1.0.3
GTEST_CONFIG=gtest-1.6.0/scripts/へのパス; ./configure --prefix=~/usr/local
make
make install
./snappy_unittest

./snappy_unittestの結果:

Running microbenchmarks.
WARNING: Compiled with assertions enabled, will be slow.
Benchmark            Time(ns)    CPU(ns) Iterations
                                                                                                    • -
BM_UFlat/0 62529 62448 1489 1.5GB/s html BM_UFlat/1 642063 641690 268 1.0GB/s urls BM_UFlat/2 9670 9661 20491 12.2GB/s jpg BM_UFlat/3 25616 25693 3619 3.4GB/s pdf BM_UFlat/4 252631 252638 467 1.5GB/s html4 BM_UFlat/5 25017 24965 8210 939.8MB/s cp BM_UFlat/6 10766 10768 18570 987.5MB/s c BM_UFlat/7 3401 3393 23866 1.0GB/s lsp BM_UFlat/8 959264 960528 178 1022.4MB/s xls BM_UFlat/9 197300 197522 572 734.3MB/s txt1 BM_UFlat/10 170729 171082 713 697.8MB/s txt2 BM_UFlat/11 523052 523450 340 777.5MB/s txt3 BM_UFlat/12 703270 706000 262 650.9MB/s txt4 BM_UFlat/13 301134 300784 482 1.6GB/s bin BM_UFlat/14 40605 40711 2284 895.8MB/s sum BM_UFlat/15 4486 4477 20768 900.3MB/s man BM_UFlat/16 66576 66546 3035 1.7GB/s pb BM_UFlat/17 200036 198865 543 883.9MB/s gaviota BM_UValidate/0 35930 35927 5538 2.7GB/s html BM_UValidate/1 407379 408102 490 1.6GB/s urls BM_UValidate/2 188 188 833333 627.7GB/s jpg BM_UValidate/3 12332 12298 15772 7.1GB/s pdf BM_UValidate/4 144234 144069 1388 2.6GB/s html4 BM_ZFlat/0 258448 258402 770 377.9MB/s html (23.57 %) BM_ZFlat/1 2720520 2719580 100 246.2MB/s urls (50.89 %) BM_ZFlat/2 30047 29972 6605 3.9GB/s jpg (99.88 %) BM_ZFlat/3 103862 103975 1904 865.2MB/s pdf (82.13 %) BM_ZFlat/4 1017405 1020358 195 382.8MB/s html4 (23.55 %) BM_ZFlat/5 105378 105219 1891 223.0MB/s cp (48.12 %) BM_ZFlat/6 48175 48191 4108 220.7MB/s c (42.40 %) BM_ZFlat/7 15026 15018 12383 236.3MB/s lsp (48.37 %) BM_ZFlat/8 3555570 3559460 100 275.9MB/s xls (41.34 %) BM_ZFlat/9 903509 904404 220 160.4MB/s txt1 (59.81 %) BM_ZFlat/10 772980 772084 259 154.6MB/s txt2 (64.07 %) BM_ZFlat/11 2399410 2399640 100 169.6MB/s txt3 (57.11 %) BM_ZFlat/12 3173020 3169520 100 145.0MB/s txt4 (68.35 %) BM_ZFlat/13 1005974 1004898 198 487.1MB/s bin (18.21 %) BM_ZFlat/14 166759 167058 1197 218.3MB/s sum (51.88 %) BM_ZFlat/15 20407 20352 9285 198.1MB/s man (59.36 %) BM_ZFlat/16 262335 262147 759 431.4MB/s pb (23.15 %) BM_ZFlat/17 834067 835316 237 210.4MB/s gaviota (38.27 %) Running correctness tests. Crazy decompression lengths not checked on 64-bit build All tests passed.

行列(画像)分解アルゴリズムGoDec (Zhou+, ICML2011)の実装を公開しました

つい2週間ほど前,機械学習のトップカンファレンスICMLが開催されました.その中のGoDecという行列分解アルゴリズムを実装したので公開します.このアルゴリズムは,簡単にいえば「外れ値抜き特異値分解」で,昨日のICML読み会で発表しました.論文はこれです.

厳密な版(Nai:ve GoDec)は遅いですが実装は非常に簡単です.遅い版でも,数百x数百ピクセルの小さな画像であれば,十分実用的な速度で動くので,実装して試してみた次第です.GoDecの論文では,この厳密な版(Nai:ve GoDec)が線形収束することを証明した上で,さらに,実用的に早くなるように(証明はないようですが)乱択化低ランク近似で高速化したFast GoDecも提案されています.僕のソースコードはNew BSDライセンスにしておきましたが,GoDec自体の特許関係は全く調べておりません.もしかしたら,特許が取られているかもしれませんので,商用にGoDecを用いることには注意が必要です.

ソースコードと使い方

コンパイルには,Eigenが必要です.C++0xを使っているので,gcc 4.6.0以上奨励です.

git clone git@github.com:niam/godec.git
cd godec
./waf configure
./waf build
build/src/godec -r 8 -k 762 hogehoge.png

これは,簡単にいえば,「762個までの外れ値を許容して,hogehoge.pngを行列だと思ってランク8で低ランク近似しろ」ということです.

例えば,以下の画像をhogehoge.pngとして入力したとします.

すると,X.png, L,png, S.png, LpS.pngがカレントディレクトリに出来ます.X.pngは,元の画像をグレイスケールに変換したもので,これが,分解対象の行列になります.

これをランク8の行列Lと,要素数762の行列Sに分解し,L+Sを実行したものがLpSです.

GoDecは,特異値分解に加えて外れ値の数kの自由度を持つので,特異値分解と比較するためには,自由度の数を合わせないとフェアではありません.200x180行列において,1ランク分の情報を保存するには200次元の特異ベクトルu+180次元の特異ベクトルv+特異値σの381個の自由度が必要になります.従って,2ランク分の情報がちょうど762自由度なので,単純な特異値分解で上位8+2=10で低ランク近似した場合と自由度がおなじになります.元の画像を特異値分解で上位10で低ランク近似したものが,下の画像です.

GoDecでは,頭の上の木漏れ日による光の点が外れ値として処理されハッキリ表示できているのに対し,特異値分解(SVD)では周りの色に紛れてしまっていることが分かります.

SVDと比較する場合,この自由度を合わせる作業は面倒くさいので,自動的に自由度を合わせてくれるモードも作りました.

build/src/godec -r 10 --rdiff=2 --comparesvdmode hogehoge.png

とやると,上の実験設定と同じ実験設定になります.つまり,特異値分解の方のランクが10で,それと2ランク差がつくようにGoDecのランクと,許容する外れ値の数kを自動的に設定してくれます.comparesvdmodeでは,10ランクで近似したsvd.pngも作られます.

アルゴリズムの解説

行列を,小さいランクの行列で近似することを低ランク近似と言います.低ランク近似の代表的な方法が特異値分解です.しかし,特異値分解は,行列内に,外れ値のような他の要素よりずっと大きい/小さい要素があったりすると,その要素を上手く再現出来なかったり,また,その要素に引きずられて,その要素の周りがにじむように近似誤差が大きくなってしまったりすることが知られています.

このアルゴリズムは,簡単に言えば「外れ値抜きの特異値分解」です.

  1. 外れ値を抜いて特異値分解でランクrで低ランク近似する
  2. もとの行列と低ランク近似した行列を見比べて,要素の差(の絶対値)が大きい「外れ値」要素を差の大きい方からk個見つけてきて新しい外れ値にする

を繰り返したものが,線形収束することを示しているのがこの論文です.

僕の発表資料はこちらです.解釈に間違いがある可能性がありますので,詳細は論文本体を参照されることをお勧めします.