前回の宿題
"merged .sca file from HKL2000"を読み込んで異常分散差を表示する。
output.sca処理の対象となるファイル情報
4行目以降の各行（1行に書いてある）には以下の情報が含まれている H,K,L,I(+),sigI(+),I(-),sigI(-) ただし、(-)が無い場合には情報は「ファイルに書き出されていない」。 I(+)が無い場合、I(+)=0, sigI(+)=-1が入るみたい。
プログラム作成
- コマンドラインからファイル名を指定
- ファイルを開いて行を読み込んでいく
- 3行目にある格子定数を文字列でも何でもよいから取り込んで後で表示でき るようにしておく
- 4行目以降では読み込んだカラム数が５の場合、異常分散差は計算できない ので無視、カラム数が７の場合、以下の式を計算していく。
delano=|I(+)-I(-)| / |sigI(+)とsigI(-)の平均値| （ただし、分母が0の場合は無視します） （調べているひまが無かったので式については別の機会に議論することにします）
- delanoという値を反射ファイルが続く限り足し続け、最終的に平均値を表示 する
宿題の答えと解説
#!/usr/bin/python import sys import math # check a number of arguments if len(sys.argv)!=2: print "Usage: readAnoSca SCAFILE" sys.exit() else: print "Processing %s"%sys.argv[1] # File open infile=open(sys.argv[1],'r') # reading all of lines linelist=infile.readlines() # initialization of parameters nline=0 # line counter sum_diffi=0.0 # summation of |I(+)-I(-)| sum_avesig=0.0 # summation of average(sigI(+),sigI(-)) nok=0 # number of reflections with anomalous # acquiring cell parameters at the 3rd line cell_string=linelist[2] # treating each line for line in linelist: # counting a number of line nline+=1 # initialization strings=line.split() # counting a number of columns ncolumn=len(strings) # data processing # ncolumn=7: anomalous signals can be obtained # ncolumn=5: anomalous signals cannot be obtained if ncolumn==7 and nline>4: # flag for judging is true at the first stage iflag=1 # for I(+)s i_plus=float(strings[3]) sig_plus=float(strings[4]) # case when I(+) was not observed -> skip this reflection if i_plus==0 and sig_plus==-1.0: continue # for I(-)s i_minus=float(strings[5]) sig_minus=float(strings[6]) # difference between I(+) and I(-) if iflag != -1: nok+=1 sum_diffi += math.fabs(i_plus-i_minus) sum_avesig += (sig_plus+sig_minus)/2.0 # Final calculation of average ave_diff=sum_diffi/nok ave_sig=sum_avesig/nok print "Process has been finished." print "<d''>/<sigI>=%8.3f"%(ave_diff/ave_sig)
つこてる関数
sys.argv
コマンドラインでプログラムを実行したときにその引数を取得する場合に利用。 プログラム名をttttとした場合、以下の例を見ましょう。
% tttt output.sca scale_out.plot
- 引数の数 len(sys.argv)は3である⇒len関数は配列の長さを返す
- 引数文字列はそれぞれ
- tttt -> sys.argv[0] コレに注意。
- output.sca -> sys.argv[1]
- scale_out.plot -> sys.argv[2]
open ファイルオープン
ファイルを開ける。テキスト・バイナリ問わずファイル処理を行なう場合はまずやること。
infile=open(sys.argv[1],'r')
infile: ファイル変数, open()：ファイルオープン関数
open(ファイル名,'オープンのモード')
ファイル名：文字通り開けたいファイルの名前
オープンモード: 読み込み、書き込み、テキスト、バイナリなどの指定をここで行う
| r | 書き込み専用 |
| w | 書き込み専用 |
| a | 追加(ファイルを開く時、ファイルの参照点は自動的にファイルの末尾に移動) |
| r+ | 更新(読み取りと書き込み) |
| w+ | ファイルを空にした後、そのファイルを読み書き用に開く |
| a+ | ファイルを読み書き用に開く、参照点は自動的にファイルの末尾に移動 |
| b | 他のオプションと一緒に指定すると、ファイルをバイナリモードで開く |
readlines （ファイルオブジェクト関数）
全ての行をリストへ格納。テキストファイルに対して適用する関数（だと思う）。 改行文字によって区切られた文字列のリストを返すので文字列リスト変数で受けておく。例だと、linelistで受けている。pythonは明示的に変数を定義しないので分かりにくいが、ここでのlinelistは
linelist=[]
という定義が抜けている。テキストファイルを処理するなら、open->readlinesという流れは定石のように思う（間違ってたらご指摘下さい）。
for文
for line in linelist:
この文章の意味は、linelist配列の中身がある限り、一つずつの要素をlineに入れて以下の処理をしてくださいってこと。このプログラムに限って説明すると１行ずつ読み込んで以下の処理をしてくださいってことになる。ちなみに行のインデントによってforやif文が区切られるのでこのfor制御は# Final calculation of averageの直前までかかっている。
例えばC言語などに使用されるようなforループは以下のように書く。
for i in range(1,15): print i
range関数によって指定された範囲にある数(1〜15)を順にiに入れて、それぞれのiに対して処理を行なう。上記例だと、printが処理部分に相当する。
文字列のカラム数を数えるには？ 文字列のsplit()
上記、.scaファイルのアノマラス信号が書いてある反射、そうでない反射は行にあるカラムの数で判定をしている。そのため、まずは行に書いてあるカラムの数を数える必要がある。
このとき役に立つのが文字列(string)のsplit()関数。
この関数は入力文字列を空白によって区切り、区切られた文字列を配列としてリストにキープする、というシロモノ。区切りに使う文字は
line.split()
とすれば、空白だし、
line.split(",")
とすればコンマになるので、適宜使うように。ただし気をつけないといけないのは、スペースによって区切られていない文字列はもちろんスプリットされない」ということである。例えば、以下のようなファイルがあったとしたら
13 34 34.3 98.8 14 39 32.9102.8
人間だったら下の行の最後のカラムと最後から２番目のカラムをかろうじて分けることが出来るかも知れない。しかし、プログラムは「スペース文字で区切ってください」といわれて正直にその処理のみを実行するわけだから、下の行で102.8というのをちゃんと認識することは出来ない。
さて本題に戻る。１行をsplit()したら、そのリストに含まれる文字列の数を数えたら、カラムの数に相当するはずなので、
ncolumn=len(strings)
としている。
文字列を浮動小数へ変換する
これまでの処理で.scaファイル中反射情報を構成する文字列を取得できたので、それらを数値として処理しましょう。
i_plus=float(strings[3])
この行で代表して書きますが、文字列のsplit()関数でゲットしたものはやはり文字列なので、数値として以下のような式を書くことが出来ないです。
i_plus=strings[3]/2.0
これをやるとstrings[3]は文字列なので式が成立しない。こいつを数値として扱いたければまずは変換が必要。
そこで出てくるのが float()関数。いわずもがなで文字列を浮動少数に変換することが出来ます。使うとしたらあとint()で整数型への変換が可能です。
