Numpyの活用

Matplotlibの基本コード

Scipyの活用

scikit-learnの活用

tensolflowの活用

Python

化学工学、数値計算に役立つTipsを紹介します。

まだコンテンツは少ないですが、Scipy、Numpy、Matplotlibの基本的な使い方に関するTipsを順次追加しています。

事前準備

※2017.06.10追記:なにはさておきAnacondaPycharmをインストールしてください。numpyやscipyのインストールで悩む必要がなくなります。

個人的にPythonが有難いのは以下の点です。

①行列計算や多次元配列の取扱いが簡単(numpy)

②Scipyの豊富な計算ライブラリ(最適化関数やODEのSolver等)

③グラフ描画が比較的容易で綺麗に描ける(Matplotlib)

④複素数計算が可能(Bode線図やNyquist線図を作る際に便利)

⑤外部ファイルの入出力が行いやすい

⑥変数宣言をせずに使える(コードが短く書ける)

Pythonでの数値計算の恩恵を受けるためにいくつか便利なライブラリを以下に紹介します。

・SciPy、NumPy、matpoltlib:最適化計算や逆行列演算、グラフ描画等に便利なライブラリです。

・pandas:データ処理に便利なライブラリです。

・PySide:GUIの作成に便利なライブラリです。

ちなみに私はインストールで何度もはじかれて苦戦しました。

インストールについてはゆくゆく整理していきたいと思います。

※2016.05.08追記:Scipy、Numpy、Matplotlibのインストールについては、こちらがよくまとまっています。

Numpyの活用

※Numpyのインストールが必要です。

・配列(array)、逆行列(linalg.inv)、転置(transpose)、積(dot)

・3×3行列の例

$$x=\begin{bmatrix} a & b & c\\\ d & e & f \\\ g & h & i \end{bmatrix}$$
x=np.array([[a,b,c],[d,e,f],[g,h,i]])

・上記3×3行列の逆行列

$$x=\begin{bmatrix} a & d & g\\\ b & e & h \\\ c & f & i \end{bmatrix}$$
x=np.linalg.inv(np.array([[a,b,c],[d,e,f],[g,h,i]]))

・3×1行列の例

$$a=\begin{bmatrix} 1 \\\ 2 \\\ 3 \end{bmatrix}$$
a=np.array([[1],[2],[3]])

また上記の行列を転置した以下の行列はa=np.transpose(a)のように書けます。

$$a=\begin{bmatrix} 1 & 2 & 3 \end{bmatrix}$$

行列同士の積は、例えば行列aと行列bの積であれば、np.dot(a,b)と書くだけです。

例えば上記のシンプルな行列\(a=\begin{bmatrix} 1 & 2 & 3 \end{bmatrix}\)と、その転置行列の積をPythonで計算する場合、非常に短いコードで済みます。

コード例            出力例

SimpleMatrix

このようにNumpyを使うと、転置行列の作成、行列の積が簡単にできます。

・csvファイルの読み込みについて

numpyのloadtxtを使うと簡単に読み込むことができます。

import numpy as np
a = np.loadtxt("sampledata.csv",delimiter=",")

numpyのloadtxtを使うと簡単に読み込むことができます。

csvファイルを読み込んで、データの統計処理を行うのも簡単です。

最大値
print(np.amax(a))
最小値
print(np.amin(a))
算術平均
print(np.mean(a))
中央値
print(np.median(a))
標準偏差
print(np.std(a))
分散
print(np.var(a))

Matplotlib基本コード

こちらのサイトが完璧です。

よく私がよく使う機能だけを以下まとめておきます。

import matplotlib.pyplot as plt
plt.plot(x,y,'r--', label='line1', linewidth=3.0, color="blue", linestyle="--", alpha=0.3) #linewidth:線の太さ、color:以下"色の変更"参照、linestyle:以下"線の種類変更"参照、alpha:透明度
plt.xlabel('time') #ラベル名の設定
plt.legend(loc='best',ncol=3, fancybox=True, shadow=True) #凡例の設定(詳細は"凡例の設定"参照)
plt.rc('grid', linestyle="--", color='gray')  # Gridのカスタマイズ
plt.grid(True) # Gridの表示
plt.show()
色の変更

・線の色の設定例

color = "blue"

色の種類はこちらにまとめられています。

基本的な色は以下テーブルのように、そのまま色の名称を記載するだけです。以下にdark、lightをつけることで色の深みを変更することも可能です。

bluegreenredyellowblack
凡例の設定

・凡例の設定例

plt.legend(loc='best',ncol=3, fancybox=True, shadow=True)

こちらこちらこちらが参考になります。

Gridのカスタマイズ

・Gridの設定例

plt.rc('grid', linestyle="--", color='gray')
plt.grid(True) # Gridの表示

こちらが参考になります。

・2dplot出力について

詳細はゆくゆくまとめていきたいと思いますが、matplotlibのplotはsubplotをまとめて一つの形式で出力できるため、非常に便利です。

例えば、NRTLで水-EthanolのXY、PXY、活量係数をまとめてプロットしたいという場合、以下のようにまとめて表示させることが可能です。

NRTLplot

以下のように各subplotの配置位置を2dgridで記載し、各suubplotの詳細を記載していきます。

subplots

・3dplot出力について

Qiitaのブログに詳細記載しましたのでご参照ください。

例として、以下のように理想気体と実在気体のPVT図を重ね書きする方法をまとめています。時間をみつけてもう少し詳細に書きたいとは思っていますが、、、

MatplotlibGraph

・ヒストグラムの出力について

基本コード

plt.hist(データ,bins=ビン数)

テキストファイルの読み込み

・物性基本情報を記載したテキストファイルの読み込みについて

例えば、以下のVapor pressureのデータ(データはperryのhandbook参照)をテキストファイルに保存しておき、各係数を必要に応じて読みだすことを考えます。

$$ln(P)=C1+C2/T+C3ln(T)+C4T^{C5}$$

C1:73.649、C2:-7,258.2、C3:−7.3037、C4:4.1653E−06、C5:2

ただし、圧力の単位はPa、適用温度範囲は、273.16~647.096[K]となっています。

上述のC1~C5パラメーターを各行毎に区切ってテキストに入力しておけば以下のような形で読み込むことができます。

VPreadCode

①f=open('VPwat.txt')で同階層にあるVPwat.txtを開き、

②lines2=f.readlines()でテキストファイルを1行ずつlines2に格納していきます。

※readlines()以外にも、read()、readline()があります。使い分けは以下のリンクが分かりやすいので、ご参照ください。

  • read()
  • ①f.close()でテキストファイルを閉じます。

    ④他の蒸気圧推算法でも対応できるように、VPの係数はMax7個まで格納できるようにしています。この際、Cに格納するデータが空だとエラーが出るため、空の部分には0を代入しています。

    ⑤C1~C7までデータを格納していきます。この際、リストからそのままデータを引っ張ってくると改行コードが邪魔であるのと、数値として認識されなかったため、.strip()で改行コードを除き、float()でそれぞれを数値化しています。

    このように、係数は分離して別置きしておけば、推算式が変わっても呼び出し元の関数を使い分けるだけで対応することができます。

    ・物性基本情報を記載したExcelファイルの読み込みについて

    Excelファイルの読み取り

    まずxlrdを読み込みます。

    import xlrd

    例えば、物性データを「PhysicalData.xlsx」として保存し、ここから必要な情報を読み取る場合、以下のように設定します。

    Excelデータ

    ExcelData

    Pythonコード

    xlrdread

    簡単に各行のコードをご紹介します。

    Excelファイルの読み込み

    book = xlrd.open_workbook('PhysicalData.xlsx')

    シートの読み込み(index(0)がsheet1)

    sh1 = book.sheet_by_index(0)

    セルの値の読み込み(例えば、sheet1のB2セルの値の読み取り)

    sh1.cell(2, 2).value

    Scipyの活用

    ・interpolate(補間)

     

    こちらのリンクに使い方がまとめられています。

      Cubic Spline

    CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)

    ・bc_type:境界条件の決め方

    natural:natural cubic spline(端点の2次微分が0)

    not-a-knot:両端点の1つ手前の点において3次微分が連続であること(Default設定)

    ・extrapolate:外挿範囲の設定

    periodic or None

    lagrange(x,y) Akima1DInterpolator(x,y)  

    こちらのリンクで紹介されている関数:\(sin(20x)+exp(5x/2)\)について各補間法を試してみます。

    補完

    ・optimize(leastsq関数)

    こちらのリンクを参考にさせていただきました。

    1次遅れプロセスのシステム同定

    1次遅れプロセスのシステム同定をleastsq関数で求めます。

    同定する1次遅れプロセスを、以下の伝達関数と仮定し、ステップ信号を与えた場合の時間応答データから、各パラメーターK、τ、Cを求めていきます。

    $$G(s)=\frac{K}{\tau s + C}$$

    leastsq関数は、leqstsq(最小二乗和誤差関数, パラメーターの初期値)の形で用います。

    まず、適当なパラメーターについて、伝達関数のステップ時間応答のデータをツールで出力します。

    以下、K=1、τ=5、C=1とした場合の出力結果です。

    1次遅れ時間応答出力テーブル

    出力結果を一旦エクセルにコピペし、乱数を与えて数値を少しズラした以下のデータをサンプルデータとして与えます。(csv形式で保存し、配列にデータを読み込みます)

    1次遅れ時間応答出力テーブル2

    初期値を、K=1、τ=1、C=1として、各パラメーターをleastsq関数で求めたプロットが以下になります。

    1次遅れ時間応答同定後

    データのバラつきをあまり強くしなかったので、サンプル点と殆ど重なって見えていますが、推算パラメーターは問題無さそうです。

    leastsq関数で求めた各パラメーターは、以下の通り。

    K=197.039744388、C=197.2462628、 τ= 983.033708335

    プロットからも明らかですが、本来の解:K=1、C=1、τ=5とほぼ同等の値が得られており、時定数における出力値が1-1/eとなっていることが確認できます

    NRTLのバイナリパラメーター回帰

    もう少し複雑な関数の最適化事例としてNRTLのパラメーター推算をこちらに整理しています。

    ・odeint,ode(常微分方程式のsolver)

    こちらのリンクを参考にさせていただきました。

    上記のリンクによると、odeinitを使うと、硬い問題かどうかによって、採用手法を自動で切り替えてくれるとのこと。

    odeを使うと、解法を自由に選択できるようなので、まずはodeを試してみます。

    主要な解法を以下に記載します。

    BDFの呼出し:ode(function).set_integrator('vode', method='bdf')

    adamsの呼出し:ode(function).set_integrator('vode', method='adams')

    lsodaの呼出し:ode(function).set_integrator('lsoda')

    ode45の呼出し:ode(function).set_integrator('dopri5')

    (functionの部分に関数定義します。)

    例題として、stiffな問題の方が差異が確認しやすいので、Van Der Pol方程式で各手法を比較してみます。

    比較前提

    Time step:0.1[sec]、積分時間:2000[sec]、μ:1000

    van der pol方程式は以下の通り。

    $$\frac{d^2x}{dt^2}-\mu \left(1-x^2\right)\frac{dx}{dt}+x=0$$

    Adamsは今回の条件では解くことがでなかったため、その他の計算手法の計算速度を比較してみると以下のようになります。

    ODE計算時間比較

    BDF:19.71[sec]、lsoda:0.14[sec]、ODE45:13.73

    続いてodeintを試してみます。

    odeint(function, X0, t)の形で定義します。(X0が初期値、tが時間定義になります)

    odeinit

    計算時間わずか0.0164[sec]。

    基本的にはまずodeintを使うという戦略は妥当そうです。

    なお、上記VDPの解をmatplotlibで作図すると以下のようになります(BDF)。

    VDPplot

    ・signal(Bode線図の作成、時間応答)

    こちらのリンクを参考にさせていただきました。

    伝達関数の定義(lti関数)

    例によって、以下の1次遅れプロセスについての定義

    $$G(s)=\frac{1}{5 s + 1}$$
    Gs = signal.lti([1], [5, 1])
    Magnitude、Phaseの計算(bode関数)
    freq, mag, phase = signal.bode(伝達関数)
    Matplotlibにおけるx軸対数プロットの作成(semilogx)
    plt.semilogx(freq, phase)

    Plot結果は以下の通り、正しい結果が得られています。

    Bode plot

    時間応答(Impulse、Step関数)
    time,y=signal.impulse(伝達関数)
    time,y=signal.step(伝達関数)

    Plot結果は以下の通り、正しい結果が得られています。

    Time response

    scikit-learnの活用

    ・datasetの生成(sklearn.datasets)

    回帰用データ(内装データ)

    ・ボストンの住宅価格

    load_boston

    参照)The Boston house-price data of Harrison, D. and Rubinfeld, D.L

    Number of Instances506
    Number of Attributes 13 numeric/categorical predictive
    Target Median Value (attribute 14)

    各パラメーターの詳細

    CRIMper capita crime rate by town
    ZNproportion of residential land zoned for lots over 25,000 sq.ft.
    INDUSproportion of non-retail business acres per town
    CHASCharles River dummy variable (= 1 if tract bounds river; 0 otherwise)
    NOXnitric oxides concentration (parts per 10 million)
    RMaverage number of rooms per dwelling
    AGEproportion of owner-occupied units built prior to 1940
    DISweighted distances to five Boston employment centres
    RADindex of accessibility to radial highways
    TAXfull-value property-tax rate per $10,000
    PTRATIOpupil-teacher ratio by town
    B1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
    LTSTAT% lower status of the population
    MEDVMedian value of owner-occupied homes in $1000's

    ・糖尿病患者データ

    load_diabetes

    ・linnerデータ

    参照)Tenenhaus, M. (1998). La regression PLS: theorie et pratique. Paris: Editions Technic.

    こちらのサイトが詳しいです。

    load_innerud
    Number of Instances20
    Number of Attributes 3
    回帰用データ生成関数

    ・make_regression

    make_regression(n_samples=100, n_features=100, n_informative=10, n_targets=1, bias=0.0, effective_rank=None, tail_strength=0.5, noise=0.0, shuffle=True, coef=False, random_state=None)

    生成例

    make_regression(n_samples=100, n_features=2, n_informative=1, n_targets=1, tail_strength=0.5, noise=50.0, random_state=None)
    makeregressiondata

    ・回帰

    linear_model

    Support Vector Regression

    こちらのサイトを参考に、無駄時間1次遅れプロセスの表現が可能か試してみました。

    Excelで1次遅れ&無駄時間プロセスの時間データを与えて、SVRで確認してみます。

    カーネル関数はRBFとPolynominalを試し、上記リンク同様、以下のパラメーターを採用しています。

    svr_rbf = SVR(kernel='rbf', C=1e3, gamma=0.1)
    svr_poly = SVR(kernel='poly', C=1e3, degree=2)

    C:罰則項、gamma : 係数、degree : 次数(default : 2)

    FOPDT

    雰囲気は出ていますが外挿範囲も含め、さすがに無理がありますね。

    ・もう少しトレースしやすそうな例としてScipyの活用でご紹介した1次遅れプロセスの表現をRBFでも行ってみます。

    FirstOrderLag

    tensorflowの活用

    ・非線形式のパラメーター

    例1)1次遅れプロセスの表現

    こちらのサイトを参考にさせていただきました。

    ・1次遅れプロセスの時間データを学習の入力値とし、1次遅れ時間、プロセスゲインを出力値として学習させます。

    どのくらいのデータが学習に必要か分からないので、1次遅れ時間を変更したケースを10ケース、ゲインを変更したケースを10ケース、それぞれについての時間データを20個用意し、学習ケースとして用意しました。

    データの生成は以下のようにエクセルで行いました。(for文で回しながら順番にcsv出力)

    traindata

    各csvデータを入力として書き込みます。

    例2)無駄時間1次遅れプロセスの表現

    おまけ

    ・絶対に勝てない石取りゲーム

    息子用に遊びで作ったコードです。意味は無いですが、対話型コードの参考として見てみてください。

    While,If文の用法も確認でき、意外とこうした遊び感覚で作るプログラムの方が勉強になったりします。

    ※1 Shellでの入力要求はinput()の形で書きますが、そのままでは文字列として認識されてしまいます。

    このため、数値として使用するためには、int(input())として一旦数値に変換する必要があります。

    逆にprintで出力するには、文字列に変換する必要があるため、str()を使用します。

    当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

    Page Top