GAN備忘録
GAN概要
潜在変数をGeneratorの入力として受け取り、画像データを生成する．その画像データをDiscriminatorが入力として受け取り本物であるか偽物であるかの真偽値を確率として連続な値で返す．
プログラムの大枠
- inputの定義
- discriminatorモデルの定義
- generatorモデルの定義
Generatorモデル
Generatorの役割としては一様乱数から画像を表現するための配列を生成することである． そのため、はじめのinputの次元数は潜在変数の次元数で出力の次元数は画像をの解像度、チャンネル数を示した次元数となる．
def build_generator(self): noise_shape = (self.z_dim,) model = Sequential() model.add(Dense(256, input_shape=noise_shape)) model.add(LeakyReLU(alpha=0.2)) model.add(BatchNormalization(momentum=0.8)) model.add(Dense(512)) model.add(LeakyReLU(alpha=0.2)) model.add(BatchNormalization(momentum=0.8)) model.add(Dense(1024)) model.add(LeakyReLU(alpha=0.2)) model.add(BatchNormalization(momentum=0.8)) model.add(Dense(np.prod(self.img_shape), activation='tanh')) model.add(Reshape(self.img_shape)) model.summary() return model
***generator*** _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= dense_4 (Dense) (None, 256) 25856 _________________________________________________________________ leaky_re_lu_3 (LeakyReLU) (None, 256) 0 _________________________________________________________________ batch_normalization_1 (Batch (None, 256) 1024 _________________________________________________________________ dense_5 (Dense) (None, 512) 131584 _________________________________________________________________ leaky_re_lu_4 (LeakyReLU) (None, 512) 0 _________________________________________________________________ batch_normalization_2 (Batch (None, 512) 2048 _________________________________________________________________ dense_6 (Dense) (None, 1024) 525312 _________________________________________________________________ leaky_re_lu_5 (LeakyReLU) (None, 1024) 0 _________________________________________________________________ batch_normalization_3 (Batch (None, 1024) 4096 _________________________________________________________________ dense_7 (Dense) (None, 784) 803600 _________________________________________________________________ reshape_1 (Reshape) (None, 28, 28, 1) 0 ================================================================= Total params: 1,493,520 Trainable params: 1,489,936 Non-trainable params: 3,584
詳細
1層目
まず初めの以下の初めの層について見てみる．入力は潜在変数の次元数 *1で、出力は256次元となっている．
noise_shape = (self.z_dim,)
model = Sequential()
model.add(Dense(256, input_shape=noise_shape))
dense_4 (Dense) (None, 256) 25856
ここでなぜパラメータ数が25856になっているかメモして置く．
ここの変換では100次元を256次元に変換させているので100*256=25600個のパラメータなんじゃないかって勘違いしやすいが、バイアスであるBの256次元分のパラメータもあるので25856個のパラメータとなっている．
LeakyReLU(2層目)
活性化関数なので次元数に関しては変化はない.イメージとしては以下の感じ．
BatchNormalization
各層でのアクティベーションの分布を、適切な広がりをもつように調節することができる．また必ずアクティベーション層の後に挿入する必要がある．
利点 - 学習速度が上がる - 重みが初期値に依存しなくなる - 次元数の変化はない
最終ノード
最終ノードの手前でnp.prodを使うことで3次元配列である画像の要素数分の全要素数を算出し出力ノードの次元数に当てている．その後画像に転用できるようにreshapeを行なっている．
Discriminatorモデル
def build_discriminator(self): img_shape = (self.img_rows, self.img_cols, self.channels) print("***build_discriminator***") model = Sequential() model.add(Flatten(input_shape=img_shape)) model.add(Dense(512)) model.add(LeakyReLU(alpha=0.2)) model.add(Dense(256)) model.add(LeakyReLU(alpha=0.2)) model.add(Dense(1, activation='sigmoid')) model.summary() return model
***build_discriminator*** _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= flatten_1 (Flatten) (None, 784) 0 _________________________________________________________________ dense_1 (Dense) (None, 512) 401920 _________________________________________________________________ leaky_re_lu_1 (LeakyReLU) (None, 512) 0 _________________________________________________________________ dense_2 (Dense) (None, 256) 131328 _________________________________________________________________ leaky_re_lu_2 (LeakyReLU) (None, 256) 0 _________________________________________________________________ dense_3 (Dense) (None, 1) 257 ================================================================= Total params: 533,505 Trainable params: 533,505 Non-trainable params: 0
詳細
第1層目
画像の次元数を1次元に変換した要素数分が初めの層のノード数となる．（やっていることとしてはGeneratorの最終層の真逆）
最終層
1次元に変換、この数値は本物か偽物かを判断する0 ~ 1の範囲で連続な確率値を算出する．そのため活性化関数にはsoftmax関数が用いられている．
モデルの結合
結合したモデルに関してはdiscriminatorの重みは固定させた上でgeneratorを学習させる必要がある．
def build_combined1(self): self.discriminator.trainable = False model = Sequential([self.generator, self.discriminator]) return model
訓練
def train(self, epochs=3000, batch_size=128, save_interval=100): print('Load Start') # mnistデータの読み込み (train_images, train_labels),(test_images, test_labels) = mnist.load_data() print('Load Done') # 値を-1 to 1に規格化 train_images = train_images.astype('float32') train_images = train_images/127.5 train_images -= np.ones((train_images.shape)) train_images = np.expand_dims(train_images, axis=3) half_batch = int(batch_size / 2) for epoch in range(epochs): # --------------------- # Discriminatorの学習 # --------------------- # バッチサイズの半数をGeneratorから生成 noise = np.random.normal(0, 1, (half_batch, self.z_dim)) gen_imgs = self.generator.predict(noise) # バッチサイズの半数を教師データからピックアップ #0 ~ train_images.shape[0]未満の乱数をhalf_batch分の配列を用意 idx = np.random.randint(0, train_images.shape[0], half_batch) imgs = train_images[idx] # discriminatorを学習 # 本物データと偽物データは別々に学習させる d_loss_real = self.discriminator.train_on_batch(imgs, np.ones((half_batch, 1))) d_loss_fake = self.discriminator.train_on_batch(gen_imgs, np.zeros((half_batch, 1))) # それぞれの損失関数を平均 d_loss = 0.5 * np.add(d_loss_real, d_loss_fake) # --------------------- # Generatorの学習 # --------------------- noise = np.random.normal(0, 1, (batch_size, self.z_dim)) # 生成データの正解ラベルは本物（1） valid_y = np.array([1] * batch_size) # Train the generator g_loss = self.combined.train_on_batch(noise, valid_y) # 進捗の表示 print ("%d [D loss: %f, acc.: %.2f%%] [G loss: %f]" % (epoch, d_loss[0], 100*d_loss[1], g_loss)) # 指定した間隔で生成画像を保存 if epoch % save_interval == 0: self.sample_images(epoch)
前処理
ループを回す前までに学習画像に対して前処理をかける．
- 画素値に対して0.0~255.0の値のレンジを-1.0~1.0へ変換する．
- channel数を次元に追加している．(*expand_dims)
補足expand_dimsについて
expand_dimsでは第2引数に指定した場所の直前にdim=1を挿入するメソッドである
print(train_images.shape)#(6000, 28, 28) train_images = np.expand_dims(train_images, axis=3) print(train_images.shape)#(6000, 28, 28, 1)
Discriminatorの学習
batchの半数分の100次元乱数をGeneratorに与えたことで生成されたfake画像と、訓練画像の正しい画像ランダムに選択(ミニバッチ学習)を元に学習を行う．
(*batchとは1epochあたりの1まとまりの入力データのことを指す)
Generator
generatorは単体で学習しないのでコンパイルは必要ない
結合モデルの学習
Generatorによって生成されたデータを正しい結果として結合モデルには渡す． Generatorの目的としては偽物を正しいと認識させることで損失関数を少なくする方向へ進めることができるので、生成した画像に対するラベルは1にしている．
疑問点
- train_on_batch
- batch_normalization
参考
*1:, 100
vector配列(C++)
動的にサイズを変更したい場合
vector< vector<int> > arr; //サイズ変更 arr.resize(20); // ()内の数字が要素数になる for( int i=0; i<20; i++ ){ arr[i].resize(10); }
初期化
要素数+初期値(2次元)
// int型の2次元配列(3×4要素の)の宣言 vector<vector<int>> data(3, vector<int>(4));
ユニークな値で初期化
C++11から直接初期化することが可能になった．
std::vector< std::vector<int> > v{{ 1, 1, 1 },{ 2, 2, 2 }};
macではg++ -std=c++11 ~~~.cppみたいな感じでコンパイルするとC++11が有効になってくれる．
参考
C++のvectorまとめ - Qiitaqiita.com
要素の総入れ替え
swapメソッドを使うか、代入でも要素がいれ変わる
std::vector<int> a{1, 2, 3}; std::vector<int> b{10, 9, 8, 7, 6}; a.swap(b);
参照渡し
呼び出し元のvector配列にも反映されるようになる．
void func(vector<int>& v){ for(size_t i=0; i<v.size(); i++){ cout << v[i] << endl; } }
void func2(vector<vector<double>>&A) { for(int i = 0; i < 2; i++){ for(int j = 0; j < 3; j++){ cout << A[i][j] << ", "; A[i][j]+=1; } cout<<endl; } }
newを用いた2次元配列の確保
int size_x = 5; int size_y = 10; double **p = new double*[size_x]; for(int i = 0; i < size_x; i++){ p[i] = new double[size_y]; } for(int i = 0; i < size_x; i ++){ for(int j = 0; j < size_y; j++){ p[i][j] = j + size_y * i; cout << p[i][j] << endl; } }
流体シミュレーション研究ログ No.2
はじめに
今回の記事では障害物があった場合でのポテンシャル線がどのように表示されるかを検証したことについてまとめる．
結果はこんな感じ↓(紫：速度)
場の初期化
検証するフィールドの広さや、障害物の大きさを決定する．
また今回のプログラムではresは離散化した際のブロック数に当たる．（なので格子点の数は+1になるので注意）
obs.size = ofVec2f(150.0, 300.0); obs.mid = ofVec2f(300.0, 0.0); obs.left = obs.mid - ofVec2f(obs.size.x/2.0, 0.0); float eps = 0.000001; rect.size = ofVec2f(800.0, 400.0); rect.res = ofVec2f(100.0, 50.0); rect.left = -1.0 * rect.size / 2.0; rect.delta = rect.size / rect.res;
それぞれの格子点を操作してその点がどのような状況（境界であるのか、それとも障害物内なのかなど）の点であるのかを確認し初期化を行う．
for(int i = 0; i <= rect.res.y; i++){ for(int j = 0; j <= rect.res.x; j++){ int id = i * (rect.res.x + 1) + j; phi.push_back(0.0); psi.push_back(0.0); vel.push_back(ofVec2f(0.0, 0.0)); type.push_back(Type::INSIDE); if(i == 0) type[id] = Type::BOTTOM; if(i == rect.res.y) type[id] = Type::TOP; if(j == 0) type[id] = Type::INLET; if(j == rect.res.x) type[id] = Type::OUTLET; if(j == nMeshX1 && i < nObsWY) type[id] = Type::OBS_LEFT; if(j >= nMeshX1 && j <= nMeshX2 && i == nObsWY) type[id] = Type::OBS_TOP; if(j == nMeshX2 && i < nObsWY) type[id]= Type::OBS_RIGHT; if(j > nMeshX1 && j < nMeshX2 && i < nObsWY) type[id] = Type::OBSTACLE; } } lineNum = 20;
ポテンシャル計算
境界について
http://i.riken.jp/wp-content/uploads/2015/06/secure_4650_080121_toda.pdf
void Environment::calc() { for(int i = 0; i <= rect.res.y; i++) { for(int j = 0; j <= rect.res.x; j++) { int id = j + (rect.res.x + 1) * i; //境界面 if(type[id] == BOTTOM || type[id] == OBS_LEFT || type[id] == OBS_TOP || type[id] == OBS_RIGHT){ psi[id] = 0.0; }else if(type[id] == TOP){ //psi[id] = rect.size.y;//UW p81 psi[id] = 1.0; }else{ float delta = 1.0 / rect.res.y; //psi[id] = i * rect.delta.y;//INLET, OUTLET, INSIDE psi[id] = i * delta;//INLET, OUTLET, INSIDE } if(type[id] == INLET){ phi[id] = 0.0; }else if(phi[id] == OUTLET){ //phi[id] = rect.size.x; phi[id] = 2.0; }else{ float delta = 2.0 / rect.res.x; phi[id] = delta * j; } } }
ラプラスの方程式
連続の式
これは流出、流入が等しいことを示す．（ベクトル演算の内積は流出、流入を示す） $$ {\nabla}{\cdot}\vec{v} = 0 \tag{1} $$
\(\vec{v}\)は(u,v)で構成されている． $$ u = \frac{\partial{\phi}}{\partial{x}}, v = \frac{\partial{\phi}}{\partial{y}} \tag{2} $$
(1)式に(2)式を代入するとラプラス方程式を得る $$ {\nabla}^2{\cdot}{\Phi} = \frac{{\partial}^2{\Phi}} {{\partial}{x}^2} + \frac{{\partial}^2{\Phi}} {{\partial}{y}^2} = 0 \tag{3} $$
差分法
コンピュータ上で数値解析を行う場合、微分演算を差分に変換して解く必要がある．
2階微分の近似
$$ \frac {d^{2}f} {dx^{2}} = \frac{ \frac{df}{dx}_{i+1} - \frac{df}{dx}{_i}}{{\Delta}x} $$
分子の前半は前進差分、後半は後退差分を用いて表すと
$$ \frac{ \frac{df}{dx}_{i+1} - \frac{df}{dx}_i}{{\Delta}x} = \frac{\frac{f_{i+1} - f_i} {{\Delta}x} - \frac{f_i - f_{i-1}} {{\Delta}x}}{{\Delta}x} = \frac{f_{i+1} - 2f_i + f_{i-1}} {{\Delta}x^2} $$
これで差分法をもちいた二階微分の近似方法がわかった．
ラプラス方程式(差分法)
流れ関数に対応するラプラス方程式は
$$ {\nabla^2}{\Psi} = 0 \\ \frac{ {\partial^2}{\Psi} }{ {\partial} {x^2}} + \frac{ {\partial^2}{\Psi} }{ {\partial} {y^2}} = 0 \tag{4} $$
2階微分の部分に対して差分法を適応していくと(4)式は以下のように書き換えることが出来る $$ \frac{{\Psi}_{i+1, j} - 2{\Psi}_{i, j} + {\Psi}_{i-1, j}} {{\Delta}x^2} + \frac{{\Psi}_{i, j+1} - 2{\Psi}_{i, j} + {\Psi}_{i, j-1}} {{\Delta}y^2} = 0\tag{5} $$
\({\Psi_{i,j}}\)を左辺に移動させると
$$ {\Psi_{i,j}} = \frac{ ({\Delta}y)^2 ({\Psi_{i-1,j}} + {\Psi_{i+1,j}} ) + ({\Delta}x)^2 ({\Psi_{i,j-1}} + {\Psi_{i,j+1}} ) } { 2({ ({\Delta}x})^2 { ({\Delta}y})^2) } $$
サンプルコード
ラプラス方程式を具体的に解いたサンプル．注目している点の周囲4点を用いてポテンシャルの値を算出している．
//FDM int iteration = 5000; float tolerance = 0.000001; int cnt = 0; float error = 0.0; float maxError = 0.0; float dx2 = rect.delta.x * rect.delta.x; float dy2 = rect.delta.y * rect.delta.y; float pp; while (cnt < iteration) { maxError = 0.0; for (int i = 1; i < rect.res.y; i++) { for (int j = 1; j < rect.res.x; j++) { int id = j + i * (rect.res.x+1); if(type[id] != INSIDE) continue; //周囲４つ int id1 = j + (i-1) * (rect.res.x + 1);//y-1 int id2 = j + (i+1) * (rect.res.x + 1);//y+1 int id3 = (j-1) + i * (rect.res.x + 1);//x-1 int id4 = (j+1) + i * (rect.res.x + 1);//x+1 //laplatian //psi pp = dy2 * (psi[id3] + psi[id4]) + dx2 *(psi[id1] + psi[id2]); pp /= 2.0 * (dx2 + dy2);
反復法(ガウスザイデル法)
反復法では許容誤差以下になるまで計算を繰り返し行う．
$$
|f{_i}{^n} - f{_i}{^{n-1}}| {\leq} {\epsilon} (i=1,2,N-1)
$$
流体が流れる範囲の全てのポテンシャルを計算し、前のタームとの誤差がその中でも最大となる値をmaxerroeとしている．この値が誤差の許容範囲以下に収まるまで計算を繰り返す．(任意のiで許容誤差内に収まらないといけない)
//さっきのコードの続き //１つ前のψと比較 error = abs(pp - psi[id]); if (error > maxError) maxError = error; psi[id] = pp; } } if (maxError < tolerance) break; cnt++; }
移流方程式とその数値解法(風上差分、CIP法)
移流方程式
自然界には波に関連する現象が多くあります．水面に伝わる波や携帯電話から出る電磁波、人の声(音)、など様々な波の形態で伝搬が行われます．この伝搬（何かしらの物理量が伝わっていくこと）を表す以下のような方程式を移流方程式と言います．
(f : 物理量, u : 流速)
導出
今回は流速が一定の場合を考えていきます． 以下のようにf(x, t)のときの物理量はΔt秒後(次の時間でのサンプリング点)にf(x+cΔt, t+Δt)となります．
（※関数fはt,xの2変数関数なので、図に示す時は以下のように時間を決めたスナップショットのようになる．）
以上より
次の時間での物理量を求める形に変形
一次の項までテイラー展開を行うと
よって
初期条件（時刻 t = 0 のときの値）が与えられれば、この方程式の解（位置と時間から求まる物理量）を求めることができる
数値解法
移流方程式を数値解析的に解くためのいくつかのアルゴリズムをあげるてみると
- 1次風上差分
- Lax-Wendroff法
- CIP法
などがあげられる．
数値拡散
以下の画像のように波が進んでいってそれぞれの離散点だけをみて物理量の値を線形補間すると実際のグラフとの誤差が生まれる．このことを数値拡散という．
(a)の黒い波が進むと青い波に位置に進む．
1次風上差分
導出1
移流方程式
時間微分に関しては前進差分を用いて、位置微分に関しては後退差分を用いる.位置微分に後退差分を用いることで出現するi-1が上流を示すため風上と呼ばれる．
(2)式を変形すると
cをクーラン数という
導出2
黒い波が1ブロック分移動した際の新しい波を青で示している．
地点iにおいて、新しい値を算出したい場合、前の時間での波（黒）の波を元に考えることが出来る．下の図でピンクの直線の方程式は以下のようになる
(2-1)式を元に地点iでの現在の値は先ほど求めた関数Fに移動した分だけ(今回で言うと1ブロック分)戻った位置の値と等しくなるので
導出2の補足
例では微少時間として、移動分は1ブロックとしているがもし移動量が100ブロックなどになった場合上の方程式で得られる値はとても大きな値となり振動してしまう． これを防ぐためクーラン数は1以下でないとう条件が存在する．
CIP法
上記に示したように風上差分だと数値拡散(なまり)が発生してしまうので誤差がつもり、本来の波の形状からどんどん形が変形してしまう．この問題の解消の方法として2点をただ線形補完するのではなく多項式で補完する手法が提案されている．イメージとしてはこのような感じ．
一次風上差分では青のような曲線になるところを多項式補完することで赤の曲線になる．
CIP法では今までの位置だけでなく微分の情報を加えて伝搬させる．以下のようにサンプリング点での粒子の動き（緑の矢印）の情報をク加えて伝搬させてあげることで元の波形に近づけようとしている．
移流方程式(3.1)を微分してもわかるように、微分した値も流速uによって情報が伝搬されている．
補間用の曲線生成
CIP法は3次多項式によるエルミート補間をおこなうので、i-1, i点を３次多項式で表すと
次にi, i-1での値を求める(F{_i}(x)はi, i-1を元に作った曲線のこと)
以上の4つの関係式から係数a, bを求める
次時間n+1の算出
ここまでで補間方程式を求めることが出来た．これを用いて次時刻n+1での値を求めるにはこれまでのように進めた、
、
を求めれば良い．
サンプルコード
void ofApp::methodCIP(){ int ip; double Fp, Gp; double c3, c2, x, s; //流速の方向を元に符号を決定 if(speed > 0.0) { s = 1.0; }else{ s = -1.0; } for(int i = 0; i <= nMesh; i++) { ip = i-(int)s;//上流点 if(ip < 0.0 || ip > nMesh){ Fp = 0.0; Gp = 0.0; }else{ Fp = f0[ip]; Gp = g0[ip]; } double dx = -s * deltaX; double dx2 = dx * dx; double dx3 = dx2 * dx; c3 = (g0[i] + Gp) / dx2 + 2.0 * (f0[i] - Fp) / dx3; c2 = 3.0 * (Fp - f0[i]) / dx2 - (2.0 * g0[i] + Gp) / dx; x = -speed * deltaT;//uΔt //こっちの方がわかりやすいかも //f1[i] = ((c3*x+c2)*x+g0[i])*x + f1[i]; //g1[i] = (3.0 * c3 * x + 2.0 * c2) * x + g1[i]; f1[i] += ((c3*x+c2)*x+g0[i])*x; g1[i] += (3.0 * c3 * x + 2.0 * c2) * x; } //値を更新 for(int i = 0; i <= nMesh; i++) { f0[i] = f1[i]; g0[i] = g1[i]; } }
流体シミュレーション研究ログ No.1
ポテンシャル流れ
『渦度0』の理想的な流れをポテンシャル流れと呼ぶ．このポテンシャル流れを解く為には速度ポテンシャルφと流れ関数ψの２つのスカラ関数が必要になる．
流れ関数
$$ \frac{\partial {\Psi}}{\partial y} = u \tag{1.1} $$
$$ \frac{\partial {\Psi}}{\partial x} = -v\tag{1.2} $$
流れ関数とは
流れ場が\(u(x, y, t), v(x, y, t)\)であるとする．この時ある地点\((x,y)\)について考える． ここに浮きをおいた場合、微小時間後、\((dx, dy)\)だけ移動する． それぞれの移動量は流れ場\(u(x, y, t), v(x, y, t)\)と平行なので微小時間を算出して見ると
$$ \frac{dx}{u} =\frac{dy}{v} $$
つまり $$ -vdx+udy = 0 $$
実際の微少時間はわからないので具体的な移動距離\((dx, dy)\)を知ることはできないが、移動した方向だけ全微分を使うことで推測することが出来る．
先ほどの式を全微分の形に変形していく
$$ d{\Psi} = \frac{{\partial}{\Psi}}{ {\partial}{x} } dx + \frac{{\partial}{\Psi}}{ {\partial}{y} } dy \tag{1.3} $$
つまりこの勾配ベクトルが流速となるスカラー関数\({\Psi(x, y)}\)が流れ関数である .
流れ関数と速度
流れ関数の勾配ベクトルと速度の内積を考える
$$ \boldsymbol{U} {\cdot} {\Delta {\Psi}} = u \frac{\partial {\Psi}}{\partial x} + u \frac{\partial {\Psi}}{\partial y} $$
(1.1), (1.2)より $$ \boldsymbol{U} {\cdot} {\Delta {\Psi}} = -uv + uv \\ = 0 $$
よって流れ関数の勾配方向（等高線にする垂直な方向）と速度ベクトルが垂直、つまり流れ関数の等高線(流線)と速度ベクトルは平行
流線
流れ関数の一定な値を結んだ線(式(1.3)が0となる線上)、また全ての位置での接線は速度ベクトルに一致する
速度ポテンシャル定義
これを満たせば速度ポテンシャルであると言える． $$ \frac{\partial {\Phi}}{\partial x} = u \tag{1.4} $$
$$ \frac{\partial {\Phi}}{\partial y} = v \tag{1.5} $$
速度ポテンシャルと速度
速度ポテンシャルの勾配ベクトルが速度ベクトルに相当するため、速度ポテンシャルの等高線と速度ベクトルは直行する．
速度ベクトルの向きは速度ポテンシャルが増加する方向である．(イメージとは逆なので注意)
また速度ポテンシャルが存在するとき、「回転運動をしない」、「流体中に速度ゼロの所がない」
、「流体全体が動いてる
流れ関数と速度ポテンシャル
上記の導出より、流れ関数と速度ポテンシャルが直行関係にあることがわかる．
複素ポテンシャル
φ,ψは直行しているので、これらの変数で構成されている複素ポテンシャルは以下のようになる．
$$ W = {\Phi} + i{\Psi} $$
一様流れ
一様流れとは速度ベクトルが一定の流れのことを言う．x方向の一定速度を\(U\)とする
$$ U = \frac{d{\phi}} {dx} $$
$$ U = \frac{d{\psi}} {dy} $$
$$ {\phi} = Ux, {\psi} = Uy $$
$$ W = {\phi} + i{\psi} = Ux + iUy = U(x + iy) = Uz $$
ポテンシャル線描画アルゴリズム
目標
すでにグリッド上の点にはポテンシャルの値が算出されている．そのグリッド上でポテンシャルの値が等しい座標を辿って線を引いて行きたい．
到達へのイメージ
四角形の中で出来るだけ滑らかな曲線を引くためには、まず四角形を三角形4つに分解して考える必要がある．
そして『湧き出し』、『吸い込み』などの条件が加わってくるとgrid上のx, y座標それぞれが同じでもポテンシャルの値は変わってくるので一般的に考えると四角形を構成する４点のポテンシャルの値は異なります．
そこで三角形に分解したあとにこの三角形の領域ないでどう線分を引いたら求めたいポテンシャルの値をたどることが出来るかを考える. ①はa,bの比例関係で、②はa,cの比例関係で求まる．これを4箇所で繰り返す感じ．
- まずはじめに表示するポテンシャルの値のmax, minを決める, またポテンシャル線の本数、間隔の値を算出する．
int lineNum = res - 1; float range = 1.0; float minPotential = flowVelocity * -1.0 * range; float maxPotential = flowVelocity * range; //ポテンシャル線一本ごとに変化する値 float dp0 = (maxPotential - minPotential) /(lineNum+1);
- まわり4点のポテンシャルの値をまずはじめに求める．その後ポテンシャルの大小関係を比べて比例式に持ち込んでいる．
for (int k = 0; k < lineNum; k++) { pp = minPotential + (k+1) * dp0; for (int i = 0; i < res-1; i++) { for (int j = 0; j < res-1; j++) { pot[0] = phi[j + i * res ]; pos[0] = ofVec2f(j * rectDelta.x, i * rectDelta.y); pot[1] = phi[j + (i + 1) * res]; pos[1] = ofVec2f(j * rectDelta.x, (i + 1) * rectDelta.y); pot[2] = phi[(j + 1) + (i + 1) * res]; pos[2] = ofVec2f((j + 1)* rectDelta.x, (i + 1) * rectDelta.y); pot[3] = phi[(j + 1) + i * res]; pos[3] = ofVec2f((j + 1)* rectDelta.x, i * rectDelta.y); pot[4] = pot[0]; pos[4] = pos[0]; //center pot[5] = (pot[0] + pot[1] + pot[2] + pot[3]) / 4.0; pos[5] = ofVec2f((pos[0] + pos[1] + pos[2] + pos[3]) / 4.0); for (int m = 0; m < 4; m++) { pos1 = pos2 = ofVec2f(-10, -10); if ((pot[m] <= pp && pp < pot[m + 1]) || (pot[m] > pp && pp >= pot[m + 1])) { pos1 = pos[m] + (pos[m + 1] - pos[m]) * (pp - pot[m]) / (pot[m + 1] - pot[m]); } if ((pot[m] <= pp && pp <= pot[5]) || (pot[m] >= pp && pp >= pot[5])) { if (pos1.x < 0.0) { pos1 = pos[m] + (pos[5] - pos[m]) * (pp - pot[m]) / (pot[5] - pot[m]); } else { pos2 = pos[m] + (pos[5] - pos[m]) * (pp - pot[m]) / (pot[5] - pot[m]); if (type == 0) { ofSetColor(0.0, 0.0, 255.0); }else if (type == 1) { ofSetColor(0.0, 255.0, 0.0); } ofDrawLine(rectStart + pos1, rectStart + pos2); } } if ((pot[m + 1] <= pp && pp <= pot[5]) || (pot[m + 1] >= pp && pp >= pot[5])) { if (pos1.x < 0.0) { pos1 = pos[m + 1] + (pos[5] - pos[m + 1]) * (pp - pot[m + 1]) / (pot[5] - pot[m + 1]); } else { pos2 = pos[m + 1] + (pos[5] - pos[m + 1]) * (pp - pot[m + 1]) / (pot[5] - pot[m + 1]); ofDrawLine(rectStart + pos1, rectStart + pos2); } } }//m }//j }//i }//k
参考
流体方程式の解き方入門 www.amazon.co.jp
『流れ』と『波』のシミュレーション www.amazon.co.jp
PythonとKerasによるDeepLearning個人的なまとめNo.03(CNN)
全結合での内部の行列演算
Denseの第1引数は変換あとの次元数に相当する．
inputs = Input(shape=(28,28)) # shape =(28,28) x1 = Flatten()(inputs) # shape (28,28)=>(784) x2 = Dense(512, activation='relu')(x1) # shape (784)=>(512) x3 = Dense(512, activation='relu')(x2) # shape (512)=>(512) y = Dense(10, activation='softmax')(x3)
行ベクトルが１つのサンプルを構成する特徴ベクトル
CNN
畳み込みでの行列計算
Conv2D第１引数は深さの次元数、第2引数は窓の大きさを表している．
inputs = Input(shape=(28,28,1)) # shape =(28,28,1) x1 = Conv2D(32, (3, 3), activation='relu')(inputs) # shape (28,28,1)=>(26,26,32) x2 = Conv2D(64, (3, 3), activation='relu')(x1) # shape (26,26,32)=>(24,24,64)
入力に対して得られたマップのことを応答マップと言う．
パラメータの数の算出の仕方
畳み込みフィルターサイズを(L,L)(L,L)、入力チャンネルをMM、出力チャンネルをNNとすると直方体変形には(L⋅L⋅M⋅N+N)(L⋅L⋅M⋅N+N)がパラメータ数
Poolingと畳み込み
最大化Pooloingでは基本的にスライド=2, ウィンドウ=2x2 それに対して畳み込みではストライド=0、ウィンドウ=2x2 . またPoolingでは最大値が残る．
最大化Poolingを行うことでのメリット
パラメータを削減出来る
flattenで平滑化した後の全結合の部分でのパラメータ数が全然違うことがわかる．これはMaxPoolingしなかった場合、平滑化した際のベクトルの次元が30976になりそれを64次元に変換するために30976*64+64のパラメータが必要になり膨大、、、MaxPoolingすることでこれを防ぐことが出来る．（ちなみに活性化関数はReLu）
圧縮性
MaxPoolingしていないモデルでは初めの入力が(28, 28)で最終的に(22, 22)になっているので、(28,28)に(7,7)のウィンドウをかけているような物である．それに対してMaxPoolingしている方では(28,28)が(11, 11)になっていて(18,18)のウィンドウで見ていることになる．どれだけの情報が最後の層に圧縮されているかイメージしてみれば後者の方が優れていることがわかる．
出力フィルターサイズの計算の仕方
畳み込みを行なった後のFilterのサイズの計算の仕方 - In : 入力サイズ - Pad : Paddingのサイズ - F : 畳こみのフィルターサイズ - S : ストライドのサイズ(畳み込みを行う際どれだけfilterを移動させるか)
$$ \frac {In + 2Pad - F}{S} $$
参考
すごく参考になった↓ https://qiita.com/nishiha/items/bfd5dfcd7fffd3c529bc
PythonとKerasによるDeepLearning個人的なまとめNo.02(回帰分析)
回帰分析
モチベーション
13の特徴量を持つデータの住宅の住宅価格の中央値を求めたい
データロード
訓練データは404sampleで特徴は13次元存在する．この特徴を元に中央値が決定される
from keras.datasets import reuters from keras.utils.np_utils import to_categorical import numpy as np import keras from keras import models from keras import layers from keras.datasets import boston_housing (train_data, train_targets), (test_data, test_targets) = boston_housing.load_data()
model構築
初めの層のinputは特徴の数分の13次元用意しておいてあとはどんなサンプル数がきても対応できるように空白にしておく
def build_model(): model = models.Sequential() model.add(layers.Dense(64, activation='relu', input_shape=(train_data.shape[1], ))) model.add(layers.Dense(64, activation='relu')) model.add(layers.Dense(1)) model.compile(optimizer='rmsprop', loss='mse', metrics=['mae']) return model
k分割法
sampleをk等分にする.
#sampleをk等分にする k = 4 num_val_samples = len(train_data) // k for i in range(k): print('processing fold', i) #val_dataの要素数はnum_val_samples val_data = train_data[i * num_val_samples : (i+1) * num_val_samples] val_targets = train_targets[i * num_val_samples : (i+1) * num_val_samples] #data内のval_data以外の部分に当たる prac_train_data = np.concatenate( [train_data[:i * num_val_samples], train_data[(i+1) * num_val_samples:]], axis=0) prac_train_targets = np.concatenate( [train_targets[:i * num_val_samples], train_targets[(i+1) * num_val_samples:]], axis=0) model= build_model() model.fit(prac_train_data, prac_train_targets, epochs=num_epoch, batch_size=1, verbose=0)
