Style File と Data File を分離した QDP-like な ROOT での作図

単なる遊びです。Style file にデータの表示形式を書き、データ辞退は Data file に書いて ROOT で図を作る簡易スクリプト。

#!/usr/bin/env python
# Teruaki Enoto 2012-02-03
from optparse import OptionParser
import array
import shlex
import ROOT

stylefile_def = "/Users/enoto/work/program/root/example/PlotXY/stylefile.txt"

class PlotXY():
    def __init__(self, data_array, stylefile, outputpdf):
        self.data_array = data_array
        self.stylefile  = stylefile
        self.outputpdf = outputpdf

    def show(self):
        print "data_array = ", self.data_array
        print "stylefile = ", self.stylefile

    def init_par(self):
        ### initialize parameters 
        self.pars = {"Title" : "str",
                     "Xmin" : "float",
                     "Xmax" : "float",
                     "Xtitle" : "str",
                     "Xlog" : "bool", 
                     "Ymin" : "float",
                     "Ymax" : "float",
                     "Ytitle" : "str",
                     "Ylog" : "bool", 
                     "Format" : "list", 
                     "Style" : " list",
                     "Can_X0" : "int",
                     "Can_Y0" : "int",
                     "Can_Xwid" : "int",
                     "Can_Ywid" : "int",
                     "TitleBorderSize" : "int",
                     "TitleAlign" : "int",
                     "TextFont" : "int",
                     "FrameLineWidth" : "int",
                     "EndErrorSize" : "float"
                     }
        for par in self.pars:
            exec("self.%s = %s('%s')" % (par, self.pars[par], '0'))

    def set_style(self):
        ### open file 
        try:
            f = open(self.stylefile)
        except IOError, err:
            print 'Error: cannot read stylefile (%s):' % (self.stylefile,), err
            quit()

        ### get file content 
        lex = shlex.shlex(f,posix=True)
        lex.whitespace += "="
        lex.whitespace_split = True

        ### read and set parameters
        while True:
            s = lex.get_token()
            if s!=lex.eof:
                val = lex.get_token()
                print s, val
                if s == "Style" or s == "Format":
                    cmd = "self.%s = %s(%s)" % (s, self.pars[s], val)
                elif s == "Xlog" or s == "Ylog":
                    if val == "True":
                        cmd = "self.%s = %s(%s)" % (s, self.pars[s], True)
                    elif val == "False":
                        cmd = "self.%s = %s(%s)" % (s, self.pars[s], False)
                    else:
                        print "Error:: Xlog/Ylog must be bool (True or False)"
                else:
                    cmd = "self.%s = %s('%s')" % (s, self.pars[s], val)
                exec(cmd)
            elif s==lex.eof:
                break
            else:
                print "Error: unexpected error in reading the parfile."
        f.close()

    def read_all(self):        
        self.graph = []
        id = 0 
        for file in self.data_array:
            if self.Format[id] == 0:
                x, y = self.read(file)
                self.graph.append(ROOT.TGraph(len(x), x, y))

            elif self.Format[id] == 1:
                x, y, xe, ye = self.read_err(file)
                self.graph.append(ROOT.TGraphErrors(len(x),
                                                    x, y,
                                                    xe, ye))

            elif self.Format[id] == 2:
                x, y, xeh, xel, yeh, yel = self.read_asymerr(file)
                self.graph.append(ROOT.TGraphAsymmErrors(len(x),
                                                         x, y,
                                                         xeh, xel,
                                                         yeh, yel))
                
            else:
                print "Error:: Format error in style file"
                quit()
            id += 1

    def read(self, file):
        print "...reading file %s..." % file
        x = array.array('f')
        y = array.array('f')
        for line in open(file):
            col = line.split()
            x.append(float(col[0]))
            y.append(float(col[1]))
        return x, y

    def read_err(self, file):
        print "...reading file %s..." % file
        x = array.array('f')
        y = array.array('f')
        xe = array.array('f')
        ye = array.array('f')
        for line in open(file):
            col = line.split()
            x.append(float(col[0]))
            y.append(float(col[1]))
            xe.append(float(col[2]))
            ye.append(float(col[3]))
        return x, y, xe, ye

    def read_asymerr(self, file):
        print "...reading file %s..." % file
        x = array.array('f')
        y = array.array('f')
        xeh = array.array('f')
        xel = array.array('f')
        yeh = array.array('f')
        yel = array.array('f')
        for line in open(file):
            col = line.split()
            x.append(float(col[0]))
            y.append(float(col[1]))
            xeh.append(float(col[2]))
            xel.append(float(col[3]))
            yeh.append(float(col[4]))
            yel.append(float(col[5]))
        return x, y, xeh, xel, yeh, yel

    def plot(self):
        can = ROOT.TCanvas("c1","A Simple Graph with error bars",
                           self.Can_X0, self.Can_Y0, self.Can_Xwid, self.Can_Ywid)
        print self.Xlog, self.Ylog
        can.SetLogx(self.Xlog)
        can.SetLogy(self.Ylog)

        # --- Optional Gstyles ---
        ROOT.gStyle.SetTitleBorderSize(self.TitleBorderSize)
        ROOT.gStyle.SetTitleAlign(self.TitleAlign)
        ROOT.gStyle.SetTextFont(self.TextFont)
        ROOT.gStyle.SetFrameLineWidth(self.FrameLineWidth)
        ROOT.gStyle.SetEndErrorSize(self.EndErrorSize)

        id = 0
        for graph in self.graph:
            if id == 0:
                graph.Draw("AP")
            else:
                graph.Draw("P")

            graph.SetMarkerColor(self.Style[id][0])
            graph.SetLineColor(self.Style[id][0])
            graph.SetMarkerStyle(self.Style[id][1])  
            graph.SetMarkerSize(self.Style[id][2])
            id += 1

        # --- Range & Title Set ---
        self.graph[0].SetTitle(self.Title)
        self.graph[0].GetXaxis().SetLimits(self.Xmin, self.Xmax)
        self.graph[0].GetXaxis().SetTitle(self.Xtitle)
        self.graph[0].SetMinimum(self.Ymin)
        self.graph[0].SetMaximum(self.Ymax)
        self.graph[0].GetYaxis().SetTitle(self.Ytitle)

        can.Print(self.outputpdf)
        
    def run(self):
        self.show()
        self.init_par()
        self.set_style()
        self.read_all()
        self.plot()

#############################
if __name__=="__main__":
#############################

    # --- Interface to command line
    parser = OptionParser()
    parser.add_option("-s", "--stylefile", dest="stylefile", default=stylefile_def, 
                      type="string", action="store", help="style file name")
    parser.add_option("-o", "--outputpdf", dest="outputpdf", default="out.pdf",
                      type="string", action="store", help="output pdf file name")
    (opt, arg) = parser.parse_args()

    # --- run main ---
    if len(arg) == 0:
        print "Error:: No input data_array."
        quit()
    else:
        plot = PlotXY(arg, stylefile=opt.stylefile, outputpdf=opt.outputpdf)
        plot.run()
    quit()

ほんで、Style file

Title = "Title"
Xmin = 0.8
Xmax = 10.0
Xtitle = "Xtitle"
Xlog = True
Ymin = 0.1
Ymax = 4.0
Ytitle = "Ytitle"
Ylog = True
#-------- Format file ---------
# Format 
# 0 : xy
# 1 : xy with x, y errors
# 2 : xy with x, y asimerr
# Style [color, marker, size]
Format = [2,1]
Style = [[2,20,1.1],[3,22,1.1]]
#--------- Other Format ----------
Can_X0 = 200
Can_Y0 = 10
Can_Xwid = 580
Can_Ywid = 450
TitleBorderSize = 0
TitleAlign = 13
TextFont = 132
FrameLineWidth = 1
EndErrorSize = 0.8

Data file は

3.0 2.0 0.2 0.4 0.2 0.1 
5.0 2.0 0.1 0.8 0.1 0.2 
3.4 2.5 0.1 0.7 0.2 0.1
5.6 0.3 0.1 0.8 0.3 0.2

で、図

パラメータファイルを読み込ませるスクリプトの定型

次のようなパラメータファイルを読むスクリプトの定型をメモしておきます。

par1 = 3
par2 = 0.4

from optparse import OptionParser
import shlex

class MainClass():
    def __init__(self, parfile, verbose=0, logfile="logfile.txt"):
        self.parfile = parfile
        self.verbose = verbose
        self.logfile = logfile

    def show_input(self):
        print "parfile = %s" % self.parfile
        print "verbose = %d" % self.verbose
        print "logfile = %s" % self.logfile

    def init_par(self):
        ### initialize parameters 
        self.pars = {"par1" : "int",
                     "par2" : "float"
                     }
        for par in self.pars:
            exec("self.%s = %s('%s')" % (par, self.pars[par], '0'))

    def read_par(self):

        if not hasattr(self, "pars"):
            print 'Error: init_par should be defined.'
            quit()

        ### open file 
        try:
            f = open(self.parfile)
        except IOError, err:
            print 'Error: cannot read parfile (%s):' % (self.parfile,), err
            quit()

        ### get file content 
        lex = shlex.shlex(f,posix=True)
        lex.whitespace += "="
        lex.whitespace_split = True

        ### read and set parameters
        while True:
            s = lex.get_token()
            if s!=lex.eof and self.pars.has_key(s):
                val = lex.get_token()
                cmd = "self.%s = %s('%s')" % (s, self.pars[s], val)
                if self.verbose > 0: print cmd
                exec(cmd)
            elif s==lex.eof:
                break
            else:
                print "Error: unexpected error in reading the parfile."
        f.close()

    def run(self):
        print "---run---"
        self.show_input()
        self.init_par()
        self.read_par()
        print "---done---"

#############################
if __name__=="__main__":
#############################

    # --- Interface to command line
    parser = OptionParser()
    parser.add_option("-p", "--parfile", dest="parfile", default="par.txt",
                      type="string", action="store", help="parameter file name")
    parser.add_option("-v", "--verbose", dest="verbose", default=0,
                      type="int", action="store", help="verbose level: int")
    parser.add_option("-o", "--logfile", dest="logfile", default="logfile.txt",
                      type="string", action="store", help="logfile name")
    (opt, arg) = parser.parse_args()

    # --- run main ---
    obj = MainClass(opt.parfile, verbose=opt.verbose, logfile=opt.logfile)
    obj.run()
    quit()

カイ二乗分布を描く。

特に何ということもないんだけれど、自由度 k=４のカイ二乗分布を ROOT と python で描いてみた。自由度 k=4 のカイ二乗分布は、平均値 k = 4、分散 2k = 8、最頻値 k - 2 = 2、中央値 k-2/3+4/27k -8/729k^2=3.4 で分布する。周期性を探す Zn-test (Raileigh Test) でよくお世話になる。

#!/usr/bin/env python

import ROOT
import pylab
from matplotlib import rc

def drange(start, stop, step):
    r = start
    while r < stop:
        yield r
        r += step

x = []
y_pdf = []
y_cdf = []
for i in drange(0,15,0.1):
    x.append(i)
    y_pdf.append(ROOT.Math.chisquared_pdf(i, 4))
    y_cdf.append(ROOT.Math.chisquared_cdf(i, 4))

rc('text', usetex=True)
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica'], 
             'size':'14'})

fig = pylab.figure(figsize=(7,7))
ax = fig.add_subplot(111)

ax.set_xlabel('Chi-square value')
ax.set_ylabel('Probability')
ax.set_title('Chi-square Distribution (Degree of Freedom = 4)')
ax.plot(x, y_pdf, 'r-', x, y_cdf, 'b--')
ax.legend(('Probability Distribution', 'Cumulative Distribution'), 'center right', shadow=False)

fig.savefig("chi2_dof4.pdf")

Palm Road 正面の写真

もう一個ついでに、Palm Road 正面を見た写真。

Stanford Green Library

Stanford 大学の Green Library。メイン図書館。
真面目に勉強している人がいっぱいる素敵なところ。
図書館の前には Cafe があるのが、なおよい。

Arduino と SHT71 センサーで実験室の温湿度ログを Mac から取得する。

実験室で温度・湿度を記録する、あるいは長時間でのログを取りたいがあります。この目的に合うStrawberry Linux 社の温度・湿度計モジュール USBRH が市販されていますが、Windows の driver しか公開されおらず、有志の開発したソフトで Linux でまでなら対応しているようです(先輩の例)。ただ、Mac User なので、できれば直接 Mac から制御したいところです。

そこで、Arduino とSENSIRION 社製のデジタル温湿度センサー SHT71 を使って、直接 Mac でデータ収集してみました。SHT シリーズは比較的高性能で、上記のモジュールにも使われています。今回は、PIN の足がついている SHT71 を用います。日本国内だと Strawberry Linux からも SHT71 や SHT11 を購入できます。自分は、アメリカで Newark 経由で購入して $33 でした。

SHT71 は、公称値で、 25℃ において湿度の精度が±3.0%、温度の精度が±0.4℃となっています。 詳しくは、SHT71 dataheet (pdf) をご覧ください。この SHT71 はデジタル値で出力してくれるので、PIN の足は順番に

1. SCL : Serial Clock, input only (同期用)
2. VDD : Source Voltage
3. GND : Ground
4. DATA : Serial Data, bidirectional (シリアルデータ通信)

となっています。

Arduino との配線は、10 kOhm のプルアップ抵抗をつけて、以下のようにしました。

肝心のスケッチですが、ライブラリを利用できます。SHT1X や SHT7X のライブラリが、サイトから sht1x.tgz としてダウンロードできます。あるいは、GitHub にもあるようです(ref)。ライブラリは以下に収容します。
~enoto/Documents/Arduino/libraries/sht1x

なお、Arduino Uno を使っており、sht1x_hw.h の中身で

# define SHT1X_CLOCK_LINE PD2
# define SHT1X_DATA_LINE  PD3

という箇所を

# define SHT1X_CLOCK_LINE PORTD2
# define SHT1X_DATA_LINE  PORTD3

のように修正しないとコンパイルが通りません。

その上で、スケッチは以下の通り。データシートにある補正式のままです。

#include <sht1x.h>
#include <sht1x_hw.h>

unsigned int h;
unsigned int t;
double temp;
double hlin;
double htrue;

void setup()
{
  Serial.begin(9600);
  delay(15);
  sht1x_init();
  Serial.println("Start");

  if (sht1x_read_humidity() < 0) {
    Serial.println("Reset");
    sht1x_reset();
  }
}

void loop()
{
  h = sht1x_read_humidity();
  t = sht1x_read_temp();

  while (millis() % 1000);

  temp  = -39.66 + 0.01 * t; // 14 bit
  hlin  = -2.0468 + 0.0367 * h - 1.5955e-6 * h * h; // 12 bit
  htrue =  (temp - 25.0) * (0.01 + 8e-5 * h) + hlin; // 12 bit

  Serial.print(temp, 3);
  Serial.print("C ");
  Serial.print(htrue, 3);
  Serial.println("% ");
}

大体、仕様の誤差範囲で、正しい温度のようですが、較正はまた確認してみようと思います。

参考サイト

• Arduino と Ruby で温湿ネットロガー
• 温湿度センサーSHT71のライブラリとエラー対策
• 2線式 温・湿度センサSHT11 アーカイブ (エレキジャック)
• センサー道場のプルアップ・プルダウン抵抗
• Arduinoで湿度・温度を計測する SHT-71使用
• Sensirion Humidity and Temperature Sensors (Arduino 本家)

Palm Road のヤシの木

Stanford University のメインストリート Palm Road のヤシの木。自転車やバスでこの横を通り抜ける。