2020-06-06

多次元データ解析ライブラリ xarrayで複素数を利用する際の注意点　（つまずいた点)

はじめに

どうも，こんにちは．かわつあんです．

これまで私はpandasモジュールを使ってデータの管理を行ってきましたが，最近になってxarrayモジュールを使い始めました．個人的にはpandasよりもデータの管理がしやすく気に入っています．

xarray.pydata.org

特にpandasと比べて優れている点は，

標準で2つ以上の軸に対応していること
軸自体に値を代入できること

だと個人的に思っています．

xarrayに関する導入についてはxarray へのcontributorの方が日本語でレビューしてくださっていますので，以下のリンクより確認ください．

qiita.com

私はデータ解析では主にフーリエ変換などを利用するために複素数のデータが多いのですが, xarrayで複素数を扱おうとした際に少しつまずいたことがありました．

この記事でこのつまずきについて共有しようと思います．

ちなみにこの記事ではバージョン0.15.1を利用しています．

DataArrayを複素数配列でインスタンス化

xarrayのDataArrayを使って複素数配列でインスタンス化してみます．今回は簡単のため3x3の複素数データを扱います．

import numpy as np
import xarray as xr
complex_array1 = np.array([
    [1 + 5.j, 2 + 4.j,3 + 3.j],
    [2 + 4.j, 3 + 3.j,4 + 2.j],
    [3 + 3.j, 4 + 2.j,5 + 1.j],
])

da1 = xr.DataArray(
    complex_array1,
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
)
print(da1)

出力結果

<xarray.DataArray (t: 3, x: 3)>
array([[1.+5.j, 2.+4.j, 3.+3.j],
       [2.+4.j, 3.+3.j, 4.+2.j],
       [3.+3.j, 4.+2.j, 5.+1.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

このように，インスタンス化については直感的に利用できます． DataArrayは出力結果が綺麗でいいですね．特に軸のデータも一緒に見れるのが好きです．

DataArrayに複素数を代入する時につまずいた点)

インスタンス化するのはいいのですが，代入したりする時は注意が必要です．というのもインスタンス化する時に，DataArrayの中に入れるデータの型を明示しないといけないようです．

ダメな例

da2 = xr.DataArray(
    np.zeros((3,3)),
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
)
print(da2)

da2[0,0] = 1 + 1j
print("---after---")
print(da2)

出力結果（ダメな例）

<xarray.DataArray (t: 3, x: 3)>
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---after---
<xarray.DataArray (t: 3, x: 3)>
array([[1., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

要素が0の配列でインスタンス化したのちに[0,0]成分に1 + 1jを代入していますが，代入されたDataArrayを出力しても実数として反映されちゃっています．

これに気づかずに解析したら地獄をみますね...

良い例: dataの型を明示する

da3 = xr.DataArray(
    np.zeros((3,3)),
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
).astype("complex")       #<======  データの型を明示
print(da3)

da3[0,0] = 1 + 1j
print("---after---")
print(da3)

出力結果(良い例)

<xarray.DataArray (t: 3, x: 3)>
array([[0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---after---
<xarray.DataArray (t: 3, x: 3)>
array([[1.+1.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

このように代入した値もきちんと反映されています．インスタンス化してすぐのデータもきちんと複素数型になっています．

ちなみに，以下のように書いても問題なく複素数を扱えます．

da4 = xr.DataArray(
    np.zeros((3,3), dtype="complex"),  #<======  データの型を明示
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
)       
print(da4)

da4[0,0] = 1 + 1j
print("---after---")
print(da4)

出力結果

<xarray.DataArray (t: 3, x: 3)>
array([[0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---after---
<xarray.DataArray (t: 3, x: 3)>
array([[1.+1.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

np.zerosで要素が0の配列を作る時点でデータ型を複素数にしてしまえば，当然代入でも複素数は使えます．

要はインスタンス化の時点でDataArrayに入れるデータ型が決まってしまうということですね．

以上のつまずいた点はpython以外ではメモリの確保等のために型を明示するため常識だと思いますが，pythonで楽しようと思ってぽけーっと解析していたのでこのようにつまずきました．．．危なかった．．．

おまけ1（複素数の操作）

複素数といえば，複素共役や実部虚部を取り出したいときありますよね．これらの操作はnumpyのndarrayとほぼ同様にDataArrayで利用できます．

complex_array2 = np.array([
    [1 + 5.j, 2 + 0.j,3 + 3.j],
    [2 + 0.j, 3 + 3.j,4 + 0.j],
    [3 + 3.j, 4 + 0.j,5 + 1.j],
])                              

da5 = xr.DataArray(
    complex_array2,
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
)

print(da5)
print("---complex conjgate ( .conj ) ---")
print(da5.conj)
print("---complex conjgate ( .conj() ) ---")
print(da5.conj())
print("---imaginary part---")
print(da5.real)
print("---real part---")
print(da5.imag)

出力結果

<xarray.DataArray (t: 3, x: 3)>
array([[1.+5.j, 2.+0.j, 3.+3.j],
       [2.+0.j, 3.+3.j, 4.+0.j],
       [3.+3.j, 4.+0.j, 5.+1.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---complex conjgate ( .conj ) ---
<bound method _method_wrapper.<locals>.func of <xarray.DataArray (t: 3, x: 3)>
array([[1.+5.j, 2.+0.j, 3.+3.j],
       [2.+0.j, 3.+3.j, 4.+0.j],
       [3.+3.j, 4.+0.j, 5.+1.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1>
---complex conjgate ( .conj() ) ---
<xarray.DataArray (t: 3, x: 3)>
array([[1.-5.j, 2.-0.j, 3.-3.j],
       [2.-0.j, 3.-3.j, 4.-0.j],
       [3.-3.j, 4.-0.j, 5.-1.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---imaginary part---
<xarray.DataArray (t: 3, x: 3)>
array([[1., 2., 3.],
       [2., 3., 4.],
       [3., 4., 5.]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1
---real part---
<xarray.DataArray (t: 3, x: 3)>
array([[5., 0., 3.],
       [0., 3., 0.],
       [3., 0., 1.]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

このようにndarrayと同じようにDataArrayでメソッドや属性(attribute)で呼べば問題なく利用できます．

ただし，複素共役を扱う時は少し注意が必要です．ndarrayの場合ndarray.conjのように属性として複素共役を呼ぼうとしたらエラーがでますが，DataArrayはDataArray.conjを出力しようとするとエラーはでず，しかも複素共役ではなく元のデータが呼ばれてしまうようです．これはバグかもしれないので，githubで報告してみます． ==> すみません勘違いでした!

numpy.ndarary.conjもxarray.DataArray.conjもそれぞれ属性は存在せずメソッドのみ存在しているようです．出力だけではエラーはでず，numpy.ndarary.conjの場合はのようにメソッドですよという内容が出力される一方で，xarrayはメソッドですよって知らせてくれる内容に加えてDataArrayの中身も同時に出力してくれるみたいです．これですこしこんがらがってしまいました．ちなみにnumpy.ndarary.conjやxarray.DataArray.conjに何かしらの演算をしようとするとちゃんとエラーが出ますのでご安心ください．

おまけ2（複素数を含むデータの保存）

DataArrayでデータを保存する時はnetCDFという形式のファイルで保存するそうです．拡張子は.ncです．実数の場合は，以下の通り書きます

da6 = xr.DataArray(
    np.zeros((3,3)),
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
)
da6.to_netcdf("da6.nc")

複素数の場合，以下の通り少し工夫が必要なようです．

da7 = xr.DataArray(
    np.zeros((3,3)),
    dims = ["t", "x"],
    coords = {
        "t":[0, 1, 2],
        "x":[-1, 0, 1],}
).astype("complex")

# da7.to_netcdf("da7.nc")
# >>ValueError: NetCDF 3 does not support type complex128

da7.to_netcdf("da7_1.nc", engine="h5netcdf",invalid_netcdf=True)

netCDFにもバージョンがあって，対応するバージョンを指定しないといけないようですね．これにはpythonのモジュール"h5netcdf"が必要なのでpipなどでインストールしてください．詳しくは公式ドキュメントを参照してください．

xarray.pydata.org

.ncデータを読み込む時はopen_dataarrayを使います．保存した2つのデータを読み込んでみます．

da6_read = xr.open_dataarray("da6.nc")
da7_read = xr.open_dataarray("da7.nc", engine="h5netcdf")
print(da6_read)
print("------")
print(da7_read)

出力結果

<xarray.DataArray (t: 3, x: 3)>
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
Coordinates:
  * t        (t) int32 0 1 2
  * x        (x) int32 -1 0 1
------
<xarray.DataArray (t: 3, x: 3)>
array([[0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j]])
Coordinates:
  * t        (t) int64 0 1 2
  * x        (x) int64 -1 0 1

読み込みの際にも，オプションでengine="h5netcdf"を指定しないといけないのが注意点ですね．

さいごに

xarrayは日本語の情報は比較的少ないですが，すごく便利ですのでどんどん広まって欲しいと思っています．英語が読める方は公式ドキュメントも是非参考にしてくださいね．

xarray.pydata.org

ではまた(✿╹◡╹)

2020-01-20

動的モード分解(Dynamic Mode Decomposition)を実装してみた

python DMD データ解析

はじめに

こんにちは．かわつあんです．最近，私の研究の一つで動的モード分解(Dynamic Mode Decomposition, 以下DMD)というものを使っています．DMDは元来，流体解析などの高次元数をもつシミュレーションデータの次元削減を目的として発展してきた，データ駆動の解析の一種です．最近現れた解析で，理論的な研究がすすみ，私の流体解析分野や非線形力学系，そして私の研究分野のプラズマ物理におけるシミュレーション・実験への応用が徐々に行われてきています．

DMDの利点は教師となるような関数をひつようとしないデータ駆動的な計算で，主要モードとその周波数および減衰率まで計算できる点です．(フーリエ解析はsin波を，ウェーブレット解析ではマザーウェーブレットを教師としていますよね)．また，フーリエ解析と異なり時間分解能と周波数の最小幅が対応していません．一方で，DMDでは変調された構造や周波数が時間変化するような系ではその情報を取り出せないというデメリットもあります．

今回の記事では，DMDとそのアルゴリズムの説明とDMDの実装をいます．なお，コーディングは以下のサイトを参考にしています． www.pyrunner.com

私の理解に誤りがありましたら，是非コメントでご指摘等してくださると幸いです．

動的モード分解(DMD)

まずは普通のDMDを実装します． DMDの歴史や思想などは極力省略して，アルゴリズムの説明に重きを置いていますことをご留意ください．

各時間 $t_k, ( k=0,1, \cdots, m )$ 毎に得られた長さ $n$ の列ベクトルを $\bf x_k$ とします．これを普通に並べると

$\mathbf{X'} = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{x}_{0} & \mathbf{x}_{1} & \cdots & \mathbf{x}_{m} \\ \vert & \vert & & \vert \end{bmatrix}$

となると思います. DMDでは，はじめに次のように並べた２つ行列 $\bf X, Y$ を考えます．

$\begin{align} \mathbf{X} = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{x}_{0} & \mathbf{x}_{1} & \cdots & \mathbf{x}_{m-1} \\ \vert & \vert & & \vert \end{bmatrix}, \quad & \mathbf{Y} = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{x}_{1} & \mathbf{x}_{2} & \cdots & \mathbf{x}_{m} \\ \vert & \vert & & \vert \end{bmatrix} \end{align}$

これはすなわち， $\bf X$ は最後の時間のデータ以外を並べた行列であり， $\bf Y$ は最初の時間のデータ以外を並べた行列になっており， $\bf X$ と $\bf Y$ は時間的に１ステップだけズレています．実験やシミュレーションで得られるデータは非線形的で動的なものがほとんどですが，DMDでは動的なデータについて線形であると近似します（詳くは [J. H. Tu (2014)]などを参考）．

すなわち

$\mathbf{Y = AX}$

と書けると考えます． $A$ はある行列であり線形作用素です． $A$ を求めるためには

$\mathrm{min} \bf || AX - Y||_F$

となる $A$ を求めればよいはず．ここで， $\bf ||P||_F$ はフロベニウスノルム( $\bf P$ はなにかしらの行列)です．一般に $\bf X, Y$ は正方行列ではないので，一般的な行列の逆行列として定義される擬似逆行列を考えます．擬似逆行列 $\bf X^+$ とすると

$\bf A \approx YX^+$

として作用素 $\bf A$ を定められます．DMDでは作用素 $\bf A$ に対する固有値方程式を解き，得られた固有ベクトルからDMDモードを計算し，固有値からDMDモードの周波数および減衰率を求めます．

一般的なDMDでは擬似逆行列 $\bf X^+$ を求めるために特異値分解を行います．特異値分解は以下のように書けます．

$\bf X = U \Sigma V^H$

ここで， $\bf X$ は $n \times m$ 行列でしたが， $\bf U$ は $n \times n$ のユニタリー行列， $\bf \Sigma$ は $n \times m$ の対角行列， $\bf V$ は $m \times m$ のユニタリー行列です．また， $\bf V^{H}$ は $\bf V$ のエルミート共役です．エルミート共役とはある行列を転置して複素共役をとったもので，実数行列における転置行列に対応しています．ユニタリー行列はある行列とエルミート共役の積が単位行列になる場合，その行列はユニタリー行列であるといい（つまり $\bf U^{H} U = U U^{H} = E$ ， $\bf E$ は単位行列），実数行列の直行行列に対応しています．

以上のように特異値分解があるすると擬似逆行列 $\bf X^+$ は

$\bf X^+ = V \Sigma^{-1} U^H$

とかけるので， $A$ の射影行列 $\tilde{A}$ を

$\bf \tilde{A} = U^H A U = U^H YX^+U = U^HYV\Sigma^{-1}$

として求め， $\bf \tilde{A}$ についての固有値方程式

$\bf \tilde{A} w = \Lambda w$

を解きます．得られた $\bf \tilde{A}$ の固有値は $\bf A$ の固有値でもあるため，DMD アルゴリズムの目的である $\bf A$ の固有値方程式が解けたことになります．

standard-DMDの場合，固有ベクトルは以下の通り計算すると，

$\begin{align} \mathbf{\tilde{A}w} &=\mathbf{w\Lambda} \\ \mathbf{U^HAUw} &=\mathbf{w\Lambda} \\ \mathbf{A(Uw)} &=\mathbf{(Uw)\Lambda} \\ \end{align}$

となるため，固有ベクトルは $\bf Uw$ となります．

一方で，現在一般的に使われているexact-DMDの場合， $A$ に対応する固有値ベクトルは，以下のように書き直されます．

$\bf \Phi = YV\Sigma^{-1}w$

確認のため， $\bf A \Phi$ を計算すると，

$\begin{align} \mathbf{A \Phi} &=\mathbf{A Y V\Sigma^{-1}w } \\ &=\mathbf{(YV\Sigma^{-1}U^H )Y V\Sigma^{-1}w} \\ &=\mathbf{YV\Sigma^{-1}\tilde{A} w}\\ &=\mathbf{YV\Sigma^{-1}w\Lambda}\\ &=\mathbf{\Phi\Lambda}\\ \end{align}$

となり，確かに $\bf \Phi$ が $\bf A$ の固有ベクトルであることができました．以上がDMDの基本的なアルゴリズムですが，基本的なアルゴリズムにもいろいろあるそうなので，参考まで．

DMDを実装

では，先ほど説明したDMDを実装してみましょう．まずテストデータを作成します．

import numpy as np
import scipy.linalg                          #scipyの線形数学モジュール．　略記は非推奨らしい
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid

#時空間を定義
x = np.linspace(-10,10,60)
t = np.linspace(0,20,80)
Xm, Tm = np.meshgrid(x, t)

#テストデータ生成
D1 = 5 * (1/np.cosh(Xm/2)) * np.tanh(Xm/2) * np.exp((0.8j)*Tm) #主要モード
D2 = 0.2 * np.sin(2 * Xm) * np.exp(2j * Tm)                    #サブモード
D3 = 0.1 * np.random.randn(*Xm.shape)                          #ノイズ
D = (D1 + D2 + D3)
fig = plt.figure(figsize = (12,3))

#データをプロット
#Primary mode
sp1 = plt.subplot(1,4,1)
sp1.contourf(x, t, D1, cmap="jet", levels = np.linspace(-3,3,50))
sp1.set_aspect(1)
sp1.set_xlabel("x")
sp1.set_ylabel("t")
sp1.set_title("primary mode")
sp1.set_yticks([0,5,10,15,20])
#Secondary mode
sp2 = plt.subplot(1,4,2)
sp2.contourf(x, t, D2, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp2.set_aspect(1)
sp2.set_xlabel("x")
sp2.set_ylabel("t")
sp2.set_title("secondary mode")
sp2.set_yticks([0,5,10,15,20])
#noise
sp3 = plt.subplot(1,4,3)
sp3.contourf(x, t, D3, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp3.set_aspect(1)
sp3.set_xlabel("x")
sp3.set_ylabel("t")
sp3.set_title("noise")
sp3.set_yticks([0,5,10,15,20])
#created data
sp4 = plt.subplot(1,4,4)
cont=sp4.contourf(x, t, D, 50, cmap="jet", levels = np.linspace(-3,3,50))
sp4.set_aspect(1)
sp4.set_xlabel("x")
sp4.set_ylabel("t")
sp4.set_title("created data")
sp4.set_yticks([0,5,10,15,20])
plt.tight_layout()

実行結果

f:id:kawats1:20200120011449p:plain — テストデータ

このテストデータに対してDMDモードを計算してみます．先ほどのプログラムに追記の形で以下のように実装します．

def dmd(data, truncation = None):
    #step1　データから行列を生成
    X = data[: ,1:]
    Y = data[:,:-1]

    #step2　特異値分解
    U, Sigma, V_H = scipy.linalg.svd(X, False) #scipy.linalg.svdでは，特異値は特異値のリストとして出力される
    
    #step3　射影行列計算　（ランクの切り捨ても行う．）特異値のランクの切り捨ては引数truncationで指定
    r = len(Sigma) if truncation is None else truncation
    
    U = U[:,:r]
    V = V_H.conj().T[:,:r]
    U_H = U.conj().T
    Sigma = np.diag(Sigma)[:r,:r] #特異値のリストを対角行列に．
    Sigma_inv = scipy.linalg.inv(Sigma)
    #step4 A_tildeの固有値問題を解く．
    A_tilde = np.dot(np.dot(np.dot(U_H, Y),V), Sigma_inv)
    Lambda, W = scipy.linalg.eig(A_tilde) 
    #step5 固有値に対応するAに固有ベクトルを求める
    Phi = np.dot(np.dot(np.dot(Y,V), Sigma_inv), W )
    return Lambda, Phi

# 
# 空間パターン
eigen_val, modes = dmd(D,truncation = 3)
fig1 = plt.figure()
plt.title("Spatial Pattern")
plt.plot(x, modes)
plt.legend(["mode1","mode2","mode3"])
plt.xlabel("x")
fig1.savefig("SP")

# 時空間パターン
omega=np.log(eigen_val)/0.1              # 複素角振動数
fig2 = plt.figure(figsize = (6,3))
plt.suptitle("Spatialtemporal Pattern")
sp5=plt.subplot(1,2,1)
sp5.contourf(x, t, np.outer(W_[:,0],np.exp(omega[0] * t)).T, 50, cmap="jet" ) 
sp5.set_aspect(1)
sp5.set_xlabel("x")
sp5.set_ylabel("t")
sp5.set_yticks([0,5,10,15,20])
sp6=plt.subplot(1,2,2)
sp6.contourf(x, t, np.outer(W_[:,1]*(W_[:,1]),np.exp(omega[1] * t)).T,50, cmap="jet" ) 
sp6.set_aspect(1)
sp6.set_xlabel("x")
sp6.set_ylabel("t")
sp6.set_yticks([0,5,10,15,20])
plt.tight_layout()

実行結果

f:id:kawats1:20200120005602p:plain — DMDモードの空間パターン

f:id:kawats1:20200120005635p:plain — DMDモードの時空間パターン

無事，DMDを使って２つのモードの抽出に成功しました．（ノイズについてはあまり関係ないので，今回は省略しています．）なお，固有値方程式を解いて得られた固有値 $\lambda_j$ から，複素角振動数 $\omega_j$ を次の通り計算しています．

$\omega_j = \mathrm{log}(\lambda_j)/ \Delta t$

もちろん実周波数 $f_j$ ，成長率 $\nu_j$ はそれぞれ以下のように計算できます．

$f_j = \mathfrak{Im}(\omega_j)/2\pi,\quad \nu_j = \mathfrak{Re}(\omega_j)/2\pi$

ここで注意が一つあります．実は今回実装したのはDMD解析の中でも空間パターンの基底関数やその時間発展を求めただけです．この他にもDMDモードの振幅を求めたり，DMDモードの重ね合わせで元のモードを再現できるかどうかのvalidationなども必要ですが，今回は割愛しています．また，詳しい理論などにも触れておらず申し訳ありません．

そこらへんの内容は需要があればまとめますが，基本的には原著論文などを参考にしてください

＝＞参考

以上のコードはGistでまとめていますので参考にどうぞ．(すこしMarkdownが文字化け？してますが．．．)

gist58c04e83be338c5786317c114fdde62e

さいごに

今回はDMDのアルゴリズムを説明してプログラムを実装しました．実は，この記事には続きがありまして，本当はスパースDMDを勉強しながら実装しようと思っていました．しかし，その前提であるDMDの説明が長くなりそうでしたので，今回はDMDの説明にとどめています．

次回こそ，スパースDMDについても記事を書く予定なので，よろしくお願いいたします．

質問，コメント，指摘，お待ちしています．

ではまた(✿╹◡╹)

2020-01-15

Sympyでギリシャ文字を定義する．(メモ)

はじめに

こんにちは．かわつあんです．最近pythonで方程式などの文字式を計算するモジュールsympyを使っていろいろ遊び始めました．文字式を扱うとなるとアルファベットに加えてギリシャ文字を定義することが多いのではないのでしょうか．そこでこの記事では，ギリシャ文字をsympyで使う方法をまとめました．

直接入力によるギリシャ文字の定義

タイトルの通り，直接ギリシャ文字を入力して定義します．入力方法は例えばギリシャ文字(ミュー) $\mu$ の場合，日本語入力できる状態で"myu-"と入力すれば，変換候補に出てきます．

これをそのまま定義してしまいます．sympyでは文字の定義はSymbolsクラスを使うと以下の通りにコードがかけます．

import sympy as sm
sm.init_pringing()

mu = sm.Symbol("μ")
mu

実行結果

$\mu$

となります．

ここで注意ですが私の場合，sympyを使う時はJupyter Lab（Notebook）を使っています．これを使えば，コードを書く際に一行sm.init_printing()を加えると出力の時にLatex(後述)のようにきれいに表示できるからです．

また，変数をセルの最後に書いておけば，自動的に出力されます．（普通のpythonでは最後の行はprint(mu)としないと出力されない）以下のスクリーンショットのようにきれいな出力になるため，おすすめです．

f:id:kawats1:20200115110121p:plain

Latex記法を使ったギリシャ文字の定義

次にLatex記法をsympyで使う方法についてまとめます．sympyではLatex記法でも文字を定義できます．

Latexは大学生の理工系を卒業している人の半数以上は卒業論文等で使ったことがあると思います．簡単にいうと数式をきれいに書くための道具です．

例えば，Latex記法で $\mu$ を書きたい場合は,\muのように入力すれば，きれいなギリシャ文字およびそれらを組み合わせた数式などが記入できます．この記事では，ギリシャ文字のLatex記法の知識を前提としています．わからない人はこちらなどを参考にしてください＝＞参考

それでは早速，Latexの記法を使った記法で例として $\mu$ および $\theta$ について定義してみます．

import sympy as sm
sm.init_printing()

mu    = sm.Symbol("\mu")
theta = sm.Symbol("\theta")
mu, theta

出力結果 $( \mu, \quad heta )$

Jupyter Labのスクリーンショットは以下の通り． f:id:kawats1:20200115111544p:plain

ダメでしたね．\thetaの部分が正しく定義されていません. というのも，文字列の中の"\t"はpythonではタブ文字として定義されているからです．

これの解決策として一番簡単なのは，文字列"hoge"の部分をr"hoge"のようにかえれば良いです．

import sympy as sm
sm.init_printing()

mu    = sm.Symbol("\mu")
theta = sm.Symbol(r"\theta")# "\theta"から変更
mu, theta

実行結果

$( \mu, \quad \theta )$

スクリーンショット

f:id:kawats1:20200115112217p:plain

なぜ正しく定義できたのかというと，""の前にrを書くとraw stringsすなわち，そのままの文字列として認識し，予約された文字もただの文字として認識してくれるようになるそうです．詳しくは以下を参考にしてください． blog.pyq.jp

ちなみに，お気づきの方もいるとは思いますが,直接入力にたいして，Latex記法を使うと斜体になってギリシャ文字が定義されるようです．以下に，Latex記法で斜体を普通の形に直そうとしたものを比べてみましたが．．．

スクリーンショット f:id:kawats1:20200115113938p:plain このように，Latex記法では斜体を直せないようです．また，直接入力とLatex記法での出力結果はやはり違うようです．

なんでだろう．．．

とりあえずは，特に気にならなけらばLatex記法を普通に使うほうがいい気はします．．．

さいごに

今回はsympyでギリシャ文字を定義する方法についてまとめてみました．私は，Latex記法でギリシャ文字を定義するのをお勧めします．若干謎も残りましたが，まあいいでしょう．．．（よくない）

ではまた(✿╹◡╹)

かわつあん奮闘記

多次元データ解析ライブラリ xarrayで複素数を利用する際の注意点　（つまずいた点)

目次

はじめに

DataArrayを複素数配列でインスタンス化

DataArrayに複素数を代入する時につまずいた点)

ダメな例

良い例: dataの型を明示する

おまけ1（複素数の操作）

おまけ2（複素数を含むデータの保存）

さいごに

動的モード分解(Dynamic Mode Decomposition)を実装してみた

目次

はじめに

動的モード分解(DMD)

DMDを実装

さいごに

Sympyでギリシャ文字を定義する．(メモ)

目次

はじめに

直接入力によるギリシャ文字の定義

Latex記法を使ったギリシャ文字の定義

さいごに