2021-03-27

変分法のオイラー方程式の変形

Rudan, Massimo. Physics of Semiconductor Devices. 2015. Springer. にあった変形

オイラー方程式 $\displaystyle \frac{\mathrm{d}}{\mathrm{d}\xi} \frac{\partial g}{\partial \dot w} = \frac{ \partial g }{ \partial w }$ (1.4) から(1.5)の導出が分からなかった．

$\displaystyle \frac{\mathrm{d}}{\mathrm{d}\xi} \frac{\partial g}{\partial \dot w}$ の微分の順番を入れ替えて，（無条件で入れ替えられるかは不明） $\displaystyle \frac{\partial }{\partial \dot w} \frac{\mathrm{d} g}{\mathrm{d}\xi} = \frac{ \partial g }{ \partial w }$ (1.4)'

ここで $g(w,\dot w, \xi)$ に多変数関数の偏微分の連鎖律を適用して $\displaystyle \frac{\mathrm{d}}{\mathrm{d}\xi} g(w(\xi),\dot w(\xi), \xi) = \frac{\partial g}{\partial w} \frac{d w}{d \xi} +\frac{\partial g}{\partial \dot w} \frac{d \dot w}{d \xi} +\frac{\partial g}{\partial \xi} \frac{d \xi}{d \xi} = \frac{\partial g}{\partial w} {\dot w} +\frac{\partial g}{\partial \dot w} {\ddot w} +\frac{\partial g}{\partial \xi}\cdot 1$

$\displaystyle \frac{\partial }{\partial \dot w} \frac{\mathrm{d}g}{\mathrm{d}\xi} = \frac{\partial }{\partial \dot w} \left( \frac{\partial g}{\partial \dot w} {\ddot w} + \frac{\partial g}{\partial w} {\dot w} +\frac{\partial g}{\partial \xi} \right) = \frac{ {\partial}^2 g}{ {\partial \dot w}^2} {\ddot w} + \frac{ {\partial}^2 g}{\partial w \partial \dot w} {\dot w} +\frac{ {\partial}^2 g}{\partial \xi \partial \dot w}$ これを(1.4)'の左辺として

$\displaystyle \frac{ {\partial}^2 g}{ {\partial \dot w}^2} {\ddot w} + \frac{ {\partial}^2 g}{\partial w \partial \dot w} {\dot w} +\frac{ {\partial}^2 g}{\partial \xi \partial \dot w} = \frac{ \partial g }{ \partial w }$ (1.5)

2021-03-08

プログラムが終了したらメールを送信するpythonプログラム

プログラミング

windows 7や10で動作するはず．

使い方

"監視する実行ファイル名.exe"がすでに実行されているときに，以下のプログラムを開始する．

time.sleep(10)のとおり，10秒おきにプログラムが終了したか確認する．

プログラムが起動していなければ（終了していれば）メールを送信し，このpythonプログラム自体も勝手に終了する．

import smtplib
import subprocess
import time

while True:
    #　以下は https://a-zumi.net/python-windows-tasklist/ を参考にした
    cmd = 'tasklist'
    proc = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE)
    time.sleep(10)
    for line in proc.stdout:
        s=str(line)
        if "監視する実行ファイル名.exe" in s: 
            break
            pass
    else:
        break
        pass
    continue
    pass

print("finish")

# 以下は https://blog.ikappio.com/python-send-mail-connecting-gmail-using-smtplib/ を参考にした

# Import the email modules we'll need
from email.message import EmailMessage

# Create a text/plain message
msg = EmailMessage()
msg.set_content('メール本文')

msg['Subject'] = 'メールタイトル'
msg['From'] = '送信元のアドレス'
msg['To'] = '送信先のアドレス'

# Send the message via our own SMTP server.
s = smtplib.SMTP('smtp.gmail.com', 587)

s.starttls()
# 2つ目の引数に生成したアプリパスワードを指定する
s.login("送信アカウントのアドレス", '送信アカウントのアプリパスワード')

# 送信
s.send_message(msg)
s.quit()

2021-01-01

格子ボルツマン法の局所平衡分布関数を自動で展開するプログラム

プログラミング

格子ボルツマン法(Lattice Boltzmann Method, LBM)の局所平衡分布関数(local equilibrium distribution)を自動で展開するプログラム．

数値計算のループを実行する際に，展開された（ベクトル，行列を使わない）式を使うことで計算の高速化が期待できるかも？

格子ボルツマン法の局所平衡分布関数は下式で与えられる．下式は[Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]による．

この式に，決められた $w_i, c_{i\alpha}, c_s$ を代入して，整理された式を出力するプログラムを示す．

$\displaystyle{ f^{\mathrm{eq}}_{i} = w_{i} \rho \left( 1+ \frac{ c_{i\alpha} u_{\alpha} }{ c_{s}^{2} } + \frac{ u_{\alpha} u_{\beta} (c_{i\alpha} c_{i\beta} - c_{s}^{2} \delta_{\alpha \beta} ) }{ 2 c_{s}^{4} } \right) \tag{3.54} }$

記号の簡単な説明として，

$w_i$ は粒子の速度ベクトルの大きさ要素
$c_{i \alpha}$ は粒子の速度ベクトルの方向要素
$c_s$ は音速
$\delta_{\alpha \beta}$ はクロネッカーのデルタ．

添え字については，

$i$ は粒子の速度ベクトルの番号
$\alpha, \beta$ は速度方向の座標の番号（0,1,2がそれぞれx,y,zに対応）．

pythonプログラムのソースコードを以下に示す．

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import sympy as sym


# In[2]:


rho = sym.symbols("\\rho")
alpha, beta = sym.symbols("\\alpha \\beta")
cs = sym.sqrt(sym.Rational(1, 3))
u = sym.IndexedBase("u")
c = sym.IndexedBase("c")

w_D1Q3 = sym.Array([sym.Rational(2, 3), sym.Rational(1, 6), sym.Rational(1, 6),])
c_D1Q3 = sym.Array([[+0], [+1], [-1]])

w_D2Q9 = sym.Array(
    [sym.Rational(4, 9)] + [sym.Rational(1, 9)] * 4 + [sym.Rational(1, 36)] * 4
)
c_D2Q9 = sym.Array(
    [
        [+0, +0],  # 移動なし
        [+1, +0],  # +x方向
        [+0, +1],  # +y方向
        [-1, +0],  # -x方向
        [+0, -1],  # -y方向
        [+1, +1],  # +x, +y方向
        [-1, +1],  # -x, +y方向
        [-1, -1],  # -x, -y方向
        [+1, -1],  # +x, -y方向
    ]
)


w_D3Q15 = sym.Array(
    [sym.Rational(2, 9)] + [sym.Rational(1, 9)] * 6 + [sym.Rational(1, 72)] * 8
)
c_D3Q15 = sym.Array(
    [
        [+0, +0, +0],  # 移動なし
        [+1, +0, +0],  # +x方向
        [-1, +0, +0],  # -x方向
        [+0, +1, +0],  # +y方向
        [+0, -1, +0],  # -y方向
        [+0, +0, +1],  # +z方向
        [+0, +0, -1],  # -z方向
        [+1, +1, +1],  # +x, +y, +z方向
        [-1, -1, -1],  # 以下同様
        [+1, +1, -1],
        [-1, -1, +1],
        [+1, -1, +1],
        [-1, +1, -1],
        [-1, +1, +1],
        [+1, -1, -1],
    ]
)


# In[3]:


def f_eq(i, w, c):
    Q = c.shape[1]
    _term0 = 1
    _term1 = sym.Sum(c[i, alpha] * u[alpha] / (cs ** 2), (alpha, 0, Q - 1)).doit()
    _term2 = sym.Sum(
        sym.Sum(
            u[alpha]
            * u[beta]
            * (c[i, alpha] * c[i, beta] - (cs ** 2) * sym.KroneckerDelta(alpha, beta))
            / (2 * cs ** 4),
            (alpha, 0, Q - 1),
        ),
        (beta, 0, Q - 1),
    ).doit()
    return sym.simplify(rho * w[i] * (_term0 + _term1 + _term2))
    pass


# In[4]:


f_eq(5, w_D2Q9, c_D2Q9)

jupyter labで実行時の出力結果：

$\displaystyle{ \displaystyle \frac{\rho \left(3 {u}_{0}^{2} + 9 {u}_{0} {u}_{1} + 3 {u}_{0} + 3 {u}_{1}^{2} + 3 {u}_{1} + 1\right)}{36} }$

式の $\rho$ は密度， $u_0, u_1, u_2$ はそれぞれ，x, y, z方向の速度を表す．

この結果は2次元の9速度モデル(D2Q9モデル)の「5番目の式」： $f^\mathrm{eq}_5$ を表し， [Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]の式(3.65)に含まれる $f^\mathrm{eq}_5$ と同じ意味を表す．

f_eq(1, w_D1Q3, c_D1Q3)

を実行時の結果：

$\displaystyle{ \displaystyle \frac{\rho \left(3 {u}_{0}^{2} + 3 {u}_{0} + 1\right)}{6} }$

この結果は1次元の3速度モデル(D1Q3モデル)の「1番目の式」: $f^\mathrm{eq}_1$ を表し， [Timm Krüger, et. al., The Lattice Boltzmann Method, Springer, 2017.]の式(3.64)に含まれる $f^\mathrm{eq}_1$ と同じ意味を表す．

2020-12-20

JupyterLabでコード整形（フォーマッタはautopep8）を導入

プログラミング

下記のサイトのコードだと"black"とかいうフォーマッタを使った例になってるらしい（一般名詞で大変紛らわしい）

https://jupyterlab-code-formatter.readthedocs.io/en/latest/how-to-use.html

{
    "shortcuts": [
        {
            "command": "jupyterlab_code_formatter:black",
            "keys": [
                "Ctrl K",
                "Ctrl M"
            ],
            "selector": ".jp-Notebook.jp-mod-editMode"
        }
    ]
}

自分はフォーマッタとしてautopep8を使おうと思ったから， cmdで"conda install -c anaconda autopep8"を実行した．

よく分からないまま色々なページを参照しながら， conda installの他にも別のコマンドを実行したから，そういったことも必要かもしれない．

上記のコードの"black"は"autopep8"に変更し，ショートカットキーはPyCharmのコード整形に合わせて，"Ctrl Alt L"にした．

{
    "shortcuts": 
    [
        {
            "command": "jupyterlab_code_formatter:autopep8",
            "keys": [
                "Ctrl Alt L"
            ],
            "selector": ".jp-Notebook.jp-mod-editMode"
        },
        
    ] 
}

上記のコードを，下図の右側にある"User Preferences"の欄に入力した．下図はJupyterLab Version 2.2.9の実行画面． f:id:kazmus:20201220200028p:plain

その後Ctrl+Sで設定を保存し，一応JupyterLabを再起動した．

Codeセル内を"Mode: Edit"で選択した状態（セル内にテキストカーソルが点滅している状態）で Ctrl+Alt+Lを同時押しするとセル内コードの整形が実行された．

2020-12-15

scipyのsolve_ivpでローレンツ方程式を解く

数学プログラミング

以下のページを参考にした． https://org-technology.com/posts/ordinary-differential-equations.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def lorenz(t, state, p, r, b):
    x, y, z = state
    return [-p * x + p * y, -x * z + r * x - y, x * y - b * z]
    pass


p = 10
r = 28
b = 8 / 3

v0 = [0.1, 0.1, 0.1]
t_end = 100
t = np.arange(0, t_end, 0.01)

sol = solve_ivp(lorenz, [0, t_end], v0, args=(p, r, b), dense_output=True)

z = sol.sol(t)

plt.plot(t, z.T)
plt.show()

2020-12-01

フィボナッチ数列の兎の個体数増加をグラフで

プログラミング数学

フィボナッチ数列的な個体数の変化を，下記リンクの図のように表したかった． https://commons.wikimedia.org/wiki/File:FibonacciRabbit.svg

生成される図は下記のようになる． f:id:kazmus:20201201193004p:plain

上記の図は以下のプログラムで生成した．

import matplotlib.pyplot as plt
from matplotlib.pyplot import plot
from collections import deque


def plotter(x_pos, y_poss):
    """
    (x0, y0) = (x_pos, y_poss[0])の点と，
    (x0, y0) = (x_pos+1, y_poss[1])の点を結ぶ直線を引く
    y_possは2つ以上の要素を持つ
    """
    plot([x_pos, x_pos + 1], [y_poss[0], y_poss[1]], "o-")
    pass


qs = deque([(0, 0, 0)])  # qsの要素のタプルは，(1世代前のy座標, 自身のy座標, 自身が分裂するかのフラグ)
print("qs ", qs)
max_loop = 5  # あまり大きくするとプログラムの実行に時間がかかる
for i_loop in range(max_loop):
    ql = len(qs)  # 今の世代の個体数
    qs_new = deque([])  # 新しい世代
    cnt_made_e = 0  # 新しい世代を何個体作ったか，カウントしておく
    for i_ql in range(ql):  # 今の世代をすべて見ていく
        e = qs.popleft()
        if e[2] == 0:  # 「自身が分裂するかのフラグ」が0なら，分裂せずに自身は分裂できるようになる（＝「分裂のフラグ」を1にする）
            # 自身を，「分裂のフラグ」を1にして新しい世代に生成

            # 新しく追加する要素のタプルは，(eの（＝1世代前の）y座標, i_loopのループで生成された個体の履歴, 「分裂のフラグ」)
            qs_new.append((e[1], cnt_made_e, 1))
            plotter(i_loop, qs_new[-1])
            cnt_made_e += 1
            pass
        elif e[2] == 1:  # 「自身が分裂するかのフラグ」が1なら，自身と，分裂できない個体（＝「分裂のフラグ」が0の個体）を生成する
            # 自身を新しい世代に生成
            qs_new.append((e[1], cnt_made_e, 1))
            plotter(i_loop, qs_new[-1])
            cnt_made_e += 1

            # 自身から分裂した個体を生成
            qs_new.append((e[1], cnt_made_e, 0))
            plotter(i_loop, qs_new[-1])
            cnt_made_e += 1
            pass
        else:
            # 「分裂のフラグ」におかしい値があれば途中で終了
            print(qs)
            exit(-1)
            pass
        pass
    qs = qs_new
    print("qs ", list(qs))
    pass

plt.show()

2020-11-25

CNNでのBatch Normalizationの平均・分散

機械学習

畳み込みニューラルネットを使ってるときにバッチ正規化したら， BN層のパラメータが(ch数×4)になった．

これはchごとに正規化していて，各chに含まれる各要素（各画素や各時刻の電圧）をそれぞれ正規化しているわけではない．

BNの解説では，ミニバッチに含まれるデータの「対応する要素」を全体として，その全体の平均・分散を取る．

しかし，チャンネルごとに正規化するなら，どこを「全体」とするかがわからなかった．

CNNでのBNはchごとに正規化するというのは正しい．

また，BNの平均・分散を計算する際の「全体」は下記のページのソースコードが示すように，

http://d2l.ai/chapter_convolutional-modern/batch-norm.html#implementation-from-scratch

#When using a two-dimensional convolutional layer, calculate the
#mean and variance on the channel dimension (axis=1). Here we
#need to maintain the shape of `X`, so that the broadcasting
#operation can be carried out later
mean = X.mean(axis=(0, 2, 3), keepdims=True)
var = ((X - mean) ** 2).mean(axis=(0, 2, 3), keepdims=True)

引用部分の"mean = X.mean(axis=(0, 2, 3), "が示すように，ミニバッチ方向とあるchに含まれる全要素を「全体」として，平均・分散を計算する．

雑感等

音楽，数学，語学，その他に関するメモを記す．

変分法のオイラー方程式の変形

プログラムが終了したらメールを送信するpythonプログラム

使い方

格子ボルツマン法の局所平衡分布関数を自動で展開するプログラム

JupyterLabでコード整形（フォーマッタはautopep8）を導入

scipyのsolve_ivpでローレンツ方程式を解く

フィボナッチ数列の兎の個体数増加をグラフで

CNNでのBatch Normalizationの平均・分散