プログラミング

実装して学ぶデータ処理システム その3 - Kabsch transform の意味

今回は XDS における最大の特徴、Kabsch transform の意味について説明する。論文の 2.3 節に対応する。Kabsch transform は、goniostat 座標系から reflection 座標系への変換である。これまでに出てきた検出器系や実験室系との変換と違い、この変換は局所…

実装して学ぶデータ処理システム その4 - 座標変換の実際

前回(本記事公開時には未公開)、Kabsch 変換の意味について説明した。今回は、実装上の問題点について述べる。DIALS のレポジトリにある James の文章 が参考になる。逆空間の一点とKabsch 変換後の反射座標系での一点は、数学的には一対一に対応している。…

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…

Jupyter notebook で cctbx などを使う

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

実装して学ぶデータ処理システム その2 - 反射が起きる条件

スポット検出の話 が間に入ってしまったが、第 1 回に引き続き、Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 の 2.2 節で述べられている反射条件を説明する。…

bash の関数呼び出しの罠

bash で定義した関数を呼び出すときには、() を付けない。プログラムを呼び出すときだって () を付けないのだから、当たり前といえば当たり前なのだが、もし () をつけてしまうとひどいことになる。 doit() { echo "I am doit()" } doit() echo 1 echo 2 doi…

実装して学ぶデータ処理システム その2 - union-find によるスポット検出

作図に手間取っているので、予定を変更して、スポット検出アルゴリズムについて説明する。元ネタは Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 の 2.5 節だ…

実装して学ぶデータ処理システム その1 - XDS の座標系

XDS の主要なアルゴリズムを手を動かして実装することを通じて、X 線回折のデータ処理について理解を深めよう。サンプルデータは Thaumatin / Diamond Light Source I04 user training を用いることにする。まず、XDS の座標系を学ぶ。nXDS は回転軸の概念が…

iPython notebook から cctbx を使う

年を取ったのでなかなか Python が身につかない――などと言っているわけにもいかなくなってきたので、頑張る。iPython と notebook が便利らしいので、cctbx もここから使いたい。iPython 自体は apt-get install ipython ipython-notebook で入るのだが、こ…

h5py で圧縮

h5py で既存の HDF5 ファイルを圧縮した。現状の HDF5 は、既存のデータセットを削除できない(論理的には削除できるが、その領域を再利用できないし、ファイルサイズが縮むこともない)ので、新たにファイルを作って、圧縮フィルタを通しながら書き込んでいく…

Python のスレッドと os.system について

os.system は、サブプロセスの実行が終わるまでブロックする。& をつけるとブロックしない。Python のスレッドは、daemon フラグが立っていないと、それが終了するまで python インタプリタが残る。daemon フラグが立っていると、他のに合わせて自動で死亡す…

bash の wait と、SHELXE の条件検討

bash のシェルスクリプト内で background job を立ち上げた場合、最後に wait を入れないと、job の完了を待たずにスクリプトが終了してしまう。したがって、job も共倒れになってしまう。例えば、TORQUE を使って 30 並列で SHELXE の条件検討を行うスクリ…

SHELXD のログから CC を抽出

Try 4989, CPU 4, CC All/Weak 33.8 / 11.4, CFOM 45.2, best 57.7, PATFOM 3.53から、CC を抽出するには、 grep "All/Weak" test_fa.lst | tr -s " " | cut -f8,10 -d" " | sed 's/,//' とする。awk だけでも書けるけど。R で 簡易プロットするには、 png("…

ジョブシステムの違い

思いつくところから増やしていく。 Sun Grid Engine (SGE) TORQUE 富士通 説明 qstat qstat pjstat -E ジョブの様子 qhost qnodes pjstat --rsc ノードの様子 要確認 qstat -q pjshowrsc キューの様子 qsub -d . job.sh qsub -d . job.sh (デフォルト) カレ…

SSL-VPN にリダイレクトする Bookmarklet

大学に、電子ジャーナルにリモードアクセスするための SSL-VPN ゲートウェイがある。http://journal/web/site/article.pdf にアクセスするには、https://gateway/server/,targetURL=journal/web/site/article.pdf にアクセスすることになっている。これは、J…

Bash の引数について

test1.sh として echo "arg1 $1" echo "arg2 $2" bash test2.sh $2 test2.sh として echo "arg1 $1" echo "arg2 $2" を用意してテスト。 #bash test.sh 1 "2 3" arg1 1 arg2 2 3 arg1 2 arg2 3 となった。本当は、「残りの引数全部」を自然に次のプログラム…

git メモ

本家をクローンして、私家版を GitHub で公開している。新しい開発環境でこの状態を復元した。 git clone https://github.com/biochem-fan/repo git checkout -b scaling origin/branch git remote add official git://official/repo git fetch git branch -…

cctbx で対称操作を列挙する

cctbx で直交座標系での対称操作を列挙する方法を訊かれて書いたコード。 from scitbx import matrix from cctbx import crystal # example is 4BWY symm = crystal.symmetry((197.725, 197.725, 562.576, 90, 90, 120), space_group="H32") uc = symm.unit_…

C で socket

さっき Python で書いたやつ を、C で書きなおした。「さっき」といっても、この記事は1年近く下書きに眠っていた!! #include <stdio.h> #include <sys/socket.h> #include <sys/types.h> #include <netinet/in.h> int main() { int sock, conn; int socket_yes = 1; struct sockaddr_in server_addr, client_ad</netinet/in.h></sys/types.h></sys/socket.h></stdio.h>…

偏微分の変数変換

これも荒いが、下書きが増えすぎているので公開してしまおう。EMAN 氏の http://homepage2.nifty.com/eman/analytic/bibun.html を Maxima でなぞってみる。なお、実用的には以前の記事の方がはるかに簡単。 define(x(r,t,p), r*sin(t)*cos(p)); define(y(r,…

Three.js のメモ

ワールド行列 Object3D.matrixWorld は、それ自身の matrix を含んだ変換である。更新を強制するには、 obj.matrixWorldNeedsUpdate = true; obj.updateMatrixWorld(); とすればよいはずだが、実際に render するまで更新されないことがある気配(詳細未確認)…

STAP細胞関係のゲノムデータを解析してみた

本記事の目的と注意 注意! 私は、NGS については amplicon sequencing の解析経験(しかも半年)しかない。本記事は、データを解析して、STAP論文(Obokata et al, Nature 2014. Article と Letter)に対して何らかの結論を導くのが目的ではない。これだけリード…

隠れマルコフモデル(HMM) のアルゴリズムまとめ

forward/backward algorithm は、HMMのパラメータ(初期確率、遷移行列、出力行列)と出力系列が与えられたときに、各ステップでの各内部状態の確率を求める。Viterbi algorithm は、HMMのパラメータ(初期確率、遷移行列、出力行列)と出力系列が与えられたとき…

cctbx における crystal.symmetry や miller.set の周辺

久しぶりにやや精神の具合がよいので、昔のメモを復習しつつ公開。こんなことをして何になるのかという気持ちは拭えないが、もうどうでもよいのだ。crystal.symmetry は、結晶の格子定数と空間群を管理するクラス。ソースコードは cctbx/crystal/__init__.py…

ラボの分析装置を hack して便利にした話――研究以外でもプログラマがどう研究に貢献できるか

以前所属していたラボで書いたプログラムの話をする。学問的な見所はほとんどないけれど、プログラマとしては愉快なタスクだったし、ラボの仕事の効率化に少なからず貢献できたと思っている。これらのプログラムは、今でもほぼ毎日使われているはずだ。エピ…

CrystFEL と XDS の連携

CrystFEL では、XDS に指数付けをさせることができる。原理としては、CrystFEL が発見したスポットリストを SPOTS.XDS として書きだして、それを JOB = IDXREF とした XDS に読み込ませて、結果を IDXREF.LP からパースしているだけ。1フレームごとに XDS を…

CrystFEL の積分に関わるデータ構造について

CrystFEL の git にある最新の実装では、論文に書いていない機能もいろいろ盛りこまれているようなので、調査したメモ。自分向けに書いているので、わかりにくい部分があればコメントで質問してほしい。GPL なソースコードを引用しているので、このエントリ…

CrystFEL の反射に関わるデータ構造について

CrystFEL における反射についてソースコードを読んだメモ。このエントリはGPLとする。コードはlibcrystfel/src/reflist.c である。反射は _refldata 型で表され、これを包む Relection 型の node が Reflist 型の赤黒木に格納されている。木には、同じ指数を…

CrystFEL のピークサーチについて

CrystFEL におけるピークサーチについてソースコードを読んだメモ。このエントリはGPLとする。実際のワークフローでは、Cheetah などのより高速なソフトウェアで画像に反射が含まれる(hit) かどうかを判定し、hit だけを CrystFEL に食わせることになる。lib…

Canvas を WebGL のテクスチャにする

WebGL では、canvas 要素をテクスチャに使用できる。three.js の fragment shader は、中で gl_FragColor = gl_FragColor * texelColor としているので、ビルボード的なものを作るときは板ポリの色を 0xffffff にしておかないとダメ。どうせテクスチャ貼って…