はじめに

興味深く思ったことを書いていこうと思う。何もかもつまらなく感じて辛くなることがあっても、今、面白いと思ったことは本当なので。その一方で、持病が悪い相にある時は、愚痴っぽいことを書いていることも多い。そんな恥をあえて晒している理由は、精神不…

行列の成分の書き出し

当たり前のことしか書いていないが、下書きに眠っていたので公開しておく。行列の計算で、成分を顕わに書き出したいときがある。まず が鉄則。行列が 3 つ以上並んだときも、隣り合う添字同士がペアになり、あぶれた両端の添字が外に出てくる。 転置があると…

Nature に筆頭で出して、英国でパーマネントの職も得たけど、やりがいがなくなったので辞めます

いわゆる退職エントリ。構造生物学のデータ解析を専門とする staff scientist としてのアイデンティティとその困難について。

"The History of the Decline and Fall of the Roman Empire" 引用集

Wikiquote に引用されている部分を翻訳。原著は public domain であり、この訳も public domain とします。Trajan was ambitious of fame; and as long as mankind shall continue to bestow more liberal applause on their destroyers than on their benef…

ベクトルを転置したものの積について

自明なことだけれど、こういう表記でつまづく人がいるのはもったいないので、ここに明記しておく。, とするとき、 は内積である。一方、 は要素ずつの積を並べた行列である。 この行列の各列は に をかけたものに過ぎないから、お互いに定数倍である。したが…

実装して学ぶデータ処理システム その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 が設定されていな…

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

指数付けの目的は、逆空間の基底ベクトルを決定することである。3 つの基底ベクトル(縦ベクトルである)を横に並べて 3 x 3 行列にしたものを方位行列 A と呼ぶ。成分を書き下すと、 である。 方位行列は 9 成分を持つから、自由度も 9 である。そのうち 6 つ…

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

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

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

何年も下書きに眠っていた原稿を公開:私も高校までは、講義を受ける時、特に積極的な学生ではなかった。むしろ引っ込み思案で、先生が「◯◯が分かる人いますか?」と訊いても、わざわざ手を挙げたりしない、典型的な"日本人型"学生だった。欧米は違う。NHK で…

parallax correction

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

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

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

おんどとりのデータを Ganglia に登録する

おんどとりという温度計測デバイスがある。サーバ監視に使っている Ganglia に測定値を送って、プロットするようにした。なお本機には、測定値をクラウドに送信したり、メールで送信したり、FTP サーバにアップする機能もある。しかし、Ganglia に登録するな…

bbcp で高速転送を行う

数 TB に及ぶ膨大な実験データをインターネット経由で転送する必要が生じた。LAN 内に比べてラウンドトリップタイム (RTT) が長いインターネット経由で scp や rsync を使うと、帯域幅を十分に活かすことができない。そこで bbcp を使ってみた。bbcp は bbcp…

DIALS でたくさんのデータセットを一気に処理する

small wedge で測定した大量のデータセットを DIALS で一気に処理してみた。まず、wedge を検出して、wedge ごとに作業フォルダを作る。 TOPDIR=/path/to/topdir/ find $TOPDIR -name "*_000001.img" -a ! -name "scan*.img" > targets.lst i=0 while read f…

R のプロットに複数の軸をつける

R

上図のように R で複数の軸をつけるには、axis の at と labels 引数を用いる。 df <- read.table("/tmp/test") # CTRUNCATE のログとかから切り出したもの par(mar = c(10, 4, 4, 2)) plot(df$V1, df$V2, type="l", main="Wilson plot", lwd=2, ylab="log(I…

R で二次元配列をなめる処理をするときの速度

R

R で、検出器上の各ピクセルの分解能を計算しようとした。最初に書いた愚直なコードは以下の通り。 lambda <- 12400 / 7000 # A distance <- 51.5 * 1000 # um pixel_size <- 50 # um buf_size <- 3000 # px buf_origin <- c(buf_size, buf_size) / 2 reso <…

NVIDIA の GPU 使用率を Ganglia で監視する

NVIDIA の GPU の使用率を、Ganglia で監視するようにした。GPU の状態は nvidia-smi コマンドで取得できる。dmon モードで起動すると、連続的に1エントリ1行の形式で出力してくれる。 nvidia-smi dmon -s pu -d 5 # gpu pwr temp sm mem enc dec mclk pclk …

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

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

RELION に mrcs 形式の particle をインポートする

公開データセットを使って RELION の練習をしている。高分解能データが触ってみたいので、 EMPIAR ID 10029 にある GroEL のシミュレーションデータを使うことにした。EMAN2 で もっと低分解能の GroEL を題材にしたチュートリアル と同じように処理したとこ…

systemd への登録

iptables の設定を行うスクリプト /etc/init.d/firewall を書いた。これを起動時に自動実行したい。Debian jessie では systemd になっているので、サービスとして登録する必要がある。以下のような /etc/systemd/system/firewall.service ファイルを作る。 …

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 は回転軸の概念が…

XDS の謎

XDS には謎が多い。ソースコードが公開されていない以上、いくら悩んでも答えが出ないことも多い。そういう疑問点を、自分が忘れないために列挙しておく。もし答え(かもしれないもの)をご存知のかたは、ぜひ教えてほしい。 GAIN.cbf と BKGINIT.cbf GAIN.cbf…

R から ADSC 形式の画像を書き出す

R

バイナリファイル書き出しの練習として。 writeADSC <- function(filename, img, lambda_in_A=1.0, orgx=ncol(img)/2.0, orgy=nrow(img)/2.0, pixel_size=0.172, camera_length_in_mm=50.0) { width <- ncol(img) height <- nrow(img) conn <- file(filename…

Rcpp 経由で miniCBF 画像を読み込む

R

Rcpp と CBFlib の練習も兼ねてやってみた。エラー処理は無視!cbflib-wrapper.cpp: #include <Rcpp.h> #include <cbf/cbf.h> using namespace Rcpp; using namespace std; // [[Rcpp::export]] NumericMatrix readCBF(string filename) { FILE *fh = fopen(filename.c_str(), "r</cbf/cbf.h></rcpp.h>…

内積や外積の偏微分

簡単だが、後の参照用にまとめておこう。つまり、積の微分の公式が成り立つ。つまり、こちらも積の微分公式が成り立つ。