en_US English ja 日本語 zh_CN 中文

pythonのmatplotlibを使って3Dの棒グラフを作って書き出す

3Dグラフって見にくく感じるたちなので基本的には好きじゃないのだけど、二つの変数に対しての変化を見たい時に、代替手段がないこともある。

最初エクセルで作ろうと思ったんだけど、意外とエクセルは3Dのグラフを作るのは便利じゃなさそうな感じだった。

じゃあまあpythonか、ってことでmatplotlibで3Dグラフを作ってみた。

最初3Dのプロットを作ろうと思って調べていたら、Axes3D.barがなんとなく使えそうな感じ。

外部リンク: matplotlib公式 mplot3d tutorial

だけどこのファンクションで棒グラフを実際に作ってみると、ぺったりした薄い棒グラフであんまり綺麗じゃない。
設定で変えれないかなと色々調べたんだけど、途中でaxes3Dには最初からbar3dっていうファンクションがあることに気づく。

外部リンク: 公式 Demo of 3D bar charts

というわけで以下はこのbar3dを使って3Dの棒グラフを作る一例。

importはmatplotlibとAxes3Dくらい。

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

まずはいつも通り、最初に適当に描画設定をする。

plt.rcParams['font.family']='Arial'
plt.rcParams['font.size']= 14
plt.rcParams['axes.linewidth']=2.5
plt.rcParams['xtick.major.width']=2.5
plt.rcParams['xtick.labelsize']=14
plt.rcParams['ytick.major.width']=2.5
plt.rcParams['ytick.labelsize']=14
plt.rcParams['figure.figsize']=(8,6)

続いて適当にX軸とY軸を設定。同じYの値に二つのデータを突っ込んで見たくてこんな感じの値を与えてみた。
多分もっといい書き方があるはずだけど、まあpyhtonは動けばいいんだよ派なので気にしない。

X = np.array([0,0.5,1,2,3])
Y1 = np.array([1.12])
Y2 = np.array([2.12])
Y3 = np.array([3.12])
Y4 = np.array([0.88])
Y5 = np.array([1.88])
Y6 = np.array([2.88])

続けてZ軸のスタートポジションを設定。
ちょっとデータを入れたくないポジションがあって、スタートポジション外して消してたんだけど、これももっといい消し方があるんと違うかな。でも動けばいいんだよ派なので……(略

Z0 = np.array([-1020,0,0,0,0])
Z01 = np.array([-1020,-1020,0,0,0])
Z = np.array([0,0,0,0,0])

続いてZ軸の高さの値を設定。

Z1 = np.array([0,19,12,8,7])
Z2 = np.array([16,9,8,6,5])
Z3 = np.array([9,8,6,5,4])
Z4 = np.array([0,0,26,16,13])
Z5 = np.array([29,14,16,12,11])
Z6 = np.array([15,16,12,5,8])

様式美からデータのプロットへ。dx,dzで棒グラフの形を決めて、dzは棒グラフの高さ。
そうそう、特定条件のpythonやら何やら使ってると、アルファ値設定しても透明にならないバグがある。
一応修正できるっぽかったけど、直すの結構面倒くさそうだったのでスキップ。そもそもpythonかmatplotlibのバージョン変えれば直るんじゃないかな。

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.bar3d(X,Y4,Z01,dx=0.19,dy=0.2,dz=Z4,shade=True)
ax.bar3d(X,Y5,Z,dx=0.19,dy=0.2,dz=Z5,shade=True)
ax.bar3d(X,Y6,Z,dx=0.19,dy=0.2,dz=Z6,shade=True)
ax.bar3d(X,Y1,Z0,dx=0.19,dy=0.2,dz=Z1,shade=True)
ax.bar3d(X,Y2,Z,dx=0.19,dy=0.2,dz=Z2,shade=True)
ax.bar3d(X,Y3,Z,dx=0.19,dy=0.2,dz=Z3,shade=True)

続けて軸やら表示の微調整。

ax.set_yticks([1,2,3])
ax.set_xticks([0,1,2,3])
ax.set_xlabel('condition1')
ax.set_ylabel('condition2')
plt.xlim(0,3)
plt.ylim(3.5,0.5)
ax.set_zlim(0,50)

最後に絵をセーブして、見せておしまい。

plt.savefig('FIGURE/'+'orig.png', bbox_inches='tight')
plt.show()

でグラフはこんな感じになる。

3Dバーの例

まあそんなに見やすくもないけど、なんとなく雰囲気は掴める3D棒グラフになりました。

最後はうちのKatja。

Katjaのご尊顔

棚の上に落ち着いたKatja

メリークリスマス!

関連記事

1. pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonのmatplotlibで縦軸二つの散布図を作って書き出す

久々にpython投稿。

ちょっと論文を書くのにy軸二つ並べて、モデルから計算した値と実験値を並べたかったので、pythonで描写。

シンプルに公式ホームページの通りに書いていけば問題はないと思う。

外部リンク: Plots with different scales

importはmatplotlibの設定くらい。

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

まずは適当に好みの図の描写設定。

plt.rcParams['font.family']='Arial'
plt.rcParams['font.size']= 14
plt.rcParams['axes.linewidth']=2.5
plt.rcParams['xtick.major.width']=2.5
plt.rcParams['xtick.labelsize']=14
plt.rcParams['ytick.major.width']=2.5
plt.rcParams['ytick.labelsize']=14
plt.rcParams['figure.figsize']=(3,4)

それから適当にプロットするモデルのデータ。

X1 = np.array([-0.5])
SLD1 = np.array([15.8])
SLD2 = np.array([6.25])
SLD3 = np.array([11.53])
SLD4 = np.array([2.9])

いつもの。

fig, ax1 = plt.subplots()

ガスガス1軸目をプロット。それから目盛りのサイズやらラベルやらの調整。

ax1.plot(X1,SLD1,marker='v',ms=11,color='orange',fillstyle='none')
ax1.plot(X1,SLD2,marker='D',ms=9,color='g',fillstyle='none')
ax1.plot(X1,SLD3,marker='x',ms=9,color='m')
ax1.plot(X1,SLD4,marker='+',ms=9,color='b')
ax1.tick_params(axis ='x', which ='both', top='off',bottom='off', pad=10)
ax1.set_xticklabels([])
ax1.set_ylabel(r"$\Delta$"+'SLD$^{2}$'+'x'+'10$^{20}$'+ ' [cm$^{-2}$]',fontname='Arial')
ax1.set_ylim(0,20)

2軸目追加。

ax2 = ax1.twinx()

2軸目用の実験データ。

X2 = np.array([0.5])
I1 = np.array([388])
I2 = np.array([187])
I3 = np.array([313])
I4 = np.array([35])

2軸目データのプロットと軸のメモリとラベルの設定。

ax2.plot(X2,I1,marker='v',ms=11,color='orange',fillstyle='none')
ax2.plot(X2,I2,marker='D',ms=9,color='g',fillstyle='none')
ax2.plot(X2,I3,marker='x',ms=9,color='m')
ax2.plot(X2,I4,marker='+',ms=9,color='b')
ax2.tick_params(axis ='x', which ='both', top='off',bottom='off', pad=10)
ax2.set_xticklabels([])
plt.xlim(-1,1)
ax2.set_ylabel('INT [cm$^{-1}$]')
ax2.set_ylim(0,400)

最後に絵を出力。figを使っても使わなくても出力は変わんない?

fig.tight_layout()
fig.savefig('FIGURE/test1.png', bbox_inches='tight')
plt.savefig('FIGURE/test2.png', bbox_inches='tight')
plt.show()

絵はこんな感じ。

y2軸のプロット

というわけで、そんなにスケールがあってなかったよ、という絵。

関連記事

1. pythonのまとめ

D

 

en_US English ja 日本語 zh_CN 中文

科学装置購入の事前情報収集のための企業訪問

ここのところ大きな予算の当たっている部署にいるため、装置購入の下見に企業訪問する機会が二度あったので、その時のメモ。

中欧の分析系装置の企業

こちらは中欧にある某分析装置関連の企業を見学しにお邪魔した時の話。買おうとしてたブツはX線散乱装置。まあそれなりの予算がつかないと買えない結構なお買い物。2泊3日の旅程で、2日目がまる1日企業訪問。

企業訪問のアポイントメント

こちらは全部同行のどっかの研究室のボスと企業との間のやり取り。
結構前からコンタクトはあったみたいだけど、予算がそろそろほんとにおりそうだから一度装置を見せてねっていう調子。さほど難しいやり取りはせずに日程だけ決めて、ついでに私ともうひとりの同行研究者が参加できる日程を決めてやり取りはおしまい。
この時点では企業内で1日を過ごすスケジュール、それから夜ご飯はご一緒しましょうねというくらいの決まりごとくらいだった。
それから自分のサンプルをテストできるとのことで、それを持って行くのを忘れないってのが重要なポイント。

企業訪問当日のスケジュール

ホテルから企業へ

朝方会社から15分ほど離れたところのホテルまで、担当の人が車で迎えに来てくれた。
結構ヨーロッパのX線散乱系の学会やワークショップに顔を出している人のようで、私以外の二人は面識があったよう。
小さな車に計4人で乗り込んで、会社までドライブ。

到着したのはなかなか立派な敷地の会社。大企業ってほどではないけれど、分析装置屋さんとしては、結構な規模の会社かな。受付はスルーしてそのまま装置のある実験部屋へ。
荷物を置いたらまずコーヒー。クッキーなどと、エスプレッソマシンとか、ジュースとかお水とかあった。こちらのコーナーは企業のお客さん用ってよりは、企業の人が使ってる場所っぽかったけど。

午前中はもっぱら実験

それはともかくとして、まずは肝心の装置の下見。
購入リストにあるまさにその装置がラボにはあった。まあ新製品の装置を買おうと思ったら大体はあるのでしょうね。
早速持って来たサンプルをみんなでテスト。こっからは結構時間がかかった。やっぱり新しい装置を買おうとしているメンツだけあって、みんなして結構ピッキーなサンプルを持参。

サンプルホルダーへのマウントから、装置の取り扱いは全て企業の人が執り行った。事故防止のためってこともあるのかな。
でも多種多様なサンプルホルダーや、特殊なサンプル環境をいくつか見せてもらった。ホルダーは液体や様々な形状の固体など全く問題ないし、最近の装置は研究室レベルでもオートチェンジャーなど足回りも良い。

一通りサンプルホルダーの取り付け方から、装置・ソフトウェアの動かし方を見せてもらいながら測定開始した。
最初は同行のボスが1時間ほどサンプルを測定。やっぱり新しい装置だけあって測定も早いしデータも綺麗。ソフトウェアも割とサクサク解析が進められそうだった。

続けて私のサンプル。
測定の足回りやデータの精度・範囲は十分良さそうなことが確認できた。だけどディテクターのギャップやサンプルの置き方の都合で、データ測定に多少は制限がありそうだった。まあシンクロトロン行ったってそんなことはあるんだから、どんな装置でも多少はしょうがない。

私のサンプルが終わったあたりでちょうどお昼時。休憩を取ることにした。

企業のレストランでランチ

この企業はヨーロッパらしいというか、素敵な社内レストランが付いていた。サラダのブッフェと、肉料理や魚料理、野菜料理のプレートを頼めるカウンター、それからパンやドリンクなど。ちなみにこちらは経費で落ちるからとのこと。

料理はタダ飯だからってわけではなく、全般にすごく美味しかった。ちょっとこの企業に就職したくなったくらい。ヨーロッパの企業らしくワインやビールもあったけど、この日はさすがにパス。

午後は実験と企業内見学

再び実験室に戻って今度は同僚のサンプル測定。
同僚は強度の変換やサンプル環境に少し思うところがあったようだけど、私たちの共通の見解としては、全体的にレベルが高いし、不利なところも少しの工夫でなんとかなるだろうというものだった。というわけでトライアルはとても好印象。

実験が終わった後は企業内見学。これはまあ全般に会社の宣伝だった。手広く分析装置を扱っているようで、いくつか興味のある製品があったので私は楽しめた。
それから企業内のマシンショップも見せてくれた。この手の企業としては普通なのかもしれないけれど、結構立派なマシンショップ。サンプルの足回りを変えたいとかは、相談してくれたらすぐ対応できるって言ってたのは、ここで作れるってことだろう。

会社の人と街中でディナー

そんな感じでおしまいになった企業訪問。ディナーは街中で会う約束をして一旦ホテルに戻ることに。企業からホテルまではタクシーで送ってくれた。こちらも経費払い。

ディナーは装置の説明をしてくれたエンジニアのオススメのレストラン。食事もそこそこ美味しかったし、ワインも美味しかった。んでもってこちらも経費。
この夕食で企業の人とはお別れになった。その後我々3人は街中のバー巡りに旅立ったわけだけど。

というわけでなかなか満足できる企業訪問だった。初めての企業訪問でどんなことするのかなとちょっとドキドキしてたのだけど、装置の説明や実際の使用・測定にほとんどの時間を費やせたのでとても良かった。
大きいお買い物なので結局入札プロセスになるのだけど、ここの装置をメインに考えて書類を作っても良いなと思えるくらい良い性能の装置だったし、アフターサービスなども問題なさそう。

アメリカのオーブン系の企業

こちらはアメリカでの学会のついでに、巨大オーブン装置界隈で割と有名な企業を訪問。直属のボスとの二人での訪問だった。
こちらは企業訪問が1日。そしてこの企業の装置を実際に使っている研究所の訪問が1日という日程。

企業訪問のアポイントメント

こちらもボスが結構前からコンタクト自体はあったんだけど、予算がついたしアメリカの学会いくしってことで、企業訪問を計画した。

ボスがメールで数回やりとりして日程などを調整したんだけど、当日の予定なんかに関しては直前まで謎だった。漠然としたスケジュールくらい言ってくれたら楽だったんだけど、アメリカの企業ってそんなものだろうか。
それからアメリカの企業訪問の場合は、事前にバックグラウンドチェックなどが入るそうで、名前やそれなり個人情報を送るように要請された。そこらへんカッチリしてるのはアメリカ。

企業訪問当日のスケジュール

ホテルから企業へ

この企業は空港のすぐそば。
というわけで前日に止まっていた空港ホテルから歩いて訪問。アメリカで大都市じゃないのに車を使わないでの移動ってかなりレア。一応歩行者信号ついてたけど、だだっ広い道路を横断するのがちょっと怖かった。

企業にたどり着いたら受付でチェックイン。
我々時間通りにたどり着いたのだけどちょっと待たされた。その後メールしてた担当者と会うことができて会議室に通される。そこでまたちょっと待たされた後に、セールスエンジニアの人が現れる。ちなみに部屋にはコーヒーと水が準備されてた。

くるはずだったもう一人のエンジニアが来れないみたいな話がありながら、向こうの予定通りのスケジュールをこなすことに。
担当のセールスエンジニアの言うところによれば、向こうが10分ほど装置などに関するスライドの発表をして、それからこちらも10分ほど研究内容のスライド発表をして基本的な情報交換。
それから適当に質疑応答しながら意見を酌み交わした。

まあ向こうがどういうことができるとか、こちらがどういう意図を持っているっていう情報を交換はできたんだけど、2vs2しかいない少人数でこういう形でやるべきだったのかっていうのはちょっと疑問だったかな。
まあ購入した場合にどういうことができるとか、どんなアップグレードができるとか、大まか必要な予算規模とかは情報交換できたのでよかったけど。

そのあとは向こうの企業の社内見学。こちらの会社の場合は最新機器を置いているってよりは、古いオーブンや、フレクシブルなセットアップが可能なカスタムオーブンみたいなのを会社の中に置いていた。

パネラのデリバリーでランチ

久しぶりのパネラブレッド。パネラブレッド好き。美味しい。企業の経費で。

企業の他のグループの人たちも混ざりながら、雑談をしつつの昼食。ちなみに朝の会議にいるはずだったエンジニアも来てた。まあアメリカの企業へのヨーロッパからのお客さんのお迎えなんてそんなものなのかもしれない。

んでお昼が終わったらディスカッションもおしまい。あっさりしているようだけど、自分のサンプル持ってきてテストができるとかってんでなければ、まあそんなに話すこともないよなってのも確か。

翌日に研究所に同行するので、この日はそのまま空港へと向かったのだった。

アメリカとヨーロッパの企業を訪問しての雑感

といわけで未経験だった装置購入の準備のための企業訪問を立て続けに2個こなした。
どちらも予算的には結構な規模のお買い物なので、装置の選択に間違いは犯したくない。もちろん最新の装置を買う場合は、そんなに問題があることはないだろうけど、
だけど実際に企業を訪問して細かく話を聞いてみると、やっぱり新たな情報や気づきがあるし、各社製品にそれなりの違いはあった。
特におっきなお買い物の場合は、多少旅費をかけてでも事前にしっかり装置や会社を見ておいた方がいいだろう。

あとは雑感としては、ヨーロッパの企業の方がお客さん感があったかな。アメリカの企業はビジネスパートナーとしてこちらを見ているって感じだった。それがお国柄の違いなのか、分析装置とオーブンっていう業種の違いによるものなのかはわからないけれど。

関連記事

1. python・科学記事のまとめ

D

en_US English ja 日本語 zh_CN 中文

フィンランドで装置付きの研究者職に応募してみた時のお話

フィンランドの某大学の装置付きの研究者に応募してみた時の経験のメモ。ポジション的には教授職、講師、ポスドク、学生のサポートや、外部の会社とかの外注サンプルの測定とか、そういうことをするポジション。

申請書類の準備と応募

準備する必要があった申請書類は以下のとおり。

  • モチベーションレター
    A4で2枚までにまとめてってこと。志望動機やら自分のバックグラウンドやらを、それなりに一般的なフォーマットを参考に作成。なんでこの大学のこのポジションをやりたいかってことをなるべく書くように注意。
  • CV
    論文リスト付きで4ページ。枚数指定がない場合は厚い方がいいという方針で、結構些細なファンディングとかアブストラクトとかも載せた。
  • 学位証明書
    最終学歴のPhDだけ。学位証明全部とか言われてなければ、マスターとか学部のはいらないはず。
  • ティーチングステートメント
    あんまり教育履歴がないので、経験と興味でA4で2枚に収まってしまった。今度集中講義とかやらせてもらうし、ここをこれから厚くしていきたいところ。ティーチングとファンディングがもう少し厚くなったら、テニュアトラックも面接にお呼ばれするようになるかなあ。
  • 推薦書かプロフェッショナルレフェレンス
    ビッグボス二人から推薦書を書いてもらった。ショートリストに入ったらお二人の影響のおかげかな? せっかくなのでCVにもレフェレンスリストはつけておいた。

リサーチステートメントがないっていうのが不思議なくらいで割と一般的な申請書類だと思う。

申請書類を出してから2週間ほどでショートリスト入りっていう連絡がくる。申請は35人ほどいたそうだ。

ショートリスト入りしてからの準備

意外なようなそうでもないようなって感じでショートリスト入りした。まあ論文の数はそこそこあるし、推薦状も効果を示してくれたのだろう。

面接の前にやっておいたことはスライド作りとX線装置・解析関係のちょっとした復習。だけど理論的なことは全然聞かれなかったから、後者は役には立たなかった。いやいくつか大学の装置を使ってやりたい実験も思いついたし、完全な無駄ではなかったけれど。

インタビューで使うスライド作り

これは昔のボスからのアドバイスで、面接にスライド用意してったらってことだった。せっかく頂いたアドバイスなので、スライド作りに励んだ。ちなみに今回の面接ではスライドを使うことはできたけど、担当の人にパソコンやプロジェクターの有無を確認しておくといいかもしれない。

job descriptionに乗ってるテクニックの確認

よくショートリスト入りしたなってくらいジョブディスクリプションに使ってないテクニックもあったので、一応それらをお勉強。だけど面接はそういう知識や経験を聞くような質問は全くなかった。

インタビュー本番に臨んで

ショートリスト入りの翌週にはインタビュー。とても早いインタビュープロセスだった。

当日会場の建物に入ったらすごく分かりにくい部屋がインタビューの会場。階段の脇にくっついてるような小さな部屋で、本当にすごく分かりにくかった。ちょっと余裕を持って会場入りしておいてよかった。

なんとかたどり着いた時には10分前。すぐに担当のHRの人がきてくれて、中に入れてくれた。
中にいたのは申請先の学部から教授、スタッフサイエンティスト2名、HR、それから他学部の研究者が一人。そんな感じのコミッティー。

まず最初に自分で準備しておいたプレゼンを発表した。これはスライドで自己紹介していいって自分から聞いたからで、向こうから求められていたわけではない。コンピュータがプロジェクターにつながらないトラブルはあったけれど、会場のコンピュータを借りて問題なくスタート。
プレゼンは問題もなくスラスラ発表できたけど、ちょっとアウトオブフォーカスだったっぽい。スライドは自分の長所で申請先に対してどんなプラスアルファができるかっていう感じで作ってた。だけど彼らの求めてるのはプラスαってよりは、彼らの手足となって職務を忠実にこなすタイプだった模様。
私としては正直そっちの方がやりやすい仕事なんだけど、欧米は基本的にいかにプラスαを出せるかっていうところをアピールする方が効くのかと思ってたんだけどね。
まあそこまでマイナスにはなってないだろうけど、特にやらなくてもよかったかなっていうプレゼンだった。

インタビューの質疑応答で聞かれたこと

1 このポジション自分の研究全くできなるけどいい? 教授系目指してるんならやめといたがいいんじゃない?

HRからの質問。これが一番本気の質問だった。絶対自分の研究はさせないぞっていう強い意志を感じた。ちょっとくらい大丈夫でしょ、みたいなことを言ったら、そんな時間ないない……みたいなことをみんな口を揃えて言ってたし。
まあ前に同じポジションにいた研究者がとても忙しそうだったのを知ってるので、入って見たら自分の研究やってる時間ないってこともあるんだろうなとは思っていた。だけど面接の時から全く時間取らせないぜって言ってくるのは想定外。
続けてファーストの論文が増えないから、テニュアトラックの道も潰れる可能性が高い。それに納得できないんだったら考え直した方がいいぜっていうだめ押し。
まあ装置に本気で集中して手足になってくれる研究者を取りたいんだろうなっていう強い気持ちを感じた。といわけで多少(せめて3割くらい?)は自分の研究したいなって思ってた私には結構な先制パンチ。ちょっと答えるのにどもった上に、その後しばらくの間ダメージが残ってしまった。今思えば申請書類にリサーチステートメントがなかったっていうのは、きっとこのへんが理由だったんだろうな。

確かに採用する側としてはどういう目的で働くかかってのは重要なんだろうな。能力自体はある程度はCVとかカバーレターとかから想像はできるしね。でもこの質問が最後だったらもうちょっとインタビュースムーズにできただろうなとは思った。叶わない願望。

2 CVとかモチベーションレターに書いたバックグラウンドの詳細

HRからの質問。装置研究系のお仕事だったので、装置のトラブルシューティングとか前の職場でしてたユーザートレーニングの詳細を詳しく聞かれた。
一方で解析の理論とかそういう詳細についてはノータッチ。個人的には外注のサンプルをもらったら実験よりも解析で詰まる人の方が多いんじゃないかなって気がしてるんだけど、向こうとしては実験重視ってことだったのかな。

この質問のあたりもうちょっと色々答えてたと思うんだけど、前の質問のダメージのせいかあんまり覚えてない。だいたいつつがなく答えてたと思うので、この質問についてはややプラスかプラスマイナス0ってとこか。

3 ティーチングの経験、フィロソフィーそれから訓練

HRからの質問。経験についてはそもそもそんなにないし、頭の中にまとまっていたことを話せたんだけど、ティーチングフィロソフィーを言葉で聞かれてちょっと詰まった。
ティーチングステートメントにそれなりに書いてはいたと思うんだけど、ちょっと思い出せずにすごく適当にぶらぶらぶらと話した。多分結構ひどかったと思う。。
これは次の就活の機会での反省点。聞かれる可能性はだいぶある質問ではあったし、もうちょっと準備しておいてしかるべきだったと。

それから大学とかでやってる生徒の教え方の講義とってる? って聞かれたので、大学とかでそういう講義があれば、とっておくと就活で役に立つかもしれない。

この教育関連の質問はややマイナスだったかな。

3 入ったら1年目・2年目のプランはどんな感じ?

スタッフサイエンティスト1からの質問。想定はしていなかった質問だけど、想定しててもよかったかなっていう質問。
まあユーザーサポートのお仕事なので、装置に慣れながら外注の仕事を受けますよっていう当たり障りのない回答。2年目のプランを答えるのがちょっと難しかった。時間があればプログラミングもしますよって言っておいたけど、ちょっと向こうの求めてる方向性とは違かったかなって気もする。

ややプラス?

4 何かお仕事を通して学びたいことはあるの?

これは僕がなんかのタイミングで、新しいこと勉強できたらいいなって言ったのに反応したスタッフサイエンティスト1からの質問。

特殊なサンプル系とかテクニックとかやったことのないの試したい的なことを言っておいた。

プラスマイナス0。

5 チームワークと個人主義どっちのが好み?

この割とどこでも見る質問。昔はチームワーク重視してるって答えが欲しいのかなって思ってたんだけど、最近ヨーロッパだとそうでもないのかなって気がしている。
個人技もチームワークも状況によって使い分けれまっせ、ていう答えを用意しておくといいのかも。

プラスマイナス0。

6 何か質問ある?

一番聞きたかった研究・ユーザー補助・教育・その他のバランスが、向こうからの質問で解決していたので、とりわけ聞きたいことはなかった。

なので、外注のサンプルがどんな感じのがあるか、申請先のグループがとりわけ求めている実験手法や解析方法はあるのか、採用されたら開始日はいつだとか、そんな当たり障りのない質問だけしておいた。

インタビューを終えての印象

終わったら60分は超えていた面接だったので、他にももう少し質問もあったと思うんだけど、ちょっとこれ以上思い出せず。
緊張してると頭が回りすぎて記憶が飛ぶ。

終わってみたらそんなに当たり障りのない感じの面接だったとは思う。ティーチングの質問にパッと答えなかったところで若干マイナスって感じ。書類でぶっちぎってたら可能性はあるかもしれないけど、ちょっと上位互換レベルの研究者が出してるのを知ってるので、彼と比べるとちょっと分が悪いだろう。

ちなみに4人ショートリストに残っていて、もしかしたらもう少し増やすかもしれないとか、不思議なことを言っていた。
まあリストを増やすってことは今のリストの人はあんまり取りたくないのかもしれないな。私やその知り合いも含めてだけど、少しは自分の研究もしたいですって感じでそれなりに業績リストの良い人が最初に集まっちゃったのかなという印象。

関連記事

1. 海外就職手続きのまとめ

2. python・科学記事のまとめ

D

en_US English ja 日本語 zh_CN 中文

サイエンスセミナーで座長をしてみよう

恥ずかしながらこの歳にして初めてサイエンスセミナーの座長を経験することになった。この投稿ではその時の経験についてまとめていきたい。

私はもともと自分の発表ですら大の苦手、というか人前に立つだけでもかなり心苦しく感じるレベルのメンタリティー。ついでに複数の物事を同時にやりくりして、上手に進行させるようなマルティタスクができない方。正直座長というのは能力的に無理かなと思っていた。

ところが数ヶ月前のある日、うちのボスから学生の年末の発表会の座長をしないかと誘われる。この手のスキル皆無だけどやってもいいの?って聞いたら、身内の学生の発表会だし、練習には一番向いてるでしょて感じのご返答。
確かにphdの学生の年末発表会のようなものなので、他でやるよりはきっと良いのでしょう、ということで好意に甘えて座長を体験させてもらうことにした。

結果としては想定通りなかなかにグダグダな経験になってしまったのだけど、そういう経験こそシェアする価値があるかなってことで、筆をとることにした。あまりこの手のが得意じゃない人の参考になれば良いかなと思いつつ。

サイエンスセミナーで座長をするのに準備したこと

というわけでまずは下調べ。
受け持つのはPhDの学生の発表というわけで、まず彼らの紹介に使うmini-bioの作成。ここら辺はオーガナイザーが紹介用のカテゴリーを作ってくれた。必要な情報は以下の通り。

  • 学生の名前
  • 修士論文のタイトル/もしくは内容
  • いつこのラボに移動して来たか
  • 博士論文のタイトル
  • ギャグ1個

発表前には気づかなかったのだけど、これが実は時間的に結構重たかった。修士論文のタイトルやPhDのテーマとかって大概にして結構長いので、これを説明するだけで45秒前後が過ぎてしまうということ。

それから英語の場合重要なのが学生の名前の発音。フィンランドの大学も大概多国籍なので発音がとても多様で難しい。
直接本人に発音が聞ければベストなんだけど、受け持つ学生が決まったのは発表の数日前。全員に聞いて回っている時間はなかったので、ちょっと発音に自信のない学生もいた。
しょうがないのでそういう名前についてはネットで検索してみた。どうやら最近はオンラインで結構マイナーな名前まで発音が調べられるようになっているよう。

それからギャグは学生も嫌だったようで、ほとんどの学生は送ってこなかったのでスキップ。
ちなみに担当の学生は6人と多いので暗記は諦めて各項目について、読む用のメモを最初から作ってしまうことにした。

迎えたセミナー当日のドタバタした話

セミナー前の準備

さて迎えた当日。

まずは朝始まる10分ほど前にセミナーのオーガナイザーとお話。
基本的にはトーク8分・質問2分の10分構成だったので、本当はこの時間の中に学生の紹介も詰め込まないといけないのだろうけど、今回はうちうちのセミナーなのでmini-bioを読み上げるところの時間は、時間から抜いちゃっていいとのこと。だからこの時点で30秒x6で3分遅れるのは許容範囲ということ。

サイエンスセミナーの開始

予定が狂い出したのはまずはビッグボスのトークから始まるオープニング。
彼の持ち時間は15分と結構長かったので、正直ここで2・3分くらいの時間の余裕ができるんじゃないかと思っていた。のだけれども実際は彼がしっかり喋りに喋った。
余裕ができるどろこか始まりですでに時間が2分オーバー。

さて一応私も名前とプロジェクトだけの簡単な自己紹介。それからセッションを一言で紹介して、経過時間お知らせのタイミングと方法だけ説明。
これが終わった時点ですでに3分くらいオーバーしていた。

というわけでこっからメインパートの始まり。あとはひたすら以下の繰り返し。

  1. 学生の紹介
    これは準備しておいた通りに紹介をするだけ。30ー45秒ほど。特に問題なし。
  2. 学生の発表
    ここでストップウォッチを押す。私の担当の学生は結構自由だった。7分で終わった人から9分30秒喋った人まで。分布広かった。次のセッションの学生たちは、だいたい8分前後で終わっていたので、ここら辺は結構運次第。学会でも同じようなものかしらね。
  3. 5分と7分で学生へ残り時間の案内
    今回は用意されていたサインを発表中の学生に見せる形式。大画面を見ながら喋るタイプの学生はなかなかサインに気づいてくれないもの。後で思ったのは、そういう学生には普通に声を出してしまえばよかったのだった。
  4. 学生の発表が終わったら拍手
    拍手ってあるタイプのセミナーとないタイプのセミナーがある?今回はオーガナイザーに確認したら、発表後・質問後・セッション後の拍手を入れて良いとのこと。
  5. 質疑応答
    拍手が終わったら質疑の開始。簡単なようで一番やっかいなところ。残り時間が2分近くある人は「Do you have questions?」「any questions?」、1分くらいで「A question?」、30秒くらいから「one quick question?」くらいの感じでやってみた。
    問題となったのは当たり前といえば当たり前なんだけど、みんな座長のことなんて気にしていないということ。クイッククエスチョンって言ってるのに、平気で長い質問が出てくるし。ちょうど時間が切れたから「Thank you」と終わりにしようとしたら、ビッグボスが急に長いコメントをつけだしたり。研究者らしくてみんなとても自由。
  6. もう1回拍手をして次の学生へ
    ここは特に問題無し。

全員終わったら発表者全員と会場に来てくれた人に感謝を述べてもう1回拍手をしておしまい。

少しずつ時間が押していったせいで、最後はちょうど10分オーバーくらい。研究室のビッグボスが質問最後に喋り出した1分半と、クイッククエスチョンて言ってるのに最後に2分くらい喋ってたどっかの先生がほとんど全てだなんけど。
まあ最初の始まる前の3分と紹介での3分が許可されていたと考えれば、実質推してたのは4分ほど。

でも最後にオーガナイザーと少し話したら、遠慮せずに質問はぱっぱと切ってもよかったと言っていたので、やはり10分はちょっと押しすぎたかな。
伸びてた時間を調整するように、臨機応変に対応できればよかったのだろうけれど、新米座長にはちょっと荷が重かった。
もう少し経験を積めばそういうこともできるのかしらね。とはいえもとから伸びてる時間で発表の時間を奪ってしまうのもなんだかあれだし、やっぱりちょっと難しいな。
今回に限って言えばスタートが遅れていたのだから、紹介の時間だけでも発表時間に含めてしまった方がよかったとは思ったけど。

サイエンスセミナーを終えて

というわけであまりうまくさばくことのできなかった初座長。たぶん次の機会があればもう少しうまくできるだろうとは思うけれど、次の機会はあるのだろうか。
やってみた感想としては、正直結構面白い体験だった。時間の配分以外はだいたいよくできたと思うので、次の時はここの計画をもうちょっとしっかりと考えて行きたいと思う。
ただこんなおっさんになってから初めてやるよりは、ある程度若いうちに経験しておいが方が上手になるタイプのものだとも思う。ポスドク始めたくらいでチャンスがあればやっておくと良いかもしれない。

次のセッションの座長を見てて思ったんだけど、残り30秒だろうと1分だろうと、1個か2個質問が出たらもういいかって感じで質問を打ち切っていた。私はなるべくギリギリまでその学生に質問の時間を与えたいと思っていたのだけど、残り30秒くらいのところはもう終わりにしてもよかったかなってのが一番の反省点。

関連記事

1. python・科学記事のまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使ってjpeg画像をまとめて圧縮する方法

久々にwordpressの投稿関連の記事。とは言ってもpython記事だけれども。

さてwordpressユーザーの方はご存知の通り、wordpressのブログにアップできる画像の容量は1Mバイトまで。というわけでデジカメの画像はもちろん、携帯のカメラで撮影した写真でも多少サイズを小さくする必要があることが多い。

数枚ならどんな方法で圧縮してもいいのだけど、これが大量の画像になってくるとそれなりに面倒臭くなっていくもの。
いくつか良いフリーソフトウェアもあるようだけど、こんな時こそpythonを使って見るのも悪く無い。書き出しのファイルサイズを多少はコントロールできたり利点もある。

そういうわけで、この投稿ではpythonを使ってjpeg画像をまとめて1M以下にサイズ圧縮する方法について書いていきたい。

Pillowを使ってjpeg画像圧縮

画像の読み込みはpillowで。画像の読み込みができるpythonパッケージはいくつかあると思うのだけど、簡単な処理をするだけなら大して変わらないんじゃ無いかと信じている。
そんなことないよ。これのが全然早いよってのがあったら教えてください。

必要パッケージのインポート

まずは必要なパッケージをインポート。ファイル名の読み込みようにtkinter各種。それから画像処理にPIL、ファイル名の処理にosをインポート。
ioは画像サイズの処理の関係で利用。詳細は後ほど。

import os
from tkinter import Tk, ttk
import tkinter.filedialog as tkfd
import tkinter.messagebox as tkms
from PIL import Image
from io import BytesIO

tkinterで読み込みファイル名の取得

取り立てて変わったこともないけれど、一応jpgでもjepgでも読み込めるようにしてはある。画像ファイルは./ORIGINALっていうディレクトリに入れておいて一気に読み込めるような設定。

tk = Tk()
tk.withdraw()
print ('select a data files')
filenames = tkfd.askopenfilenames(filetypes= [('jpg','*.jpg'),('jpeg','*.jpeg')], initialdir='./ORIGINAL')
tkms.showinfo('file paths are',filenames)
tk.destroy()

画像サイズを1M以下になるまでjpeg圧縮

というわけでメインパート。画像を縮小するのではなく圧縮する場合だと、jpeg保存のqualityを変えていけばサイズが変わる。この場合qualityを5ずつ落としていっているけれど、そこらへんは適当に調整していただきたい。
for文でqualityを下げながら一回一回ファイルを書き出し、osのstat.st_sizeでファイルサイズのチェックをしている。
ちなみに書き出しは./REDUCEDというディレクトリの中に。

for i in range(len(filenames)):
  img = Image.open(filenames[i])
  basename = os.path.basename(filenames[i])
  for j in range(10):
    img.save('./REDUCED/'+basename+'.jpeg',quality=90-j*5,optimize=True)
    if os.stat('./REDUCED/'+basename+'.jpeg').st_size < 1000000:
    break
  img.close()

しかしこのやり方だと毎回ファイルを書き出して、読み込んでっていうのが少し気持ち悪い。

というわけで以下は実際にはファイルを書き出さない例。ioのByteIOを利用することでそんなことが可能になる。初めて使ったのでこのパッケージの詳細についてはよく把握していない。

qual = [0]*len(filenames)
  for i in range(len(filenames)):
  img = Image.open(filenames[i])
  for j in range(10):
    img_file = BytesIO()
    img.save(img_file,'jpeg',quality=90-j*5,optimize=True)
    #print (img_file.tell())
    if img_file.tell() < 1000000:
      qual[i] = 90-j*5
      break
  basename = os.path.basename(filenames[i])
  img.save('./REDUCED/'+basename+".jpeg",quality=qual[i],optimize=True)
  img.close()

実際にファイルを書きださないで処理しようと思った理由は速度があげられるかなっていうところもあった。

というわけで二つのスクリプトを1.3M、3.5M、4.5Mのjpeg画像を使って速度チェック。
ちなみに書き出されたファイルのサイズは674、989、962[KB]。元のサイズ差のせいなのか、qualityを5ずつ落としているせいなのか目が荒い。

ファイルを書き出す最初のプログラム:

4.87; 4.89; 4.91; 4.87; 4.84

ファイルをメモリに保存しとく2個目のプログラム:

5.34; 5.35; 5.61; 5.22; 5.23

というわけで意外なことに、一々ファイルを書き出すプログラムの方が早かった。メモリを食う方が速度が落ちるのかな。使っているパソコンのシステム次第で変わりそうな気もするけどどうなんだろう。
大した速度の違いがあるわけではないので、個人的には2個目の方が良いかなって思ってるんだけど、どんなもんでしょうか?

変換前後

最後に圧縮前後の画像の一例を載せようと思ったら、当たり前だけど圧縮前の画像がサイズの都合でアップできなかった。なんのためにこの投稿を書いていたのやら。画像を並べたスクリーンショットを参考までに。

関連記事

1. pythonのまとめ

2. pythonでデジカメのraw画像を一括読み込みして別形式で保存する

D

en_US English ja 日本語 zh_CN 中文

pythonでデジカメのraw画像を一括読み込みして別形式で保存する

以前Raw画像を取り込んで処理するっていう投稿を書いた。

内部リンク: pythonでraw形式の画像を読み込んで処理する方法

この投稿はX線データなどの画像処理を書いたものだった。しかし後々考えるとRaw画像の取り込みっていうタイトルだと、写真画像が取り込みたくて訪れてくれてる人もいるのかなとちょっと気になっていた。

それが理由というわけでもないのだけど、この投稿では写真の画像を取り込んで処理するプログラムについて書いていきたい。
こんなプログラムを書いた理由は、容量がきつくなってきた古いパソコンのお掃除がしたかったこと。嫁が大量の写真データを保存していたので、Raw形式で保存しておく必要のないものを一括で読み込んで、jpegかなんかで小さめのファイルで書き出そうという試み。

使用した画像はNikonのデジタルカメラのRaw画像。拡張子はNEF。

rawpyを使ってRaw画像の読み込みと保存

というわけで画像の読み込みに使ったのはLibRawという様々なデジカメ画像を読み込んでくれるライブラリ。pythonの場合だとrawpyという素敵なwrapperがあるのでとても便利。installはpipで。

外部リンク1: LibRawのウェブページ

外部リンク2: rawpyのウェブページ

pip install rawpy

必要なライブラリなどのインポート

まずは必要なパッケージをインポートする。この投稿で使っているのはtkinter, rawpy, imageio, PILのImage,それとos。

from PIL import Image
from tkinter import Tk, ttk
import tkinter.filedialog as tkfd
import tkinter.messagebox as tkms
import rawpy
import imageio
import os

tkinterを使ってファイルの一括読み込み

まずはいつも通りtkinterを使ってファイルを一括読み込みする。

def main():
  tk = Tk()
  tk.withdraw()
  print ('select data files')
  filenames = tkfd.askopenfilenames(filetypes= [('raw','*.NEF')], initialdir='./')
  tkms.showinfo('file paths are',filenames)
  tk.destroy()

これでファイル名の読み込みは完了。読み込みの詳細は別投稿で。

内部リンク: pythonでファイルを読み込む時にファイル名を取得する方法のまとめ

rawpyを使っての読み込みとフォーマットを変えて保存

ファイル名の取得ができたらメインパート。

#forを使って読み込んだファイルを1個1個処理。
for i in range(len(filenames)):
  rawData = rawpy.imread(filenames[i])
  basename = os.path.basename(filenames[i])
#postprocessを使うとnumpyのarrayにデータを詰め込める。最初のオプションはホワイトバランスをカメラの設定で。2個目は明るさを勝手に調整しないっていう設定をオフにしている。
  rgb = rawData.postprocess(use_camera_wb=True,no_auto_bright=False)
#imageioを使ってjpegに書き出すパターン。
  imageio.imsave('test.jpeg',rgb)
#同じくPILを使って。
  img = Image.fromarray(rgb)
  img.save(basename+".jpeg",quality=85,optimize=True)
#保存フォーマットは適当に目的に合わせて。
  img.save(basename+".tiff")

  rawData.close()

というわけでこんな感じでrawpyを使うと、とてもあっさりと写真画像の読み込みができるということがわかった。
上のrgbはnumpyのarrayに詰まっているので、色々と処理をしたかったらここですることができる。
書き出しはまあ何を使ってもすぐ簡単にできるので問題はないでしょう。

rawpyの難点としては少し読み込みが遅いこと。
初期設定を使うと1ファイルの読み込みで4秒ほど、2ファイルをforで回して7秒ほどだった。特にforで遅くなるわけでもなさそう。
demosaic_algorithmを変えたらスピードが結構変わった。どこかウェブの片隅で早いらしいと読んだLINEARを使ったら約2倍の速度になった。何か速さの代償があるのかは知らない。

それでも1ファイル2秒ほどかかるので結構遅い。しかし一回動かしたら放っておけば良いタイプのプログラムなので、別に多少遅くても問題はないでしょう。

左:macのプレビューで書き出し 右:rawpy経由で書き出し

もう一つ問題点としては、ホワイトバランスなど画像読み込みの調整。
これは別にrawpyの問題というわけではなく、Rawデータを開くとどんなソフトだろうと起こる問題だそうな。
Rawデータから元々のデジカメがjpegなどに変換している詳細は企業秘密。カメラのホワイトバランスの値は使うことができるけど、同じ画像にしたかったら詳細をかなり頑張って調整しないといけないのだとか。
以下はそんな感じで色々と設定をいじりながら試していた時の一つ。一番よく映ったオプションというわけではない。
読み込みのオプションの詳細についてはrawpyのホームページに乗っている。

外部リンク: rawpyのオプション

rgb = rawData.postprocess(output_color=rawpy.ColorSpace.sRGB, demosaic_algorithm=rawpy.DemosaicAlgorithm.LINEAR,output_bps=8, use_camera_wb=True,no_auto_bright=True,bright = 1.8,gamma=(1.8,4.0), no_auto_scale=False)

カメラのホワイトバランスを使うuse_camera_wbか、rawpyにホワイトバランスを調整してもらうuse_auto_wbが一番劇的に結果が変わる。
基本的にはこの二つのどちらか好きな方を使ってあとは微調整をする感じだろうか。
それっぽく画像が変換できるオプションを知っている方がいたら教えていただきたいところ。

というわけでこの投稿はひとまずここでおしまい。時間ができたらもう少し他の読み込み方法など試してみたいところ。

関連記事

1. pythonのまとめ

2. pythonを使ってjpeg画像をまとめて圧縮する方法

D

en_US English ja 日本語 zh_CN 中文

pythonでファイルを読み込む時にファイル名を取得する方法のまとめ

久々のpython投稿。

さてだいぶ前にpythonでのテキストファイルの入出力について書いた。

内部リンク: pythonでのテキストデータの入出力方法いろいろ

上記の投稿はどちらかというと開いたファイルの中のテキストデータをどう処理するかということに焦点を絞って書いた。

今回の投稿はちょっとそれとかぶってもしまうところもあるのだけど、最初にファイルを開くところ、ファイル名(もしくはファイルパス)を取得するということに焦点を絞って書いてみたい。

tkinterでファイル名取得

最近私がもっぱらこれを使っているので最初はtkinterから。
自分で使うのにはいちいち窓が立ち上がるし、クリックしなきゃいけないこともあるので鬱陶しさもある。しかし他の人にプログラムを渡すときには、どうやってファイルを読み込むかをあまり説明しなくていいので楽。

また大量のファイルを一度に読み込みたい場合だと、窓でファイルを選択するときにshiftやcontrol/commandを使って大量選択が簡単にできるのも利点。

というわけで以下はtkinterを使ったファイルパス取得の例。

from tkinter import Tk, ttk
from tkinter import filedialog as tkfd

tk = Tk()
tk.withdraw()
print ('select data files')
filenames = tkfd.askopenfilenames(filetypes= [('text','*.xy'),('EXCEL','*.xls')], initialdir='./')
tkms.showinfo('file paths are',filenames)

print ('select a background file')
bkgname = tkfd.askopenfilename(filetypes=[('text','*.xy')], title = 'bkg file', initialdir='./')
tkms.showinfo('bkg file path is', bkgname)

tk.destroy()

大量のデータファイルを読み込むときはaskopenfilenamesを使って、一つのファイルの場合はaskopenfileを使う。
initialdirで窓を開いた時の最初の位置を指定できるので、適当に./DATAみたいなフォルダを作っておいてデータをそこに入れておけば、command+Aで簡単に全部データが開けるかな。

ただし本当に大量のデータを読み込むって場合は、別の方法を使った方が楽かもしれないけど。

tkinterのその他詳しい使い方については別投稿を参照していただきたい。

内部リンク: pythonのtkinterを使ってユーザー入力を取得する方法

argvsでファイル名取得

C言語の時は一番お世話になっていたargvs。

自分で2、3個のデータファイルを読み込もうとする時には一番楽な方法。
Macのターミナルでtabでの名前補完と組み合わせれば、ファイルの読み込みがサクサクとできる。
ただ大量に読み込ませたいファイルがある場合は大変。せいぜい5個くらいまでの読み込みに使う感じだろうか。

import sys
argvs =sys.argv
argc = len(argvs)
if (argc != 3):
  print ("input filename and output name")
  quit()
file = open(argvs[1], 'rb')

pythonではsysをimportしておく必要がある。

if文は読み込みたいインプットの数が正しい数であるかを確認。間違ってたら正しいインプットを要求して落とす。

一番目に読み込んだ名前はargvs[1]になるので注意が必要。argvs[0]はpythonの.pyファイルになる。argvsは他にも単純に適当な文字列を追加するのにも使える。上の例を書いた時には、argvs[2]に出力ファイルネームを入れていた。

そうしてファイル名を取得したらopenやnp.loadtxtなどで開いてあげると良い。

インプット用のファイルを作っておいてファイル名はプログラム内に

プログラムに直接ファイル名を書き込んでしまっておく。必ず何かしらパラメータの入力を受け付けるタイプのプログラムでは便利。input.txtなどのファイルを作っておいて、プログラムを走らせる前にその中のパラメータだけ書き換える感じ。

プログラムの方を書き換えないで良いのが利点と言える。

ファイル名はパスでも指定できるし、同じ階層にいるなら名前を書いておけば良い。

f = open('./input/input.txt')
f = open('input.txt')

など。

ファイルを出力するための名前を取得する方法

ファイルを書き出すのに名前を設定しておきたい場合ってのも結構ある。数少ないなら上のargvsの例のようにユーザーに名前を決めてもらうのも良い。
もちろんプログラム内に直に出力名を決めておいても良い。ただしこの場合プログラムを回すたびに上書きされるのでそれだけ注意が必要。

たくさんファイルをまとめて読み込んだ時に便利なのは、入力ファイル名に何かしら名前を足して出力ファイルの名前に変えてしまう方法。

以下はargus[1]で入力ファイル名を取得して、入力ファイルにsubtという文字を足して出力するという例で。Subtってのは何かしらデータファイル。

name = 'subt'+argvs[1]
np.savetxt(name, Subt)

これを例えばtkinterで大量に読み込んだファイル名に対して同じことをしてあげれば、入力ファイル名に次々と自動的にファイル名が割り振られていくことになる。

とここまで書いて、思い出した。tkinterの場合だとパス取得になってるから、このままだとおかしなことになるんだ。

basename = os.path.basename(filenames[i]+str(i)+'.txt')

こんな感じにos.path.basenameを使ってファイル名に変換してあげる必要がある。

今の所思い浮かんだ名前の取得方法はこれくらい。何か他のものを使う機会があれば後ほど追記していくことにする。

関連記事

pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使ったバックグラウンド強度の評価法いろいろ

色々と科学データを扱うとき、バックグラウンド強度の見積もりをすることは多い。たいていの場合は付属のソフトがついていて、バックグラウンドの計算をして減算してくれる機能がついている。
しかしそんな場合でも勝手に自動でやってくれるならいいのだけど、意外と1個1個マニュアルでやらなければならないこともある。10個くらいのデータならマニュアルでやっても良いのだけど、100個から1000個とかになってくるとプログラムを走らせてまとめて見積もりをさせたくなる。

まとめてファイルを読み込む方法はいろいろあるけれど、下記の投稿などを参照していただきたい。

内部リンク:pythonのtkinterを使ってユーザー入力を取得する方法

ということでこの投稿ではまとめて読み込んだデータを、まとめて処理するのに便利なバックグラウンド評価や減算方法についていくつか書いて行きたい。私は主にX線散乱に使っているけれど、ラマンなど他のスペクトルデータでもなんでも良い。

数式を使ったバックグラウンド計算

以下データはnumpyのndarrayに入った散布データということで。0列にX軸の値、そして1列に対応する強度データなどが入っているとする。下の例ではpwdという名前になっているが、なんでこんな名前をつけたかは忘れた。ndarrayへのデータの収納などに興味があれば下記投稿などで。

内部リンク:pythonでのテキストデータの入出力方法いろいろ

さてどうバックグラウンドを処理するかということだけど、一番都合が良いのはバックグラウンドが理由があるタイプのもので数式に従うはずのケース。理論式をそのまま使っても良いし、フィッティングをしても良い。
さてそれをどうpythonで引くかってことだけど、普通に式を作ってデータポイントごとに引くだけで良い。
まずは意外とよく使うコンスタントな最小値をバックグラウンドとして引く場合。これは最小値をnumpyのnanminで見つけて引くだけで良い。範囲指定(この例の場合XからY)の最小値や平均値を拾ってもいいし、どこかの定点の値を恣意的に引っ張ってしまっても大した問題はないだろう。

#XからYの最小値を探して各データポイントから引く
const = np.nanmin(pwd[X:Y,1])
pwd[:,1] = pwd[:,1]-const

続いて以下はデータをPower Lawでフィッティングをしてから引く場合の例。この例ではpower lawだけど、gaussianでもpolynomialでも好きな関数を使ってフィットすれば良い。フィッティングにはlmfitを使っている[1]。

内部リンク:pythonのlmfitを使ったカーブフィッティング

from lmfit.models import PowerLawModel

amp, exp = fit_PL(pwd)

def fit_PL(pwd):
#データファイルの44から70と574から584のエリアを使ってフィットしたかったのでnumpyのconcatenateで結合。astypeで変換してたのはデータファイル読み込みの都合だったと思う。たいていの場合なくて良いはず。
y = np.concatenate((pwd[44:70,1],pwd[574:584,1]),axis=0)
y = y.astype(float)
x = np.concatenate((pwd[44:70,0],pwd[574:584,0]),axis=0)
x = x.astype(float)

#lmfitのデフォルトのPowerLawModelを使ってフィット
PL_mod = PowerLawModel(prefix ='PL_')
pars = PL_mod.make_params()
pars['PL_amplitude'].set(10000)
pars['PL_exponent'].set(-1.2, min=-4, max = -1)
res = PL_mod.fit(y, pars, x=x)
print (res.fit_report())

#Power Lawのパラメータを取得
amp = res.best_values.get('PL_amplitude')
exp = res.best_values.get('PL_exponent')

#プロットするならmatplotlibなどで
plt.plot(pwd[:,0],pwd[:,1])
plt.plot(pwd[:,0],amp*(pwd[:,0]**exp))

#Power Lawパラメータをリターン
return amp, exp

この例の場合だと、データの最初と最後を使ってpower lawにフィットしたかったので、numpyのconcatenateを使ってそのようなデータセットを作ってからフィットした。
そしてamplitudeとexponentの二つのパラメータを関数から返す。

pwd[:,1] = pwd[:,1] - amp*(pwd[:,0]**ex)

この二つのパラメータを使って各データポイントの値を計算し、引き算をすればPower Lawの引き算がおしまい。こんな感じで式を使う場合は、簡単に引くことができる。

バックグラウンドがシンプルな場合のベースライン補正的なバックグラウンド演算

上記のような計算できるバックグラウンドならば良いのだけど、計算式によらないいろいろ混ざったバックグラウンドが引きたいってことも結構多い。以下はそんな時に使えるバックグラウンドの見積もり方について。こういった手法って結構いろいろ転がっていて、pythonだと割と簡単に導入できることが多い。

まずはバックグラウンドがどんどん下がって行くや、どんどんと上がっていくのようなシンプルなケースに使えるバックグラウンドの引き方。アイディアは少し古い1997年のクレイのX線研究の論文から[2]。たぶんオープンアクセス。
私がやったのは論文の手法と全く同じではないのかな。少しわかりにくかった場所があったのでそこらは適当に書いたので。結果はほぼ一緒になっていると思うけど。

大雑把にいうと、ある地点から直線を引く。直線の傾きを少しずつ変えていってデータのカーブと接したところでバックグラウンドラインとする。でその接点からまた直線でのバックグラウンド評価をやり直すということを繰り返す。基本的にはピーク以外のところがフラットになる操作。

装置的要因や、アモルファス的なバックグラウンドなどを考えて引くものではない。どんな実験したんだっていうような安定感のないデータ群を、とりあえずフラットなベースラインに揃えたい場合に使う感じ。
アモルファスなどの非ピーク情報を得たい場合にはあまりよろしくない手法。根元をがっつり削っちゃうので。
あとはスペクトルデータのベースライン補正などにも使えるかな。

#まずはスタート地点(init_ang)を指定する。conv_lineは値を行番号に変換する定義関数。意外と便利というかよく使う。
init_ang = 9.5
init_ang = conv_line(pwd,init_ang)
lipbkg = fit_lip(pwd,init_ang)

#バックグラウンド演算の関数fit_lip()。lipbkgにバックグラウンドのデータを詰めていく。
def fit_lip(pwd,init_ang):
  pwd = pwd.astype(float)
  plt.plot(pwd[:,0],pwd[:,1])
  step = 1
  init_loop = init_ang-1
  lipbkg = pwd[:,1]

#バックグラウンドを決める直線を引き直すためのループ。回数はデータ次第だけどある程度大きければ良いはず、適当。
  for j in range(100):
    init_loop = init_loop+step+1

#直線の開始地点がデータポイントの行数を超えたらループ終了。この場合582だったけどもちろんデータ次第。
    if init_loop > 582:
      print ('break')
      break
#バックグラウンドの直線とデータのカーブが接しているかチェックするループ。直線は単純に一次関数で定義して、傾きを10ずつ変えている。データ次第で刻みは適宜調整する必要がある。
    for i in range(40):
      slope = -(i*10)
      sec = pwd[init_loop,1] - slope*pwd[init_loop,0]
      linear = slope*pwd[init_loop:,0]+sec
      #plt.plot(pwd[init_loop:,0],linear,alpha=0.2)
      diff = pwd[init_loop:,1]-linear
      step = 1
#生データと直線の差分の最小値がプラスになったところで直線とカーブが接していると判断。argpartitionで何個目のデータポイントで接しているのかIndexを抜き出して、バックグラウンドの強度を得るのと次のループでのスタート地点を決めるのに利用。
      if np.nanmin(diff[1:]) > 0:
        step = np.argpartition(diff[1:],0)[0]
        print('i,j',i,j)
        lipbkg[init_loop:init_loop+step+2] = linear[0:step+2]
        break

  plt.plot(pwd[:,0],lipbkg[:])
  plt.show()

  return lipbkg

#データの値と行番号を変換する関数。差分の一番小さくなるIndexをargpartitionで抜き出すだけ。よく使うので比較的速度の速いargpartitionで。
def conv_line(plot,num):
  num = float(num)
  dif = np.absolute(plot-num)
  mini = np.argpartition(dif,0,axis=0)
  return mini[0][0]

というわけで、ある点から引いた直線がカーブと接するたびにその地点からバックグラウンド評価をやり直すことで、直線的なベースラインを得る手法でした。
評価範囲を限定しながら、直線の角度をもっとフレクシブルにすればデコボコしたバックグラウンドにも使えるだろう。しかしこのプログラムを書いた時は一様に減少するバックグラウンドだったのでこれだけで大丈夫だった。forループたくさんだし早くもないけれど、困るほど遅くもなかったと思う。

ピークvs滑らかなアモルファスカーブのスムージングを利用した分離法

こちらはシャープなピークとそうじゃないものをざっくばらんに分けたいという場合に有用。元々はX線のバックグラウンドで開発された手法だけれども[3]、スペクトロスコピーのベースラインの平滑化とかにも使える。
必ずしも正確にバックグラウンドを表現できる手法とはいえないかもしないが、いくつか利点がある。
比較的早い。プログラムの途中に組み込みやすい。評価手法次第だがある程度の定量性がある。演算の回数やスムージングの条件次第でブロードピークの除去から、スペクトルのベースラインを平滑化したい時まで利用できる。比較的複雑なバックグラウンドにも対応できる。などなど。
アイディアとしてはピークを含んだデータに対してスムージングを行うと、ピークは潰れるけどピークじゃないところは大体そのままになる。そしてスムージング後のデータとオリジナルデータを各データポイントごとに比べて値の小さな方を残す。これをどんどんと繰り返すことで、ピークが速やかに潰れていき最終的にバックグラウンドが残ることになる。弱点は大体どんな評価法でもそうだけどピークのオーバーラップ。

というわけで以下はスクリプト。オリジナルの論文とは違う論文のものを参考に、savgol_filterを利用したもの[4]。

#下に書いたbkg_estimate関数を使って最初にバックグラウンドを潰す。
bkg_init = bkg_estimate(Imin, Iave, pwd)
#bkg_iteration関数を一度回して1回目のバックグラウンドを得る。
BackG = bkg_iteration(bkg_init,pwd,lin1,lin2,win)
#後は好きなだけスムージングサイクルを回せば良い。
last = 599
for j in range(1,last+1):
  BackG = bkg_iteration(BackG,pwd,lin1,lin2,win)

#バックグラウンドの初期値を決める関数。のちのバックグラウンド減算の回数を減らせる効果があるのだけど、今のコンピュータだと誤差範囲。オリジナルの論文が書かれた2000年には効果があったのだろう。
def bkg_estimate(arr1, arr2, arr3):
  Imin = arr1
  Iave = arr2
  Int = copy.deepcopy(arr3[:,1])
  plt.plot(arr3[:,0],Int)
  for i in range(len(arr3[:,1])):
    if arr3[i,1] < (Iave +2*(Iave-Imin)):
      Int[i] = arr3[i,1]
    else:
      Int[i] = Iave +2*(Iave-Imin)

  plt.plot(arr3[:,0],Int)
  return Int

#ピークをスムージングでつぶしながら、バックグラウンドを計算する関数。
#引数はarr2は諸事情でプロット用に入れていたので、普通はなくても良い。それからlin1,lin2もデータ範囲の指定に使ってるだけなのでなくても良い。winは結構この引き方に致命的に効くパラメータなので、メインの関数から書き直せるようにしているので引数になっている。もちろん直接関数の中に書けば省略できる。
def bkg_iteration(arr1,arr2,lin1,lin2,win):
  bkg = arr1
  area = 0
#scipyのsavgol_filterを使ってピークをスムージングで潰す。
  Nbkg = savgol_filter(bkg,win+1,1,mode='nearest')
#あるデータ範囲でスムージング前後の小さい方の値をバックグランドとして残すことでピークがどんどん潰れて行く。
  for i in range(lin1,lin2):
    if (bkg[i]-Nbkg[i])<0:
      bkg[i] = bkg[i]
    else:
      bkg[i] = Nbkg[i]
  plt.plot(arr2[:,0], bkg)
  return bkg

状況によっては千回とか回すので、forループの中でいちいち保存するよりも自分自身(BackG)を書き換えた方が良いのかなとプログラムを書いたのだった。
しかしこういう書き方は絵面があまりよろしくない。どういう風に書くと美しいのだろうか?なんかPython的な書き方もありそうなものだけど。

いずれにしろこの手法の一番の問題点はいつスムージングを終わらせるかというところ。
特にバックグラウンドを使うわけでなければ目で見て好きなところで止めれば良いけれど、バックグラウンドを使いたい場合はある程度定量性のある止め方を考えないといけない。しかしこの投稿の本題ではないので深くは突っ込まない。

参照文献
[1] lmfit, https://dx.doi.org/10.5281/zenodo.11813
[2] Clays and Clay Minerals, 45, 132-146, 1997.
[3] Journal of applied chemistry. 33, 977-979, 2000.
[4] Carbohydrate polymers. 78, 543-548, 2009

関連記事

pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使って座標変換計算 〜回転行列を使って座標回転

座標計算はプログラムを書いておくと便利。 一度書いてしまえば使い回しもしやすいし。
無料のソフトウエアもあるけど融通が効かない場合もあるので。

ということでこの投稿ではpythonでの座標計算、主に回転行列を使って任意軸周りに座標を回転変換する方法について書いていきたい。

pythonでgenfromtxtを使った座標の読み込み

読み込みは色々とあるでしょうけれどご自由に。
xyz座標をfloatなりdoubleなりで取得しておく。 私は以下のような文字と数値の混ざった結晶pdbファイルからの読み込み。

ATOM   1  C4  N   1    -0.094  -0.316   0.079 1.000  2.63
ATOM   2  C3  N   1    -0.818   0.983  -0.299 1.000  2.63
ATOM   3  O3  N   1    -1.800   1.360   0.658 1.000  2.63

残りの文字列などを残したまま編集したかったのでgenfromtxtを使用してnumpyのarrayに投入。 ヘッダーが四行ついていたのでskip_header=4で読み飛ばす。
書き出す時ににヘッダーを元に戻したいならば、処理した座標を書き出す前に別にヘッダーを読み込んでファイルに書き込んでおけば良い。

orig = np.genfromtxt(argvs[1],dtype=['U8','int','U8','U8','int','float','float','float','float','float'],skip_header=4)

pythonでの座標の移動

とりあえず座標の移動は足し算でも引き算でもしとけば良い。
のだけれどgenfromtxtで違うデータ型を読み込んだ場合だと、少しややこしいことがある。
今上の例でorigをプリントすると(3,)というshapeを持つ、一見タプルを突っ込んだ三行のリストのように出力される。 でも実はこれはタプルではなくstructured arrayというもの。

[('ATOM', 1, 'C4','N', 1, -0.094, -0.316,  0.079, 1., 2.63)
 ('ATOM', 2, 'C3','N', 1, -0.818,  0.983, -0.299, 1., 2.63)
 ('ATOM', 3, 'O3','N', 1, -1.8  ,  1.36 ,  0.658, 1., 2.63)]

さてこれ一見タプルだけどタプルじゃない。 そこでそのまま値の変更が可能なのが良い点。
しかし同じデータ型で読み込んだ時とは異なり、普通のarrayではない。 そのためスライシングなどを使って、一気に値を変更する方法が見つけられなかった。

私の例の場合はデータの5列目から7列目の値を変更したかったのだ。 とりあえず下の例ではAdd[i][5]からAdd[i][7]に対してfor文を回して書き換えをしている。

Adcord = addcoords(orig, 6.8,1.8,42.0)
def addcoords(orig,num1,num2,num3):
 Add = copy.deepcopy(orig)
 for i in range(len(orig)):
 Add[i][5]=num1+Add[i][5]
 Add[i][6]=num2+Add[i][6]
 Add[i][7]=num3+Add[i][7]
 return Add

しかしこの投稿を書くのにいろいろ調べていたのだけど、こういう風にindexをつけてデータを抜き出してる例ってほとんどない。
というかgenfromtxtでデータ型の混ざったテキストを読み込んでそのまま処理する例自体があまり見当たらず。

きっともう少しスマートなやり方はあるのだろう。
だけどこのくらいならfor文1万回回したとしてもそんなに遅くはないのでとりあえず良しとしておく。

それからここでdeepcopyを使ってるのは、関数内で引数データに対して処理をして返すのに必要だったのだと思う。
これもベターな書き方があるんだろうけど、調べるのが面倒臭くてスキップしてしまっているものの一つ。

pythonで回転行列を使っての座標の回転

さてようやく本題の座標の回転へ。
私がやりたかったのはある2点を通る線を軸にしてある座標を回転させるということ。

まずは任意軸周りの回転行列の定義。 ロドリゲスの公式というのだったかな。 引数のAngが角度でnが任意軸の単位ベクトル。

def rotmat(Ang,n):
 R = np.array([[np.cos(Ang)+n[0]*n[0]*(1-np.cos(Ang)), n[0]*n[1]*(1-np.cos(Ang))-n[2]*np.sin(Ang), n[0]*n[2]*(1-np.cos(Ang))+n[1]*np.sin(Ang)],
[n[1]*n[0]*(1-np.cos(Ang))+n[2]*np.sin(Ang), np.cos(Ang)+n[1]*n[1]*(1-np.cos(Ang)), n[1]*n[2]*(1-np.cos(Ang))-n[0]*np.sin(Ang)],
[n[2]*n[0]*(1-np.cos(Ang))-n[1]*np.sin(Ang), n[2]*n[1]*(1-np.cos(Ang))+n[0]*np.sin(Ang),np.cos(Ang)+n[2]*n[2]*(1-np.cos(Ang))]])
 return R

上記の行列を使って回転させる関数の定義。 引数は座標と角度。 角度はどっかでラジアンに直しておく必要あり。
ちょっと大きめの座標だったのだけど、七十七行目と0行目の座標を軸(n1)に残りの座標を全部回転させている。 linalg.normでベクトルの大きさを出力し、それで割って単位行列にしてある。

for文の範囲が半分なのは座標の半分だけ回転させたかったから。 ここら辺は任意で。

一度全座標を軸ベクトルの始点で引いてから回転させている。 がちょっと前に作ったスクリプトなので、どうしてこうしていたのかは忘れてしまった。 こうしないとうまく動かなかったんだっけ?

回転自体はnp.dotを使って各座標に対して行列積を計算してあげるだけ。

def rotresi1(orig,thmo):
 new = copy.deepcopy(orig)
 Ang = (thmo*np.pi)/180
 n1 = np.array([new[77][5],new[77][6],new[77][7]])-np.array([new[0][5],new[0][6],new[0][7]])
 n1 = n1/np.linalg.norm(n1)
 ROT = rotmat(Ang,n1)
 for i in range(0,int((len(orig)/2))):
  Sc = np.array([orig[0][5],orig[0][6],orig[0][7]])
  Coord = np.array([orig[i][5],orig[i][6],orig[i][7]])
  Coord = Coord -Sc 
  new[i][5]=np.dot(ROT,Coord)[0]+orig[0][5]
  new[i][6]=np.dot(ROT,Coord)[1]+orig[0][6]
  new[i][7]=np.dot(ROT,Coord)[2]+orig[0][7]
return new

んでたくさん違う角度で回した座標とかを一気に吐き出したかったら、下のようにfor文を回してあげれば良い。 for文は避けても良いけれど、このくらいの計算量ならまいっかなってくらいの速度。
データ数が莫大ならばちょっと工夫した方が良いかもしれない。

最後は回転した座標をheaderを一行だけつけてフォーマットを揃えてsavetxtで出力。 ここら辺も任意で。

for i in range(36):
 rot = 10*(i+1)
 New = rotresi1(orig,rot)
 np.savetxt('adj'+str(rot)+'.pdb',New,fmt='%4s %6d %3s  %-4s  %2d %11.3f %7.3f %7.3f %4.3f %5.2f',header = 'CRYST1    5.130   19.280   10.660  90.00  90.00  97.00',comments='')

というわけでpythonを使った座標計算の一例でした。
最近こういう計算から離れてしまっているので、もう少し取り組んでいきたいところ。

関連記事

pythonのまとめ

D