blog

リア向け。

Python(Numpy)での多重ループの書き方

Pythonで大量のデータを使おうとすると遅い。
Numpyを使ってもforで回してたら意味ない。遅い。
numpyのmeshgridを書けばfor文いらない。速い。

試しに3重ループになる行列積を書いてみる。

#!/usr/bin/python2
#vim:fileencoding=utf8

import numpy as np

def main():
    arr34 = np.random.randint(0,100,12).reshape((3,4))
    arr45 = np.random.randint(0,100,20).reshape((4,5))

    # 行列の掛け算(行列の場合のみ使えるが一般的な多重ループには無理)
    ref1 =  np.mat(arr34)*np.mat(arr45)

    # 典型的なダメな例。forループを使う。
    ref2 = np.zeros_like(ref1)
    for i in range(arr34.shape[0]):
        for j in range(arr45.shape[1]):
            for k in range(arr45.shape[0]):
                ref2[i,j] += arr34[i,k]*arr45[k,j]

    # これもダメな例。forループを使う。
    from itertools import product
    ref3 = np.zeros_like(ref1)
    for i,j in product(range(arr34.shape[0]),range(arr45.shape[1])):
        ref3[i,j] = np.sum(arr34[i,:]*arr45[:,j])

    # 多重ループの書き方としてはこれがベストか?ただし状況によってはメモリの使用量が心配
    idx345 = np.meshgrid(np.arange(3),np.arange(4),np.arange(5))
    tmp = np.sum(arr34[idx345[0],idx345[1]]*arr45[idx345[1],idx345[2]],axis=0)

    assert np.all(ref1==ref2)
    assert np.all(ref1==ref3)
    assert np.all(tmp==ref1)

if __name__ == "__main__":
    main()

Latexで書いた英文をWordで文法チェックする (Pandoc篇)

Latexはプレーンテキストで書けて、数式も簡単に書けて、構造も把握しやすく、Windows/Mac/Linuxと環境を選ばず、無料で使えて、仕上がりも綺麗という便利ツールですが、Wordと違いスペルチェックと文章校正という機能がありません。

英文チェックだけはWordを使いたいってことで、rtfやhtmlに一度変換して、rtf/htmlからWordに変換するという方法が
spatiohack: LaTeX原稿のスペルチェック
Maraigue風。:LaTeXで書いた文章をWordでスペルチェックする - TeX4htを使う方法
などで紹介されています。

pdfファイルにしたものをhtmlに変換して、同じようにhtmlからwordに変換するという方法もあります。
pdf4htmlEX coolwanglu/pdf2htmlEX · GitHub
というツールが非常に再現度が高いのですが、出力されるhtmlは若干トリッキーです。というか理解できん。。一応ブラウザでhtmlを表示させてテキストをWordにコピペすることは可能でした。

その他にもWordを使わずLatexソースファイルのままチェックする方法として
Aspell - TeX Wiki
もあります。しかしグラマーチェックはしてくれないので、最終的にはWordに頼る人も多いはず。また、複数の英文チェッカーを使いたいという方もいると思います。

Webサービスとして使える英文チェッカーには
Ginger 英文チェッカー
1Checker - Proofreading, Grammar Check, and Smart Enrichment
なんてものがありますが、この場合もLatexソースだとLatexコマンドが邪魔になるので、やはりWordファイルに変換してから使いたいところです。

前置詞とか単語の使い方の使用方法が不安だという箇所がわかっているなら、文法チェッカーを使わずに使用例を探してみるのがいいです。共起表現検索ができるWebサービスが便利です。
WebLSD検索語入力
英語共起表現 - Weblio辞書
Incremental PubMed/MEDLINE Expression Search: inMeXes
昔はNativeCheckerというサイトもあったのですが、残念ながら使えなくなりましたね。

単純に人に聞いてみるというのもいいです。Twitterで答えてくれる人を確保しておくとか、英語の得意な先生や同僚に聞くとか。Lang-8という外国語でブログを書くと、添削してもらえるサイトもあります。
Multi-lingual language learning and language exchange | Lang-8: For learning foreign languages
Lang-8を使う場合は、日本語を勉強中の方の投稿を見たら逆に添削してあげましょう。


話を戻します。Latexで書いた英文をWordでチェックするにはLatexソースファイルをdocx形式にする必要があります。もちろん上で紹介したrtf/html/pdfを経由して変換する方法もあります。しかし、直接Latexソースファイルからdocx形式に変換するツールがあります。
Pandoc - About pandoc
Pandoc ユーザーズガイド 日本語版 - Japanese Pandoc User's Association

このPandocというツール、様々な形式の変換をサポートしていてLatexからdocx形式にも変換できます。インストール方法はこちらがわかりやすい。
HTML - 多様なフォーマットに対応!ドキュメント変換ツールPandocを知ろう - Qiita

Pandocをインストールしたら、次のコマンドでLatexからWord形式に変換できます。

$ pandoc document.tex -s -o document.docx

これでWord形式のファイルができました。あとはWordのスペルチェッカーを使うだけですね。ただし、Latexのeqnarrayやalign環境で書いた数式は正しく再現してくれないようです。$$で書いた文中の数式は変換できていました。eqnarray/align内の英文をチェックして欲しいという人はいないので、問題ないでしょう。

Latexで書いた文ををスペルチェックする方法としてオススメな選択肢の1つです。Latexが嫌になってしまった人にもMarkdownやreStructuredTextで書いてLatexに変換させることができるのでぜひ試して欲しいですね。

ssh接続先で作った図をターミナル上に表示

困っていること。

皆さんsshでサーバとかに入って、データの解析して
図を作って、その図を見るのにどうしてるでしょうか?
自分のPCにscpとかでダウンロードしますか?
その転送作業、面倒くさくありませんか?

端末。sshクライアントを変えよう。

通常、Windowsからのssh接続では
TeraTerm, Putty, Poderosa等を使います。この辺のsshクライアントがメジャーです。
最近私が使い始めたのはRLoginというクライアントですね。
http://nanno.dip.jp/softlib/man/rlogin/
タブ型クライアントとしてはPoderosaがありましたが、もう開発が止まっていて
乗り換えたいと考えてました。そこでRLoginです。詳しい説明はリンク先で。

で、どうやって図を表示させるの?

このRLoginですが、さらさらっと機能一覧を見ていると、なんと
Sixel Graphicsに対応しているではありませんか。
Sixelというのはターミナル上で画像とかを表示させるための規格です。
つまり画像をターミナル上に表示させることができるのです!

例えばmadoka.jpgとかいうファイルがあれば

$ jpegtopnm madoka.jpg| pnmquant 256 | ppmtosixel

と一度pnm/ppm形式に変換して、さらにsixelへ変換させれば
ターミナル上に画像が表示されるわけです。
簡単ですね。これだけです。
ssh接続先のサーバでグラフの図をjpgでもpngでもいいのでファイルを作ってやれば
jpegtopnmとかpngtopnmで図が出せるというわけです。

応用編

私は普段PythonでデータをゴニョゴニョしてPythonのmatplotlibでそのまま図を作ってしまいます。もちろんmatplotlibはpng/eps/pdfとか主要な形式はサポートしているので、ファイルに出力させてしまえばOKです。

が、すべてファイルとして保存しておきたいわけでもないので、Pythonスクリプトから直接sixelを吐けないかと考えるわけです。もちろんできます。

Pysixelというライブラリが公開されているようなので、これを使わさせて頂きます。
https://pypi.python.org/pypi/PySixel
インストール方法はリンク先に載ってるので省略。
PythonでのCode Exampleもリンク先に載ってるので使い方はすぐにわかりますね。
特に難しくもないのでライブラリの紹介だけですべてが終わってしまいそうです。
それではこれを書いてる意味がほとんどなくなるのでmatplotlibと連携する方法を書きましょう。そうしましょう。

まず、Windowとかは表示させないので、matplotlibのbackendをAggに変えてしまいます。

import matplotlib
matplotlib.rcParams['backend'] = "Agg"

そして何らかの図をplotします。

import pylab as pl
import numpy as np
x = np.linspace(0,2*np.pi)
y = np.sin(x)
pl.plot(x,y)

次がポイントです。グラフの図を保存するときはsavefigでOKなのですが、今回は図を保存せずに直接sixel形式として標準出力に出します。
Pysixelではファイルっぽいものを読んで表示させるので、そのファイルっぽいものを作ってそこにデータを流します。

from StringIO import StringIO
buf = StringIO()  # ファイルっぽいもの
pl.savefig(buf)   # bufに保存。ファイルじゃないよ!
buf.seek(0)       # ここポイント。無いとダメ。

seek(0)しないと図の形式がわかりませんとかいうエラーがでます。
あとはPysixelで表示させるだけですよ。簡単。

import sixel
writer = sixel.SixelWriter()
writer.draw(buf)

ただ問題として、このSixel Graphicsの表示機能とscreenを両立させることができないんですよね。。。
2013/11/20追記
kefir_/saitohaさんから解決方法を教えてもらいました。ありがたや。
> https://gist.github.com/saitoha/7546579 今週末にでも改良版を載せられたらいいのですが。。


使い方をまとめた全ソース貼っときます。

#!/usr/bin/env python2
#vim:fileencoding=utf8
import matplotlib
matplotlib.rcParams['backend'] = "Agg"
 
from StringIO import StringIO
import pylab as pl
import numpy as np
import sixel
 
def main():
    pl.figure(figsize=(6,4))
    x = np.linspace(0,2*np.pi)
    y = np.sin(x)
    pl.plot(x,y)
 
    buf = StringIO() # ファイルっぽいもの
    pl.savefig(buf) # bufに保存。ファイルじゃないよ!
    buf.seek(0) # ここポイント。無いとダメ。
 
    writer = sixel.SixelWriter()
    writer.draw(buf)
 
if __name__=="__main__":
    main()

スクリーンショット
f:id:nasing-i:20131114234956p:plain

ターミナル下部にジョブの実行状況を表示

(最初python側でsleepさせてたんですが、この環境だとうまく動かないらしい)
screenのhardstatusで表示させます。
screenのbacktickはコマンドの実行結果を取得することができるので、
ジョブ状況を取得するスクリプト書いてbacktickで読みます。

某環境向けジョブ状況取得スクリプト(使い回しのコピペ作成です)
~/.screen/llqmonitor.pyとでも保存してください
usernameは自分のものに書き換えること。

#!/usr/bin/env python
import subprocess
username = "your username"

qstat = subprocess.Popen(['llq'],stdout=subprocess.PIPE)
grep  = subprocess.Popen(['grep',username],
                         stdin =qstat.stdout,
                         stdout=subprocess.PIPE)
status = [ line.split() for line in grep.stdout.readlines() ]
#jobid,user,day,time,status,priority,category
running = [ row for row in status if row[4]=="R" ] 
waiting  = [ row for row in status if row[4]=="I" ]

if len(running)>0:
    running_names = "(%s)"%(",".join([row[0] for row in running]))
else:
    running_names = ""
print "%d:rest %d%s:running %d:waiting"%(
    len(status),len(running),running_names,len(waiting))

こちらを.screenrcに書きます

backtick 1 30 30 python $HOME/.screen/llqmonitor.py # interval 30 sec
hardstatus alwayslastline
hardstatus string '%{kC}Host:%{kC}%H %-w %{ky}%n %t%{kC} %+w %= %1` %Y/%m/%d %{kC}%C%a'

screenを起動すればジョブ状況が常に表示されるはずです。

こんなかんじ
f:id:nasing-i:20131104174101p:plain

Fortranで出したバイナリデータをC++で読む方法

たぶん(自分を含めて)誰も使わない。


Fortran code (output作成例)

program main
    integer,parameter :: N=256
    integer :: a=1,b=2,i
    real(8) :: arr(N),val=3.14d0

    forall(i=1:N) arr(i)=1.*i

    open(20,file="test.bin",form="unformatted")
    write(20) a,arr,b,val
    close(20)
end program


C++ (要Boost, C++11)

#define BOOST_RESULT_OF_USE_DECLTYPE
#define BOOST_SPIRIT_USE_PHOENIX_V3
#include<iostream>
#include<iterator>
#include<fstream>
#include<string>
#include<cassert>
#include<boost/spirit/include/qi.hpp>
#include<boost/spirit/include/qi_repeat.hpp>
#include<boost/spirit/include/support_multi_pass.hpp>
#include<boost/fusion/include/vector.hpp>

namespace spirit = boost::spirit;
namespace qi     = boost::spirit::qi;
namespace fusion = boost::fusion;

int main() {
    const std::string filename = "test.bin";
    const int N=256;
    fusion::vector<int,
        int, std::vector<double>, int, double,
        int> store;
    fusion::at_c<2>(store).reserve(N);

    std::ifstream ifs(filename, std::ios::binary);
    using iterator_type   = std::istreambuf_iterator<char>;
    auto first = spirit::make_default_multi_pass(iterator_type(ifs)),
         last  = spirit::make_default_multi_pass(iterator_type());
    const auto ret
        = qi::parse(first,last,
                   qi::dword
                >> qi::dword >> qi::repeat(N)[qi::bin_double] >> qi::dword >> qi::bin_double
                >> qi::dword
                 , store);
    assert(ret && fusion::front(store)==fusion::back(store));

    for(auto v : fusion::at_c<2>(store)) {
        std::cout << v << ", ";     // 1, 2, 3, 4, 5, .. , 256
    }
    std::cout << fusion::at_c<4>(store) << std::endl;   // 3.14

    return 0;
}