2017-04-28

RELION の classification で、クラスの割合を抽出する

RELION の 2D / 3D classification で、iteration ごとに各クラスの割合がどう変化したかをプロットしたいと言われたので、スクリプトを作った。

for i in `seq 0 20`; do 
 j=`printf %03d $i`;
 echo -n $j" ";
 relion_star_printtable run_it${j}_model.star data_model_classes _rlnClassDistribution | tr '\n' ' ';
 echo;
done > dist.dat

relion_star_printtable がクラスごとに改行してしまうので、tr を使って置換するなど、ちょっとダサいことになっている。

こうすると、run_it_0NN_model.star (NN = 00 .. 20) から、割合を切り出して行に並べたデータファイル dist.dat

000 0.200000 0.200000 0.200000 0.200000 0.200000
001 0.200000 0.200000 0.200000 0.200000 0.200000
002 0.199695 0.201143 0.198778 0.200882 0.199502
003 0.214579 0.171917 0.207016 0.180553 0.225935
004 0.228429 0.149543 0.176681 0.134442 0.310905
005 0.105184 0.321839 0.107021 0.308757 0.157199 
(後略)

ができるので、gnuplot で

plot for [i=2:6] 'dist.dat' using 1:i with line title "Class ".(i-1)

などとしてプロットする。for を使うには、比較的新しいバージョンが必要。
4.6 ではできたが、4.2 ではダメだった。

なお、ファイル数 (20) やクラス数 (5) を指定しなくて済むようにするのを、読者の練習課題とする。

2017-03-11

Jupyter notebook で cctbx などを使う

プログラミング

DIALS を bootstrap.py からビルドしたあと、base/bin/pip を使って Jupyter notebook をインストールしても cctbx を import できない。環境変数を設定する wrapper となっている dials.python などと違い、PYTHONPATH や LD_LIBRARY_PATH が設定されていないようだ。

次のように設定して解決した。

jupyter notebook --generate-config

として、設定ファイル ~/.jupyter/jupyter_notebook_config.py を生成し、そこに以下のような記述を付け足す。何を足したらいいかは、dials.python の中で当該環境変数を調べ、Jupyter notebook 上での値と比べれば分かる。

import os
os.environ['PYTHONPATH'] = "/opt/dials/modules/cctbx_project:/opt/dials/modules:/opt/dials/modules/cctbx_project/boost_adaptbx:/opt/dials/modules/cctbx_project/libtbx/pythonpath:/opt/dials/build/lib:" + os.environ['PYTHONPATH']
os.environ['LD_LIBRARY_PATH'] = "/opt/dials/build/lib:/opt/dials/base/lib64:/opt/dials/base/lib:" + os.environ['LD_LIBRARY_PATH']

2016-05-04

指数付けアルゴリズム総説

構造

指数付けの目的は、逆空間の基底ベクトルを決定することである。3 つの基底ベクトル(縦ベクトルである)を横に並べて 3 x 3 行列にしたものを方位行列 A と呼ぶ。成分を書き下すと、
$A = \begin{pmatrix}\mathbf{a^*} & \mathbf{b^*} & \mathbf{c^*} \end{pmatrix} = \begin{pmatrix}a^*_x & b^*_x & c^*_x\\ a^*_y & b^*_y & c^*_y\\ a^*_z & b^*_z & c^*_z \end{pmatrix}$
である。
方位行列は 9 成分を持つから、自由度も 9 である。そのうち 6 つの自由度は逆格子定数(a, b, c, $\alpha$ , $\beta$ , $\gamma$ ) であり、残り 3 つは結晶方位を表す。方位行列を、なんらかの基準方位(例えば a* 軸が実験室系の x 軸に一致し、b* 軸が x-y 平面上にあるとか)に合わせた基底行列と、それを実際の結晶の向きへ回転させる回転行列の積へと分解することもしばしば行われる。たとえば、MOSFLM における A = U B という分解がそれである。B が基底を表し、U が回転を表す。U は回転行列だから、det(U) = 1 である。逆格子定数から具体的に成分を求める方法は Orthogonalization matrix の導出と、格子定数による偏微分 - biochem_fan's note で説明した。

指数付けの結果得られた A 行列によって、反射の指数 $\mathbf{h} = {}^t(h, k, l)$ から、逆空間座標 $\mathbf{r} = {}^t(r_x, r_y, r_z)$ を求めることができる。すなわち、
$\mathbf{r} = \begin{pmatrix}r_x \\ r_y \\ r_z \end{pmatrix} = \begin{pmatrix}\mathbf{a^*} & \mathbf{b^*} & \mathbf{c^*} \end{pmatrix} = A \begin{pmatrix}h \\ k \\ l \end{pmatrix}$
である。逆に、観察された反射に指数をつけるには、
$\mathbf{h} = A^{-1} \mathbf{r}$
を計算すればよい。実際には、 $\mathbf{r}$ には実験誤差があるから、 $\mathbf{h}$ の成分は整数にはならない。最寄りの整数値になるように丸めるわけだが、あまりに小数成分が大きい場合は、(1) 指数付けが間違っている (2) 複数格子が存在する (3) ice ring などを spot finder が拾っている可能性が高い。

ここまでの議論から分かるように、観測された回折点の逆空間座標と指数のペアが 3 つ以上あれば、方位行列は確定する。一枚の回折画像には数十から数百の回折点があるのが普通だから、指数さえ決まれば方位行列は最小二乗法などで精密化できる。しかし、現実には、方位行列が分からないから指数がつかず、指数がつかないから方位行列が定まらないという問題に直面する。したがって、指数付けには、上の関係を超えた考察が必要となる。

続く!

周期を見つける方法。1次元 FFT (DPS), 3 次元 FFT (Denzo ほか)
差ベクトルによる方法。REFIX。
この 2 つが数学的に等価であること。Patterson 図。
方位を優先的に調べる方法。DirAx, DIALS known cell, MOSFLM prior cell
初期精密化。誤差に対して robust であるために。差ベクトル法

などを述べたい。

2016-05-04

方向ベクトルの表現について

数学

長さが 1 の方位ベクトル (x, y, z) を表現したい場合がある。長さが1だから、自由度は 2 である。しかし、x と y 成分を与えても、z の値は南半球と北半球の二通り $\pm\sqrt{1 - x^2 - y^2}$ あって、一意には決まらない。つまり、これは筋が悪い表現方法である。

非線形最小二乗法などのために、変化を局所的に線形化することがある。この場合も、 $(x + \del x, y + \del y, \sqrt{1 - (x + \del x)^2 - (y + \del y)^2}$ のように展開するのはよくない。上述の z 成分の符号の自由度があるからだ。しかし、z 成分も二乗以上の微小量を無視すると、 $(x + \del x, y + \del y, z - \del x - \del y)$ は、実際、元のベクトルと長さが一致しているので、つじつまはあっている。とはいえ、最小二乗法と相性がよいのは、極座標・球面座標系である。

REFIX 論文(W. Kabsch "Automatic indexing of rotation diffraction pattern" J. Appl. Cryst. (1998). 21, 67-72) の式 8 を例に工夫を見てみよう。ここでは、入射ビームが
$\mathbf{s_0} = \begin{pmatrix}-\sin\eta\cos\kappa/\lambda\\\sin\kappa/\lambda\\\cos\eta\cos\kappa/\lambda\end{pmatrix}$
と球面座標系で表現されている。この表現の場合、 $\kappa$ がほぼ 90 度になる状況では、 $\eta$ の値が縮退してしまう。しかし、実験室座標系は回転軸が $\mathbf{e_2}$ 方向を向くように取られているので、回転軸と入射ビームが一致するという実際にはありえない状況にならない限り、特異点の問題は発生しない。うまく考えてあるのだ。

2016-05-04

私が講義を受ける態度を学んだひと

essay

何年も下書きに眠っていた原稿を公開:

私も高校までは、講義を受ける時、特に積極的な学生ではなかった。むしろ引っ込み思案で、先生が「◯◯が分かる人いますか?」と訊いても、わざわざ手を挙げたりしない、典型的な"日本人型"学生だった。欧米は違う。NHK でも「白熱教室」として放送されたhttp://www.justiceharvard.org/: Sandel 教授の倫理学の講義を見れば分かる通り、先生が問いを投げようならば、クラス中から一斉に手が挙がる。説明に不明瞭が箇所があれば、講義の終わりを待たず、その場で手を挙げて質問をする。TV 番組だから誇張されているのだろうと思う人がいるかもしれないが、そんなことはない。私は今(執筆当時; 2013 年)イギリスの大学院で講義を受けているが、番組の通りである。それでも信用ならぬというなら、Open Courseware の動画を見ればよろしい。

私は大学に入学してから、こういう「欧米型」の積極姿勢を身につけた。そのきっかけをくれたのは、ある同級生だ。仮に A さんと呼ぶことにしよう。A さんは、いつも最前列の中央付近に座っていた。武道系の部活の出身だからか、常に背筋を伸ばしてキリっとしていた。髪の毛はごく短く服装も中性的で、失礼ながらしばらくの間、男性だと思っていたくらいだ(ただし私は人間の性別や年齢を外見から判断するのがおそろしく苦手なので、それを差し引いて考えてほしい)。最初、男だと思っていたおかげで、男子校出身の私も気後れすることなく、隣に座ることができた。

当時の A さんの振る舞いは、分かりやすく言えば、Harry Potter ハリー・ポッター 1 巻の大真面目だった頃の Hermione ハーマイオニーである。当時の私も、「ハーマイオニーだ。ハーマイオニーそのものだ」と思った。先生が現れても教室のざわめきが収まらなければ、「皆さん、静かにしてください!」と注意するのであった。学年の委員を決める時も、自分から学年代表に立候補していた。

A さんで印象に残っていることが 1 つある。それは線形代数学の講義だ。定年間近の厳しい先生で、1コマ目、鐘と同時に入ってくると、ざわめく教室に向かって「大学に入ったからといって弛んでいてはいけない」と一喝した。内容も、1 コマ目こそ対象が医学部であるということも考慮してか、線形性の理論が CT の原理(filtered back projection)に繋がっていくのだという informal な話だったが、2 コマ目からは、係数体を定義してベクトル空間の公理に進むという抽象度の高いものとなった。公理・演繹型の数学にはじめて接する我々は、目を白黒させてしまった。そこへ、A さんは手を挙げて、「すみません、そもそも、係数体とは、何のためのものなのでしょうか」と訊いた。それで、先生もはっとして、導入の動機などを説明しはじめたのだった。

入学して数週間するうちに、私は A さんをすっかり尊敬してしまった。世の中にこんなに立派な人間がいるのかと思った。こんなに真剣に講義を受けている人がいるのかと思った。だが、その頃には A さんが女性であるということを認識したので、声を掛けることができなくなってしまった。男子校出身の私は、女性にどう接していいのやら、皆目分からなかったのである。不用意に女性に声をかければヘンタイだと思われるのではないか、周りから囃されるのではないかと恐れていた。そもそもいつも隣に座るのが失礼に当たるかもしれない。下手をするとストーカーだと思われるかもしれない! そこで一つ離れた席に座ることにした。たまに A さんの方から話しかけてくることがあると、内心大歓喜で、敬語で会話したものだった。私は自分と A さんが同性でないのを悲しんだ。同性だったら、余計な気を使わずに、もっと自由に会話できるのにと。後で聞いた話では、A さんのほうも私がクソマジメに見えて、果たして休み時間に雑談をしていい相手なのか迷っていたそうである。

A さんは学年が進んで専門科目が始まったころから、おしゃれをするようになり、髪も伸ばし始めた。講義には休まず出席していたが、落第して追試になる科目も出てきた。私は、「A さんは堕落した!」と残念に思った。A さんの友人である別の女性は、「君はそういうこと言うけどさ、A さんは今まで大変だったと思うよ。やっと自由になれたというか」と言った。後に聞いた話では、A さんの兄弟姉妹は大変優秀であり、強いプレッシャーを感じていたに違いない。

そうだ、私は A さんを、「リアル・ハーマイオニー」として、その振る舞い・キャラクタ的な意味で勝手なイメージを押し付けて申し訳なかった――と思っている。A さんは、現在、臨床医として活躍されている。

2016-05-04

parallax correction

数学

検出器のセンサー層(CMOS では Si)には有限の厚みがある。これが、高角スポットに半径方向への広がり radial streak を引き起こしたり、高角強度の過大評価につながる。

検出器距離(センサの手前まで)を D, 反射角を $\theta$ , センサの厚さを d とすると、センサ層の奥までの距離は D + d。したがって、光線がセンサを横断する距離は、 $\frac{d}{\cos \theta}$ である。まだ、光線がセンサ表面に達する点は $D \tan \theta$ で、奥に達する点は $(D + d) \tan \theta$ だから、広がりは $d \tan \theta$ である。

仮に、D = 50mm, 反射角 45 度 (センサ半径 50mm に相当), 厚さを 50 um, ピクセルサイズを 50um とすると、スポットの広がりは 50 um だから、1 px 分に相当する。仮に 300 um になると、6px にもなり、なかなか無視できない大きさだ。実際には、光子がセンサー層を横断するうちにエネルギーを失っていくので、奥の方で発生するシグナルは小さくなり、重心位置はここまではズレないけれど。

2016-02-15

Globus Online でファイルを転送する

Linux

スウェーデンにある学術データベースに、数 TB に及ぶデータを登録することになった。過去に同じサーバに sftp で転送した人曰く、1-2 MB/sec しか出なかったうえ、ファイルの破損がいくつかあって大変だったそうだ。幸い、当該データベースが Globus Online の endpoint を設置したとのことだったので、初めての機会だが Globus Online でアップロードしてみた。

Globus Online の背景はややこしくてよく分からないのだが、Grid Computing というクラウド型の分散システム上で動く GridFTP なるプロトコルを、使いやすいように Web ベースのクライアント等と合わせてパッケージ化したサービスらしい。サーバ側は課金制で、クライアント側はフリー。Grid Computing 基盤や Grid FTP 自体は OSS で Ubuntu にもパッケージがあるので、自分で頑張ってインストールすれば Globus Online とは独立したグリッドを構築することもできるようだ。

Globus Online の使い方は Web にあるチュートリアルに書いてあるのだが、scp 等とは仕組みが違うので最初は混乱した。scp では、サーバとクライアントがあって、サーバに向かってセッションを張ってデータを転送する。Globus では末端は end point と呼ばれ、(サービス的にはともかく)原理的にはサーバとクライアントの区別はない。クライアント側(今回、アップロードしようとするデータがある手元のマシン)を end point にするために、globusconnectpersonal というプログラムを導入する必要がある。具体的には、Globus Online の Web からアカウントを作った後、end point 追加を選ぶ。すると、キーが発行されるので、これをクライアントに設定する。さらに、Globus Online から閲覧できるフォルダを -restrict-paths オプションで指定して起動する。

./globusconnectpersonal -start -restrict-paths /path/

このプログラムはデーモンとして常駐する。止めるときは、

./globusconnectpersonal -stop

だ。

次に Web インターフェイス、または ssh から使える CLI から、転送 job (transfer) を submit する。この操作は、手元の end point から行う必要はない。任意のマシンから可能だ(第三者転送という)。

end point に接続する時は、ユーザ名とパスワードを入力して activate する必要がある。activation には有効期限があり、それを経過すると転送が止まってしまうので、忘れずに re-activate する。

いったん job を開始したら、job を submit したマシンは止めてしまっても問題ない*1。例えば、ラボにあるワークステーションを end point にしたら、下宿にあるノート PC から Globus Online の Web にアクセスし、ラボのワークステーション上のファイルを海外のサーバへコピーせよという指示を出すことができる。そうしたら、ノート PC の接続は解除してしまってよい。仮に同じことを rsync でやる場合は、ノート PC からワークステーションに ssh し、rsync を実行する。ノート PC からの ssh が切れても転送が続くようにするには、rsync を screen や nohup の上で走らせる必要がある。それに比べると簡単だ。

job は、エラーから自動で回復する。例えば、activation が expire したとか、end point の片方が切断されたとかで転送が止まってしまっても、reactivation したり、接続が回復すれば、勝手に再開してくれる。これを fire-and-forget というらしい。

Globus Online には CLI インターフェイスがあり、scp など、見慣れたコマンドが使えるように見えるが、これは見かけが同じだけで、中身は違うことに注意されたい。例えば、end point でないマシンにあるローカルファイルを CLI インターフェイスの scp コマンドを使って転送することはできない。あくまで、転送は end point 間だけだ。

不便なことに、Web インターフェイスからは、フォルダの作成ができない。フォルダを作るには CLI インターフェイスを使う。

$ssh username@cli.globusonline.org
Welcome to globusonline.org, username.  Type 'help' for help.
$ mkdir user#server-endpoint/~/incoming/upload-dir

のようになる。rsync や scp では、user@host:/path/to/file だが、Globus では user#endpoint/path/to/file と指定する。

160907追記: @yam_cpp 氏曰く、Web インターフェイスの ☰ のアイコンからできるそうだ。

転送速度のチューニングについては、広域ネットワークにおける GridFTP による効率的なデータ転送 (KEK Report 2012, 高瀬亘) やデータ転送プロトコル GridFTP の並列TCP コネクション数調整機構 GridFTP-APT (2006, 伊藤建志) など、資料がいろいろあるが、デフォルトのままでも十分な速度が得られたので触っていない。

*1:もちろん、データが置いてある end point は起動しておかないといけない