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記法でギリシャ文字を定義するのをお勧めします．若干謎も残りましたが，まあいいでしょう．．．（よくない）

ではまた(✿╹◡╹)

2020-01-04

pythonistaがC++を使って配列を関数で扱う

C++ ctypes python

はじめに

どうも，かわつあんです．私のブログでは，これまでに"ctypesでC++の関数をpythonから呼び出したい"というモチベーションのもとpythonista目線で記事を投稿してきました．今回は，C++で行列を関数で扱ってみます．これができればやっと本題に進めます．内容はC++を知ってる方には簡単すぎると思いますので，御注意の程よろしくお願いします．

配列を引数にして関数内で処理をしてスカラーあるいは配列を返すようなことはデータ解析では日常茶飯事です． pythonだったら何も考えずにできますがCやC++では参照やポインタの概念を理解して注意して扱わないといけません．ちなみにvectorクラスを使えば，比較的pythonに近く書けるので今回はそちらを採用しています．（実は最近になって，pythonの関数も参照渡しだと知りました．pythonでも参照渡しであることに注意が必要なようですね．）

この記事では $3\times 3$ 行列の行列式を求める関数，および逆行列を求める関数をC++及びpythonで自作してこれらの比較をしてみます．なお，参照やポインタの概念の詳細はいろんな方がまとめていらっしゃいますのでしますのでご了承ください．

引数が行列，戻り値がスカラー（数値）の関数

例題として， $A = \begin{bmatrix} a_{11} ,a_{12} ,a_{13} \\ a_{21}, a_{22}, a_{23} \\ a_{31}, a_{32}, a_{33} \end{bmatrix}$ の行列式 $\mathrm{det} A$ を計算してみます． $3\times3$ 行列の行列式は数式的にはサラスの公式より，以下の通りにかけます．

$\mathrm{det}A= a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{32}a_{21} - a_{13}a_{22}a_{31} - a_{11}a_{23}a_{32} -a_{12}a_{21}a_{33}$

例えば， $A = \begin{bmatrix} 1,2,3 \\ 3,5,3\\3,8,3 \end{bmatrix}$ の行列式は $\mathrm{det}A = 18$ となります．

これをpythonで実装してみます．

import numpy as np

def det(A):
    det_A = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,1] + A[1,0]*A[0,2]*A[2,1] \
          - A[0,2]*A[1,1]*A[2,0] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2]
    return det_A 

a = np.array([[1,3,3],
              [2,5,8],
              [3,3,3]])
print("det A = {}".format(det(a)) )

出力結果

det A = 18

このように比較的に直感的にコーディングできますね．

一方で，C++ではvectorクラスを使うと以下のようにかけます．

#include <iostream>
#include <vector>
using namespace std;

int det (vector<vector<int>>A){
    int det_A;
    det_A = A[0][0]*A[1][1]*A[2][2] + A[0][1]*A[1][2]*A[2][1] + A[1][0]*A[0][2]*A[2][1]
          - A[0][2]*A[1][1]*A[2][0] -A[0][0]*A[1][2]*A[2][1] -A[0][1]*A[1][0]*A[2][2];
    return det_A ;
};

int main(){
    vector<vector<int>> a {{1,3,3},{2,5,8},{3,3,3}};
    cout << "det A = " << det(a) << endl;
}

出力結果

det A = 18

pythonと似たようにコーディングできましたね．ここまではC++初心者でも比較的理解しやすいと思います．しかしC++やCで難しいのは次のトピックで使う参照渡しですので，まだ安心しないでください．

引数が行列，戻り値が行列？の関数

さて次のトピックについては"戻り値が行列？の関数"というように？マークで強調してみましたが．．．

先に言っておくとC++(及びC)において戻り値を配列にすることはできません.

pythonばっかりやってきた僕は最初にこれみてびっくりしました．だったら配列の関数実装できないじゃん．．．て思っていましたが... ちゃんと方法はあります．それが参照渡しというやり方です．

参照渡しとは，ざっくり言うと関数内で計算した配列を記憶しているメモリの場所を引数として関数に渡すことです．関数内でメモリに格納された値を更新することで，関数処理が終えたあともメモリに計算後の値が残ります．fortranでいうサブルーチンのように使えばいいということですね．まあ，初学者の方は辛いと思いますが詳しくは他の文献などをご覧ください．

話はそれましたが，配列を計算して配列を返す関数を実装していきます．今回は簡単のため $3\times3$ 行列の足し算を行ってみます．行列の足し算は簡単で各成分ごとに足し算するだけです．具体的には，以下の通りです．

$A = \begin{bmatrix} 1,1,1 \\ 1,1,1, \\ 1,1,1 \end{bmatrix}, \quad B = \begin{bmatrix} 2,2,2 \\ 2,2,2, \\ 2,2,2 \end{bmatrix}, \quad C = A+B = \begin{bmatrix} 3,3,3 \\ 3,3,3, \\ 3,3,3 \end{bmatrix}$

まずは，これをpythonで実装していきます．numpyの配列なら何も考えず足すだけでは？というツッコミはなしでお願いします．あくまで今回は愚直に，各成分同士の足し算を成分に持つ配列を出力する関数を使います．

実装したコードは以下の通りです．

import numpy as np

def sum_matrix(A,B):
    C = np.zeros((3,3))
    for i in range(3):
        for j in range(3):
            C[i,j] = a[i,j] + b[i,j]
    return C

a = np.array([[1,1,1],
              [1,1,1],
              [1,1,1]])
b = np.array([[2,2,2],
              [2,2,2],
              [2,2,2]])
c = sum_matrix(a,b)
print("C = A+B =\n{}".format(c))

実行結果

C = A+B =
[[3. 3. 3.]
 [3. 3. 3.]
 [3. 3. 3.]]

まあ，愚直にかけばこんなもんでしょう．

C++だとどういう風にかけるのでしょうか？参照渡しを使って書くと次の通りです．

#include <iostream>
#include <vector>
using namespace std;

void sum_matrix (vector<vector<int>>A, vector<vector<int>>B, vector<vector<int>>& C){ //vector<vector<int>> C ではないことに注意
    for (int i=0; i<3;i++){
        for (int j=0; j<3;j++){
            C[i][j] = A[i][j] + B[i][j];
        }
    }
};

int main(){
    vector<vector<int>> a {{1,1,1},{1,1,1},{1,1,1}};
    vector<vector<int>> b {{2,2,2},{2,2,2},{2,2,2}};
    vector<vector<int>> c {{0,0,0},{0,0,0},{0,0,0}};

    sum_matrix(a, b, c);
    cout << "C = A+B ="<< endl;
    for (int i=0; i<3;i++){
        for (int j=0; j<3;j++){
            cout << c[i][j] << ", ";
        }
        cout << endl;
    }
}

出力結果

C = A+B =
3, 3, 3, 
3, 3, 3, 
3, 3, 3,

この場合，計算結果をいれる配列を初めに作って，その配列の参照を引数に入れます．配列Cの参照を&Cとして引数に入れ，関数内で配列Cのメモリの値を更新しています．ちなみに，もし

vector<vector<int>>&C

の部分を

vector<vector<int>>C

としてしまうと，値は更新されず出力結果が以下の通りになってしまいます．

C = A+B =
0, 0, 0, 
0, 0, 0, 
0, 0, 0,

これは，配列の参照ではなく値渡しと言う別の操作になってしまうためです．なお，ポインタ渡しでもコードを少し変えれば関数で配列計算等ができますが，ここでは省略します．

このように，参照渡し（あるいはポインタ渡し）の話は難しいので詳しくは他の文献等をみてください． Qiitaにいい感じの記事がまとまっていました．参考にどうぞ．

qiita.com

さいごに

今回は，C++言語で配列を関数で扱うために，pythonistaとしてpythonのコードと比較しながら実装しました．参照渡しという概念を意識しないといけないのでc++のコードは少し難しいですね．私のブログをみて少しでもpythonistaの方が参照渡しの雰囲気を味わえたのならそれだけで幸いです．

次回はctypesを使ってc++の関数を呼び出す内容を記事にしようと思います．

追伸いろいろ試行錯誤したのですが，現状ctypesではvectorの扱いが難しい，numpy配列のメモリの扱いがおかしい，など問題がある様なので，この後の記事はかけそうにないです．申し訳ない．．．

ではまた(✿╹◡╹)

かわつあん奮闘記

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

目次

はじめに

動的モード分解(DMD)

DMDを実装

さいごに

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

目次

はじめに

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

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

さいごに

pythonistaがC++を使って配列を関数で扱う

目次

はじめに

引数が行列，戻り値がスカラー（数値）の関数

引数が行列，戻り値が行列？の関数

先に言っておくとC++(及びC)において戻り値を配列にすることはできません.

さいごに