Python3で数理生態

プログラミングの基礎

配列，ファイル入出力，Control Structure
関数，値渡し，参照渡し（要注意！）
関数，クラス，メソッド，インスタンス、モジュール

[1] 配列、ファイル入出力、Control Structure

ここで解説する内容はtest_01.pyとしてまとめてあるのでダウンロード可能．

クラスの読込み

数値計算用のクラスとファイル入出力用のクラスをそれぞれ読み込む（無ければインストール）．

import numpy as np   #for numerical calculations

import pandas as pd  #read with pandas as dataframe

このとき "import hogehoge as hg" とかすることによって、長い名前のクラスを省略などして使える。

配列の使い方

数値計算用の配列（サイズ固定）を全要素初期値ゼロで用意する。たとえば、10行（row）2列(column)のものを用意するには、以下の通り．

v01 = np.zeros([10,2])  #note that the index starts from zero

初期値をすべて１で用意するなら、以下の通り。

v02 = np.ones([10,2])  #note that the index starts from unity

配列の各要素にアクセスするには添え字を使って、

v01[1,1]

みたいな感じ。ただし、C言語と同様に、要素はゼロから始まる。v01[0,0]は存在するが、v01[10,2]は存在しない。

v01[0,0]

v01[10,2]

ループ文（for loop）

for loopにはC言語のときと違って範囲の指定の仕方に柔軟性があるが、単純な場合を想定する。配列v1の一列目の値を入れなおすことをループにするなら、以下の感じ。重要なのはコロン（：）を打つのを忘れないこと．配列の定義と同様、rangeの右端の値10はループに含まれない．

#examle of for loop

for i in range(0,10):

    v01[i,0] = i

条件文（if）

条件文も超簡単．

#example of if statement

if v01[1,0] == 0:

    print('v01[1,0] is 0.')

elif v01[1,0] == 1:

    print('v01[1,0] is 1')

else:

    print('v01[1,0] is neither 0 or 1.')

データフレーム (DataFrame)

ファイル入出力の準備として配列をデータフレームに変換してみる．ちなみにPyCharmのGUIの右下の変数リストを確認すれば，v01は{ndarray}という型であることが分かる．pandasのDataFrame変換メソッドを用いて，以下のように変換する．

#convert ndarray into dataframe

df1 = pd.DataFrame(v01)

列に名前をつけるにはたとえば以下の感じ．

df1.columns = ['time', 'N']

データフレームの0行目から3行目まで（３行目は含まず）を表示するなら以下のようにする．

df1[0:3]

データフレームdf1の内容をcsvファイルに保存するのは簡単．

#write into csv

df1.to_csv('test01.csv')

ただし、余計な情報も書き込まれる（csvファイルのA列）．

逆にcsvファイルからの読込みは以下のようにおこなう．

#read csv into  dataframe

df2 = pd.read_csv('test01.csv', header=None)

ヘッダー無しで読み込んでいるので、df2の内容を変数リストから"View as DataFrame"で開くと、下のようになる．

ヘッダー情報を除いて読み込みなおすには、以下の通り．ちゃんとゼロ行めを指定する．

df2 = pd.read_csv('test01.csv', header=0)

まだ、先ほどcsvファイルを作ったときに混入した要らない情報（0列め）があるのでそれを除く．”：”を使うと全行（もしくは）にアクセスできる．"iloc"属性(attribute)によって必要な情報だけにアクセスできる．

df2.iloc[:,[1,2]]

この内容を新しいデータフレームにコピーする．dir()文で既に生成した変数等の一覧を確認後、delで不要になったデータフレームdf2を削除する．

#copy and delete

df3=df2.iloc[:,[1,2]]

dir()

del(df2)

dir()

[2] 関数，値渡し，参照渡し（要注意！）

各言語の関数には値渡しおよび参照渡しの独自のルールがある．Pythonの場合は以下の通り．より理論的な解説は，https://qiita.com/makotoo2/items/3b7b44baf16d487d46aa

ここのサンプルコードはtest_02.pyでダウンロード可能．

まずは、単純に二つのパラメータa, bを持ち、和を出力する関数を定義してみる．

import numpy as np   #for numerical calculations

def test_sum(a, b):

    return a + b

そこにm, nという二つの引数(argument) を代入すると3が出力される．

m = 1

n = 2

test_sum(m,n)

次に，二つのパラメータを持ち二つ目の値を一つ目に加える関数を定義してみる．

def test_add(a, b):

    a += b

    return a

これを実行すると当然3が出力される．

def test_add(a, b):

    a += b

    return a

test_add(m,n)

だがしかし，mの値は更新されない．したがって，この場合は値渡しである．

print(m)

しかし、数値配列v01を定義し，配列をパラメータとするような次の関数の場合は，引数の中身が更新される参照渡しである．ちなみにパラメータが一つの値の変数なのか配列なのか関数定義中のパラメータリストを一見しても分からないのが非常に気持ち悪い．C言語の仕様の方がいいでしょ．

v01 = np.zeros([10]) #define vector

v01.size  #check the size

#define a simple function

def test_func(vector_in, v_size):

    for i in range(0, v_size):

        vector_in[i] += 1.0

vec_size = v01.size

test_func(v01, vec_size)

v01

しかし、実はそんなに簡単じゃない（https://qiita.com/makotoo2/items/3b7b44baf16d487d46aa）．次の関数test_func_dは上と同じように機能する（参照元の値が変わる）が，ちょっと変えたtest_func_ddでは参照元の値が変わる．これは変数のidを確認すると分かる．

def test_func_d(vector_in):

    print("vector_in: {}".format(id(vector_in)))

    vector_in += 1.0

    print("vector_in: {}".format(id(vector_in)))

test_func_d(v01)

print(v01)

def test_func_dd(vector_in):

    print("vector_in: {}".format(id(vector_in)))

    vector_in = vector_in + 1.0

    print("vector_in: {}".format(id(vector_in)))

test_func_dd(v01)

print(v01)

演算子によって参照元の値が変わったり変わらなかったりするのはやっかいである．しかし計算速度の向上のためできるだけtest_funcの中で使っているようなforループは使わずに配列のまま計算したい．下のしてもダメ．つまり単純な代入演算子（＝）でも新たな変数が関数内部で生成され、参照元は変更されない．

def test_func_ddd(vector_in):

    print("vector_in: {}".format(id(vector_in)))

    vector_out = vector_in + 1.0

    print("vector_in: {}".format(id(vector_in)))

    vector_in = vector_out

    print("vector_in: {}".format(id(vector_in)))

print("v01: {}".format(id(v01)))

test_func_ddd(v01)

そんなバカな！！ここも参照（https://rcmdnk.com/blog/2015/07/08/computer-python/）．

ちなみにサイズをチェックする.sizeにどれだけ計算負荷があるかを確認するために，以下のtest_func2も定義して比較してみる．

def test_func2(vector_in):

    for i in range(0, vector_in.size):

        vector_in[i] += 1.0

timeというクラスをインポートし，v02というムダに大きな配列を準備，以下のスクリプトで比較したところ実行時間に大差が無いので，.sizeの負荷は無視できそう．

start =  time.time()

vec_size = v02.size

for i in range(0, 100):

    test_func(v02, vec_size)

elapsed_time = time.time() - start

print("elapsed_time:{0}".format(elapsed_time) + "[sec]")

start = time.time()

for i in range(0, 100):

    test_func2(v02)

elapsed_time = time.time() - start

print("elapsed_time:{0}".format(elapsed_time) + "[sec]")

[3] 関数，クラス，メソッド，インスタンス、モジュール

先にここのオイラー法のアルゴリズムを参照．オイラー法の部分をクラス、メソッド化してみる．numeric_odeという名前のクラスを定義しその中でパラメータのいくつかを自動的に初期化するためのinitメソッド(コンストラクタ）を定義し, オイラー法の一ステップ分の計算をしてくれる関数をeeという名前のメソッドで定義する．普通の関数との違いは一つ目の引数は必ずselfとすること．あと、アンダーバー二つ(__)で始まるものはprivate変数．

#definition of class for the numerical integration, including Explicit Euler and 4th-order RK

class numeric_ode:

    def __init__(self, h_interval, dim):

        self.h = h_interval

        self.dim = dim

    # definition of explicit Euler

    def ee(self, in_vec, out_vec, time, diff_vec):

        self.__k1 = np.zeros([self.dim])

        # calculate the next step vector

        diff_vec(in_vec, self.__k1, time)  # calculate k1

        for i in range(0, self.dim):

            out_vec[i] = in_vec[i] + self.h * self.__k1[i]

上の例で関数(メソッド）は関数を引数にすることが可能なことに注意．ここではeeの最後の引数diff_vecは関数である．

使い方は以下の感じ．

nv = np.zeros([1])

nv[0] = 0.1 #initial condition (initial population size)

time = 0.0 #initial condition (initial time, 0)

write_index = 1 #need for counting the numerical step

write_interval = 0.1 #with this interval, the result will be saved

#write initial condition

text=str(time)+','+str(nv[0])+'\n'

with open('result_ode3.csv', 'w') as f:

    f.write(text)

ここでクラスを実体化させたもの（＝インスタンス）をパラメータを与えて「生成」する．ここではnum_odeという名前のインスタンスを生成する．

num_ode = numeric_ode(dim=1, h_interval=h)

このインスタンスからオイラー法のメソッドをnum_ode.eeとして呼び出して使う．

for i in range(0, int(end_time/h)):

    num_ode.ee(nv, nv, time, diff)

    time += h   #update time

    ###recored the result every "write_interval" only

    if(write_index < int(write_interval/h)):

        write_index += 1

    else:

        text = str(time) + ',' + str(nv[0]) + '\n'

        with open('result_ode3.csv', 'a') as f:

            f.write(text)

        write_index = 1

#load result as a dataframe

table=pd.read_csv('result_ode3.csv', header=None)

plt.plot(table.iloc[:,0], table.iloc[:,1], color='blue', linestyle='solid')

クラス変数とインスタンス変数の違いがよくわからない．https://qiita.com/hiroykam/items/51a97b648b4f5eadf3b8
スコープの違いがよくわからない．https://qiita.com/0xfffffff7/items/6ef16e79fe9886acb3f2
モジュールとクラスの関係がわからない．

Google Sites

Report abuse