Pythonを使って散布図を台形積分

初めに

卒業研究でよく散布図を用いるのですが、散布図の積分って面倒なんです。プロット数が多いので尚更そうなんですが骨が折れるので、今回は散布図で台形積分するプログラムをpythonを使って書きました。

調べても意外と載ってないのでここにメモしておきます。 プログラム全体の説明はせず、部分的に説明していきます。 ついでにpythonでデータ整理をしたいときに便利な機能も掻い摘んで説明していきます。

※今回ご紹介する以外にも今までにpythonを用いたプログラムはいくつか書いてきましたが、体系的に勉強したことがない筆者は蛇足なプログラムが多々あると思うのでご理解ください。アドバイスを頂けるとありがたいです。

積分するプロットについて

今回、積分しようとしているのはヒステリシスループで囲まれた面積です。ヒステリシス損失と言われるもので、変圧器などで使われる軟磁性体でよく議論されているような気がします。

積分の方法は省略して、以降プログラムの各種機能・使い方を説明します。

テキストファイルを読み込む

with open(name,'r') as f :
    source = f.read()

  • open関数内にてファイル名を「name」で指定してmodeは読み取り専用モードの「r」としている。

  • 「with」を使うことでファイルを「close」する必要がない。また、「as f」とすれば「f」で指定したファイルをプログラム内で呼び出せる。

  • 「f.read()」はstr型で読み込んでsourceに格納している。

文字・数字を成形する

source = f.read().split()
  • 先ほど格納したsourceを「.split()」を使うことでスペース・改行・タブで区切ってlist型にできる。
source = f.read().split(",",2)
  • 第一引数に指定した文字で区切ったり第二引数で最大分割回数を設定できたりする。

２つ以上の指定文字で区切りたい場合

「reライブラリ」を使用する方法。詳しくは省略するがこれを利用すると特定の文字で自由自在に区切ったり省いたりできる。下に「+-。」で区切る例を載せている。

re.split('[+-。]',source)

リストを結合して文字列に戻す方法

「join」を利用。下にlist型の「source」を各要素の間に「|」(縦棒)を入れて結合した後、文字列を出力する例を載せている。

'|'.join(source)

リスト型の簡単な成形方法

データの取り出しと結合

source_y = []
source_y.extend(source[10::3])
  • 空のlistであるsource_yに「.extend」でlistの後ろから要素を追加していく。演算子「+」でも同じことができる。
  • 追加する際に「source[10::3]」とすることでlist型の「source」の10個目の要素から3つ飛ばしで最後まで追加するように指定している。
source = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]
source_y.extend(source[9:21:3]) 

>>> source_y = [9,12,15,18]
  • 「source[始まるindex:終わるindex:飛ばす要素数]」となっている。ちなみに上のプログラムでもわかるように「9以上21未満」となるので注意。

データの演算

source_y_2= [0]+source_y                  
source_y = source_y + [0]
  • 「.extend」と同様に使える。今回は積分をするので「[0]」をlistに加えることで意図的にindexをずらして後々の計算が簡単にできるようにしている。
  • 「[0]」の演算を「source_y」の前に加えるとそれに続く全ての要素のindexが+1ずれる。後に加えるとlistの最後に加えられることになるためindexは変化せず、要素の数を増やすことができる。
  • ２つのlist「source_y」と「source_y_2」の要素数が合わないと後々の計算ができないため両方に「[0]」を加えている。

numpyを利用する

numpyはpythonの行列計算ライブラリ。今回は積分する際にlistの要素ごとに計算が必要になるため利用した。numpyを利用することで要素ごとの四則演算を簡単に行うことができる。

listからndarrayへの変換方法

import numpy as np
arr_y1= np.array(source_y,dtype='float32') 
arr_y2= np.array(source_y_2,dtype='float32')
  • 「np.array」を用いてlist型「source_y」をnumpyで利用できる「ndarray」型へ変換している。
  • 「dtype='float32'」でfloat型のデータとして変換するように指定している。
  • 整数のint型や桁数を増やしたfloat64なども指定できる。

ndarray内で最大値・最小値のindexを得る

index_y1 = np.argmax(arr_y1)
index_y2 = np.argmin(arr_y1)
  • 「index_y1」にndarray内の最大値のindex、「index_y2」にndarray内の最小値のindexを代入している。

ndarray内で連続データの符号が反転するindexを得る

mult = np.sign(np.multiply(arr_y1,arr_y2))
hanten =np.where(mult<0)
print(hanten)

>>> array([132, 225], dtype=int64)

先に用意しておいた積分用にindexを１つずらしたlistと元のlistを掛け合わせることで符号が反転するindexを特定する。

  • 「np.multiply(arr_y1,arr_y2)」で２つのndarrayの要素ごとの積をとっている。
  • 「np.sign」で符号のみに着目してマイナス値を「-1」、0を「0」、プラス値を「1」としたndarray「mult」を生成する。
  • 「np.where」でmultが0より小さくなるindexをすべて得る。

出力結果を見るとわかるがこのやり方だと計算用にindexの数字だけを取り出すのが難しい。そこでこの値はプロットに異常がないかの確認用とすることにした。

index_y =np.argmin(mult)

先ほど紹介した最小値のindexを得るプログラムを利用した。
しかし、このようにすると最初の「-1」のindexしか得られない。

index_y =np.argmin(mult)
mult = mult[index_y+1:]
index_y2 = np.argmin(mult)+index_y+1

そこで、確認済みのindexまでをすべて削ったndarrayを新たに作ることで次の「-1」のindexを無理やり得られるようにした。
（これを応用してfor文でndarrayを削っていって全ての符号反転indexを得ることができるプログラムが書けそうですが今回は２つだけで十分なので書いてません。）

numpyで絶対値と四則演算

arr_y1 = np.abs(arr_y1)

「np.abs」でarr_y1の絶対値をとる。

com_y = np.add(arr_y2,arr_y1)
dis_y = np.subtract(arr_y2,arr_y1) 
mul_yy = np.multiply(com_y,dis_y)
div_yy = np.divide(com_y,dis_y)
delta = np.abs(mul_yy/2)

今回の積分は台形積分のため割り算を除いて要素ごとに四則演算が必要となる。

  • 「np.add」足し算
  • 「np.subtract」引き算
  • 「np.multiply」掛け算
  • 「np.divide」割り算

最後の「mul_yy/2」は全ての要素を2で割っている。

ndarrayの一部データ削除と合計

delta = np.delete(delta,0,0)
delta = np.delete(delta,delta.size-1,0)

行列である微小面積「delta」は最初と最後だけ「0」を足した影響で値が大きくなるので省く。

  • 「np.delete」を利用する。
  • 「np.delete(行列名,行番号,列番号)」で指定。行番号と列番号は0スタートなので注意。
  • 「delta.size」でdeltaの要素の数を取得し0スタートを考慮して「delta.size-1」すると最後の要素を消せる。
summ = delta[0:100].sum()
  • 「.sum()」で行列の合計をy1に代入している。
  • 「delta[0:100]」とすることで0から100番目の要素の合計を出せる。

テキストファイルを出力する

 with open("面積.txt",'a') as f_out :
    f_out.writelines("\n"+"面積：" +str(summ) +"\n")
    f_out.close()

open関数内にてファイル名を「面積.txt」としてmodeは書き込みモードの「a」としている。ここで書き込みモードに関して説明する。

  • 「a」：ファイルがない場合は新規作成。ある場合は続きに書き足される。

  • 「w」：ファイルがない場合は新規作成。ある場合は上書きされる。

  • 「x」：ファイルがない場合は新規作成。ある場合はエラー。

最後に「f_out.close()」でファイルを閉じている。必要ないかもしれない・・・。

ループ処理・待機・例外処理

try:
    while True:

        print("ENTERキーを押すと繰り返し計算できます。")
        key = input("半角SPACEを入力すると終了します。")
        if key == ' ':
            break

except:
    print("error")

  • 「try」と「except」で例外処理を行う。例外の場合は「error」と表示される。

  • 「while True:」でループ処理が可能。

  • 「input」を利用することで待機ができる。

  • 「input」で特定の入力を読み込み、プログラムを終了させることができる。

実行ファイル「.exe」に変換

pip install pyinstaller

「pyinstaller」を利用する。ライブラリを入れていない人はpipでインストール。

pyinstaller hello.py --onefile

コマンドプロンプトなどで階層を合わせて上記のコマンドを実行すると実行ファイルがdistフォルダの中にできる。 ファイルが多数生成されるため、実行ファイルを作る前にプログラム本体をフォルダの中に入れておくのがおすすめ。

型の確認方法

source = []
print(type(source))

>>> list

型のエラーはよく出るのでこれで確認するのがおすすめ。

終わりに

ここまで読んでいただきありがとうございます。メモ書きということで、最初から最後まで分かりにくく見にくい内容になっています。掻い摘んで色々と説明しているためまとまりもないのですがご了承ください。また、今後も修正を重ねていく予定です。