setup.py や easy_install での version 指定

easy_install で PyPI から package/module を取ってくると、引数無しでは最新版が落ちてきます。例えば以下のようにやると、2011 年 12 月現在で v1.6 が install されます。

$ easy_install networkx

古い version を指定して落としたい時には、次のように version 番号を指定すれば OK。これは v1.4 の場合です。

$ easy_install networkx==1.4

setup.py の install_requires で version 指定をする場合も以下のようにやります。この場合、Astropysics が networkx の v1.5 以降に対応していないため、v1.4 を指定しています。

from numpy.distutils.core import setup

setup(name="PyTeVCat",
      version="1.1.1",
      description="Python wrapper for TeVCat",
      author="Akira Okumura",
      author_email="oxon@mac.com",
      url='https://sourceforge.net/p/pytevcat/',
      license='BSD License',
      platforms=['MacOS :: MacOS X', 'POSIX', 'Windows'],
      packages=["tevcat"],
      install_requires=['astropysics', 'networkx==1.4'],
      package_data={"tevcat": ["img/*.png",]},
      classifiers=['Topic :: Scientific/Engineering :: Astronomy',
                   'Topic :: Scientific/Engineering :: Physics',
                   'Development Status :: 4 - Beta',
                   'Programming Language :: Python',
                   ],
      long_description='tevcat.Python interface for TeVCat (http://tevcat.uchicago.edu/)'
      )

VMWare Fusion + AutoCAD でマウスの動きがアホ遅いのを直す方法

以下の構成で、マウスの動きがめちゃくちゃに遅いので、ググってたら直しかたが分かりました。多分、この組み合わせじゃなくても、VMware + AutoCAD ってのさえ同じであれば、同様の症状が出るはずです。

AutoCAD の window で製図の部分に mouse が移動した時だけ、とんでもなく移動速度が遅くなります。非常にストレスが溜まって、もう VMware なんてやめて、Windows 専用機を買いたくなるレベル。でも解決方法が分かりました。

VMware の環境設定で、マウスの最適化を切れば良い。これだけ。ゲームなんてしないので、問題なし。


▲ Gaming: Never optimize mouse for games を選択する。

久しぶりにマンガで声出して笑ったwww『マンガ 量子力学』石川真之介 著

僕の研究仲間のお一人である石川真之介博士が新作『マンガ 量子力学』を出版されました。おめでとうございます。ご本人から献本して頂いたので、久しぶりに書評を。



1. どういう人にお薦めか

マサルさんとか、王様はロバとか、脱力感のある、シュールな笑いが好きな人向け。ただし、若干の量子力学に関する知識を持っていたほうが良いでしょう。量子力学の試験勉強の入門用だとか、物理系じゃない方が量子力学を知るための入門用としては向かないかもしれません。

じゃあ、何のための本だよ、と突っ込みたくなりますが、この本の題は『マンガ 量子力学』であって、「マンガで分かる」とか「マンガで勉強」なんてどこにも謳ってないんです。帯にこそ「『量子力学』の入門書」と書いてありますが、「ギャクマンガ」をブルーバックス*1が出すわけにはいかないので、これはただの引っかけです。

量子力学のネタを随所に練り込んだ、不思議世界の魔法少女ギャグマンガだと思って購入すれば失敗しません。そういう心積もりで手に取るなら、絶対に期待を裏切ることはないはずです。忘年会の景品やクリスマスプレゼントを考えるこの季節、懐に余裕のある理系の学生は、ネタで 1 冊、購入してみるのも良いかも。

2. なぜ『マンガ 量子力学』をあえて出版するのか

これは本人や出版社に聞いてみないと分かりませんが、量子力学の入門書やマンガ参考書が数ある中で、なぜブルーバックスからわざわざ出すのでしょう。この本の最大の特長は、原作と作画が同一人物ということです。著者の石川真之介博士は、人工衛星を使って宇宙を X 線で観測するのが専門の、バリバリの研究者です。もちろん、プロの漫画家ではなく宇宙物理の研究で飯を食っています*2。しかし自分の脳内にある物理の考え方を、自分の趣味とするマンガで世の中に伝えたいという思い、それが彼のこの作品です。

もちろんプロの漫画家ではないので、そこらへんのマンガに比べれば絵は下手くそです。たまにカイジのような関節の捻りがあったり、『進撃の巨人』に出てくるマイナーキャラのような表情を見せることがあります。でも、著者がマンガ好きの物理屋だからこそ、魔法少女量子力学を使って戦闘するなんていう、意味の分からない設定が生れるんです。

例えば、同じような系統のものとして『マンガでわかる量子力学』を見てみましょう。Amazon の「なか見検索」もできます。



こういう本は基本的に、物理屋が考えた教科書をマンガ化しただけです。なのでマンガとして面白くありません*3。どうせマンガを読んで量子力学が分かるはずないんだから、『なっとくする量子力学『なっとくする演習・量子力学』でも読んだほうがいいでしょう*4。『なっとく〜』は初学者に本当にお薦めです。僕も持っています。


話を戻して、『マンガ 量子力学』の最大の売りは、現役の宇宙物理研究者が脳内で作り上げた魔法少女の物語を、自分で絵まで描いて出版してしまったところです。物理の研究者なんて変な奴ばかりなので、自然とギャグマンガになります。著者の意図は多分違うのでしょうが、明らかにギャグマンガです。

3. あらすじ

毎日を退屈に過ごす平凡でちょっと可愛い (?) フェル美 (ふぇるみ) ちゃんが、人間語を操る謎の野良犬フォト丸 (ふぉとまる) の助けを借りて、「量子魔法少女」に変身して大活躍する物語です。古典力学を使う古典 (こてん) ちゃんとの魔法弾を使った壮絶なバトル、そして「先輩」とフェル美ちゃんの恋の行方は…。

いやもう、石川君が博士論文を書きながら、こんなわけの分からない物語を同時並行で描いてたことが凄い。時間の使い方とその発想がやばい。「お巡りさん、ここです!」のレベル。

4. 評価

  • 画力♀ ★★☆☆☆
  • 画力♂ ★☆☆☆☆
  • 画力動物 ★★☆☆☆
  • 独創性 ★★★★★
  • 笑い ★★★★☆
  • 理学書として ★★★☆☆
  • 総合評価 ★★★★☆

ということで個人的にはお薦めです。

5. 改善点

読者として、こうだったらもっと良いのにと思う点を挙げておきます。

まずは画力。これはやっぱりもっと頑張れるけど、あまり本質的でない。今くらいの画力だからこそ、笑って楽しむことができるのかも。あまりに画力がありすぎて、本当に激しい戦闘シーンや先輩との恋愛が描かれていたら、それはこの本の魅力を半減してしまうかもしれない。諸刃の剣。

次に、量子力学についての説明の多くが、スピンの量子化になってしまっていること。前半部分はほとんどこれで物理的内容を占めてしまっている。スピンが分からない人には、理解が追いつかないかもしれないかなと思った。でもスピンを既に知っている人がギャグマンガとして読む本なので、これもまあ良いか。

量子力学の「入門書」を実は目指しているなら、マンガ部分と解説部分のバランスがちょっと悪いかもしれない。後者をもう少し噛み砕いてマンガ部分に詰め込めれば、もっと (物理として) 面白いと思う。多分、アニメ化して半年くらいかけて徐々に物理の内容を説明して行くのが良い。

ブルーバックスは高い。普通のマンガが 400 円とか 500 円で買えるのに、ブルーバックスのこの定価は 820 円。もちろん数千円の理学書に比べれば安いけど、やはりマンガとして 400 円くらいで買えるほうが望ましい。iPad アプリにして、1 冊 100 円、それを 5 冊に分けて内容を充実させて配信、だったりすると良いのにな、と。

*1:ブルーバックスは文庫サイズで気軽に読め、数多くの理系分野の本が出ている、理系少年御用達の講談社のシリーズです。

*2:と言っても、マンガを出版している時点でプロの漫画家か。

*3:読んだことないけど。

*4:初出で書名を間違えてしまいました。Amazon 経由で購入されたと思われる方がいらっしゃるのですが、大変申し訳ありませんでした。

ROOT の Makefile.arch が色々と更新されちゃったので、その対処法

最近の ROOT では、Makefile.arch が 2 点変更になりました。

  1. $ROOTSYS/test/Makefile.arch から $ROOTSYS/etc/Makefile.arch へ移動した
  2. ld の引数に -undefined を使わなくなったので、明示的に使用する library を与えないといけない

そのため、古い ROOT、新しめの ROOT の両方に対応できる Makefile を書くには、いくつか変更を要します。

まず 1 点目。これまで、自分は $ROOTSYS/test/Makefile を真似して Makefile を書いてきました。同じことをしている人は多いでしょう。しかし $ROOTSYS/test/Makefile.arch が $ROOTSYS/etc/Makefile.arch に移動したため、これを Makefile の中で include しては、file が存在しないと怒られてしまいます。Revision 41594 以降でこの問題が現れます。

そこで、従来

include $(ROOTSYS)/test/Makefile.arch

としてきたものを、次のように書き換えました。

MAKEARCH	:=	$(shell find $(ROOTSYS)/test -name Makefile.arch)

ifeq ($(MAKEARCH), )
# 41594 or later
MAKEARCH	:=	$(shell find $(ROOTSYS)/etc -name Makefile.arch)
endif

include $(MAKEARCH)

このようにすることで、実際に存在するほうの Makefile.arch を読んでくれます。もう少し綺麗な書き方があると思いますが、Makefile あんまり詳しくないので、これで良いことにします。

2 点目。Revision 40240 から、OS X 上では ld の -undefined という引数が使われなくなることになりました。そのため、標準的でない ROOT の library や外部 library を使用する時には、ld をする時点で明示的にそれら library を link してやる必要があります。僕の場合は、libGeom と libGeomPainter を使う library を書いています。

$ROOTSYS/test/Makefile を見ると、40240 以降では以下のような書き方に変わっています。

ifeq ($(PLATFORM),macosx)
# We need to make both the .dylib and the .so
		$(LD) $(SOFLAGS)$@ $(LDFLAGS) $^ $(OutPutOpt) $@ $(EXPLLINKLIBS)

それ以前だと、該当箇所は以下の通りです。

ifeq ($(PLATFORM),macosx)
# We need to make both the .dylib and the .so
		$(LD) $(SOFLAGS)$@ $(LDFLAGS) $^ $(OutPutOpt) $@

Makefile.arch の中に $(EXPLLINKLIBS) という変数が追加されており、この中に必要な library が記述されるようになりました。変数の中身は、-L/usr/local/root/lib -lGpad -lHist のようなものです。

40240 より前では -undefined が引数に与えられていたため、このような library の指定がなくとも ld が動きました。しかし、40240 以降を使う場合は、以下のように libGeom などを $(EXPLLINKLIBS) に追加しましょう。

ifneq ($(EXPLLINKLIBS), )
EXPLLINKLIBS	+= -lGeom -lGeomPainter
endif

以上の変更を施した実際の Makefile は、ここで見られます。

銀河系内全天画像を扱う人は CAR Projection をちゃんと理解したほうが良い

世の中に転がっている銀河系内全天画像 (all-sky map) の FITS の header には、WCS (World Coordinate System) が間違って書き込まれているものを散見します。WCS には様々な projection (投影法) が規定されているのですが、その中で CAR projection を勘違いしたまま使用する人を何度も見てきました。そのような間違った WCS が書き込まれた FITS を解析した場合、座標系が間違っているのでとんでもない解析ミスに繋がる可能性があります。

世の中に落ちている HI、CO、Av など、様々な全天画像の FITS header に間違いが含まれています。どこからか落としてきた FITS 画像を使う場合は、使用前に信頼できる描画ソフトで座標系を描いてみることをお勧めします。

WCS とは

WCS は FITS 画像上の pixel がどのように実際の天球面に対応するかを定めた座標変換の規則のことです。通常の CCD 画像などは 2 次元の広がりを持つだけですが、現実の天体は球面座標で表される天球面上に存在しています。そのため、画像上の見かけの 2 次元座標は何らかの座標変換によって天球面上に配置されなければいけません。さらに光学系の歪みが無視できない場合、この歪みも WCS にて補正することが可能です。

詳細は以下の論文に記述されています。

  1. M. R. Calabretta and E. W. Greisen (2002) A&A, 395, 1077-1122
  2. M. R. Calabretta, F. G. Valdes, E. W. Greisen, and S. L. Allen (2004) A&A に投稿予定のまま?
  3. M. R. Calabretta and B. F. Roukema (2007) MNRAS, 381, 865-872

WCS が登場する前までは FITS の 2 次元画像をそのまま直交座標系だと見做す単純な座標変換が使われていました。次のような FITS header を用意すれば、CRPIX1/CRPIX2 を基準点とし、FITS の (x, y) 座標の点を銀経銀緯上で {\ell = (x - \mathrm{CRPIX1})\times\mathrm{CRDELT1} + \mathrm{CRVAL1}} {b= (y - \mathrm{CRPIX2})\times\mathrm{CRDELT2} + \mathrm{CRVAL2}}
で表すことが出来ます。

CTYPE1  = 'GLON'
CRVAL1  =   1.800000000000E+02
CRPIX1  =                    1
CDELT1  =  -5.000000000000E-01
CTYPE2  = 'GLAT'
CRVAL2  =  -0.000000000000E+01
CRPIX2  =                  181                                      
CDELT2  =   5.000000000000E-01

これは非常に単純で明快です。

CAR Projection とは

WCS には CAR projection という座標変換が存在します。この CAR は plate carrée のことで、本来は Cartesian ではありません*1。一見すると、上述の WCS を使わない単純な座標変換に似ているのですが、同じだと勘違いして使用すると大きな間違いを犯します。

文献 [1] には次のように書かれている為、ここだけ読むと全く同じに思えるでしょう。

5.2.3. CAR: Plate carrée
The equator and all meridians are correctly scaled in the plate carrée projection14, whose main virtue is that of simplicity. Its formulae are
x = φ, (83)
y = θ. (84)

ただし、この {\phi}{\theta} は必ずしも赤経赤緯や銀経銀緯に対応しません。{(x, y) = (\mathrm{CRPIX1}, \mathrm{CRPIX2})} に対応する天球面の場所から大円を引き、その大円上にある座標が等 {\theta} 線になります。そして、{\phi/\theta} から {\alpha/\delta}{\ell/b} への変換が再度必要になります。同じく文献 [1] に以下のような説明がさらにあります。

2.8. Comparison with linear coordinate systems
It must be stressed that the coordinate transformation described here differs from the linear transformation defined by Wells et al. (1981) even for some simple projections where at first glance they may appear to be the same. Consider the plate carrée projection defined in Sect. 5.2.3 with φ = x, θ = y and illustrated in Fig. 18. Figure 1 shows that while the transformation from (x,y) to (φ,θ) may be linear (in fact identical), there still remains the nonlinear transformation from (φ,θ) to (α, δ).

{\mathrm{CRPIX1/CRPIX2}} の場所から大円が描かれて等 {\theta} 線になるということは、例えばその場所が {(\ell, b) = (0^\circ, -30^\circ)} だった場合、2 次元画像上で x の位置だけ移動すると、{(180^\circ, -30^\circ)} に辿り着かずに {(180^\circ, +30^\circ)} に辿り着くことを意味します。

通常の CAR projection の用途では、恐らく WCS 登場前の CTYPE1 = 'GLON'、CTYPE2 = 'GLAT' と同じ使い方をしたいだけでしょう。その場合は、次のような header が最小構成です*2。これは、360 × 180 pixel の 2 次元画像で、1 pixel が 1° に相当しています。また画像の中心 (x, y) = (180.5, 90.5) が銀河中心に対応しています。

SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                   16 / number of bits per data pixel
NAXIS   =                    2 / number of data axes
NAXIS1  =                  360 / length of data axis 1
NAXIS2  =                  180 / length of data axis 2
CTYPE1  = 'GLON-CAR'           / axis type
CRVAL1  =   0.000000000000E+00 / longitude
CRPIX1  =                180.5 / ref pixel
CDELT1  =  -1.000000000000E-00 / longitude increment
CTYPE2  = 'GLAT-CAR'           / axis type
CRVAL2  =   0.000000000000E+00 / latitude
CRPIX2  =                 90.5 / ref pixel
CDELT2  =   1.000000000000E-00 / longitude increment
END

この画像を DS9 を使って座標軸を表示してみます*3


▲ CRPIX1/CRPIX2 を (ℓ, b) = (0°, 0°) に持ってきた画像の例。緑の十字が CRPIX1/CRPIX2 の場所。

間違った CAR を使った例

画像の上下中央を銀河面に設定したいにも関わらず、CAR projection を理解せずに使用した場合には座標系が歪みます。


▲ CRPIX1/CRPIX2 を (ℓ, b) = (0°, -30°) に持ってきた画像の例。CRPIX2 = 60.5、CRVAL2 = -30 に変更してある。


▲ CRPIX1/CRPIX2 を (ℓ, b) = (0°, -60°) に持ってきた画像の例。CRPIX2 = 30.5、CRVAL2 = -60 に変更してある。


▲ CRPIX1/CRPIX2 を (ℓ, b) = (0°, -90°) に持ってきた画像の例。CRPIX2 = 0.5、CRVAL2 = -90 に変更してある。

間違った FITS Header に出会った場合の修正方法

全天画像などで FITS 作成者が明らかに CAR projection を間違って使っている場合は、基本的に、CRPIX2 の場所を銀河面上に移動すれば座標系は正しいものに直ります。もしくは、'GLON-CAR'、'GLAT-CAR'、'RA---CAR'、'DEC--CAR' などを 'GLON'、'GLAT'、'RA'、'DEC' に置き換えてしまっても良いでしょう。

もし FITS 作成者と連絡が取れる場合は、FITS header に間違いがあることをちゃんと指摘し修正してもらいましょう。

*1:[1] には Although colloquially referred to as “Cartesian” という表現がされています。

*2:Comment は別に必要ありません

*3:DS9 は全天画像をうまく取り扱えず、画像の大きさによっては軸を画像に重ね書き出来ない場合があります。

Tektronix MSO/DPO 4000 Series を Python で動かす

古い機種ばっか使っていて知らなかったのですが、LAN port を備えた最近のオシロは、めちゃくちゃ簡単に使えるんですね。Baud rate とかも気にする必要なく、telnet 接続するだけで動きます。PySerial よりもずっと簡単。しかも新たに何も install しなくて良い。

Tektronix MSO/DPO 4000 Series の場合、次の Python code で一発です。XX.XX.XX.XX の箇所は実際の IP address (DHCP もしくは固定)、4000 は port 番号でこれが default の設定です。

>>> import telnetlib
>>> tel = telnetlib.Telnet('XX.XX.XX.XX', 4000)
>>> tel.write('*idn?\n')
>>> tel.read_until('\n\r')
'TEKTRONIX,MSO4054B,C010237,CF:91.1CT FV:v1.12 \n\r'

Python の標準 module である telnetlib を使えば、オシロの IP と port 番号を指定して command を送受信するだけです。telnetlib.Telnet.read_until() を使うことで、指定の delimiter が出現するまで読み出します。

こちらは socket module を使用する場合の例。

>>> import socket
>>> soc = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
>>> soc.connect(('XX.XX.XX.XX', 4000))
>>> soc.send('*idn?\n')
>>> soc.recv(100000)
'TEKTRONIX,MSO4054B,C010237,CF:91.1CT FV:v1.12 \n\r'

(2011/11/10 追記) どちらでやっても、波形の data 点数が 5690 以上だと 'CURVE?' が正常に動作しません。VXI-11 だと動くのですが。

(2011/11/11 追記) Tektronix の日本語電話サポートの情報によると、raw socket での通信に以前から bug があり、転送がうまくいかないとのこと。Protocol を terminal mode にした場合は動作し (ただし binary 転送はうまくいかない)、None にした場合は一切応答しないはずとのこと。しかし、自分とこの環境では部分的には動作しているので、この情報と整合しない。

Python 2.6 以降だと PySerial で '\r' が readline の delimiter に使えない

PySerial では Python 2.6 以降の場合には、'\r' (CR) が readline のときの delimiter に使えません。つまり、装置側の出力に '\r' が delimiter として使われている場合、次の code だと '\n' (LF) を待ってしまうために time out してしまいます。

import serial
s = serial.Serial('/dev/ttyUSB0', ... )
s.write('SOME_COMMAND?')
s.readline()

これは、Python 2.6 以降の場合に PySerial が io.RawIOBase.readline を使うようになっているからです。io.RawIOBase は delimiter が '\n' に決め打ちなため、制御する装置側で delimiter の変更ができない場合は詰んじゃいます。

入門 Python 3

入門 Python 3

Python 2.5 以前の場合だと、PySerial は serialutil.FileLike という独自の file IO を使うため、delimiter を自分で設定することができました。

import serial
s = serial.Serial('/dev/ttyUSB0', ... )
s.write('SOME_COMMAND?')
s.readline(eol = '\r')

どうしても '\r' を delimiter に使いたい場合は、単純には serial.Serial.readline を override してしまいます。

class MySerial(serial.Serial):
    """
    Wrapper for Serial
    """
    try:
        import io
    except ImportError:
        # serial.Serial inherits serial.FileLike
        pass
    else:
        def readline(self):
            """
            Overrides io.RawIOBase.readline which cannot handle with '\r' delimiters
            """
            ret = ''
            while True:
                c = self.read(1)
                if c == '':
                    return ret
                elif c == '\r':
                    return ret + c
                else:
                    ret += c

単純に '\r' を読むまで 1 文字ずつ進めているだけですが、実用上は問題ないでしょう。

もう少し io.RawIOBase.readline に近づけたい場合は、それをそのまま改造してしまうのも手でしょう。例えば OS XPython 2.6 の場合、/System/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/io.py の中で次のように io.RawIOBase.readline が定義されているので、'\n' を '\r' に置き換えてしまえば OK のはずです。

    def readline(self, limit = -1):
        r"""Read and return a line from the stream.

        If limit is specified, at most limit bytes will be read.

        The line terminator is always b'\n' for binary files; for text
        files, the newlines argument to open can be used to select the line
        terminator(s) recognized.
        """
        self._checkClosed()
        if hasattr(self, "peek"):
            def nreadahead():
                readahead = self.peek(1)
                if not readahead:
                    return 1
                n = (readahead.find(b"\n") + 1) or len(readahead)
                if limit >= 0:
                    n = min(n, limit)
                return n
        else:
            def nreadahead():
                return 1
        if limit is None:
            limit = -1
        if not isinstance(limit, (int, long)):
            raise TypeError("limit must be an integer")
        res = bytearray()
        while limit < 0 or len(res) < limit:
            b = self.read(nreadahead())
            if not b:
                break
            res += b
            if res.endswith(b"\n"):
                break
        return bytes(res)