はじめに

2017年現在CAE業界では、Deep Learningの活用はこれからという段階
- 複雑な非線形現象にNNを用いる(NNとカップリング）
- 解析結果の予測
- DQNによる最適化
今回はCFDのあるタイムステップ時の解析結果を生成することを考える
前回の無残な記事
~~CNNで定常流れの解析結果を計算はじめ数ステップから定常解を予測できなかった件~~ --> 再トライ中

CFD環境

ご挨拶

import PyFoam
import PyFoam.FoamInformation

pyVer = sys.version
pfVer = PyFoam.versionString()
ofVer =  PyFoam.FoamInformation.foamVersion()

print("Python version: %s" %pyVer[:6])
print("PyFoam version: %s" %pfVer)
print("OpenFOAM version: %d.%s" %ofVer)

Python version: 3.5.3 
PyFoam version: 0.6.8.1
OpenFOAM version: 4.1

Deep Leaning環境

TensorFlow Googleが開発しているDeep Leaningライブラリ

pip show tensorflow

Name: tensorflow
Version: 1.1.0
Summary: TensorFlow helps the tensors flow
Home-page: http://tensorflow.org/
Author: Google Inc.
Author-email: opensource@google.com
License: Apache 2.0
Location: /home/ubuntu/.pyenv/versions/anaconda3-4.3.1/envs/py35/lib/python3.5/site-packages
Requires: wheel, numpy, six, protobuf, werkzeug

TensorFlowフロントエンド環境

keras
Pythonで書かれたDeep Leaningライブラリ
バックエンドでMXNet、Deeplearning4j、TensorFlow、CNTK、Theanoが動作(from Wikipedia)

pip show keras

Name: Keras
Version: 2.0.2
Summary: Deep Learning for Python
Home-page: https://github.com/fchollet/keras
Author: Francois Chollet
Author-email: francois.chollet@gmail.com
License: MIT
Location: /home/ubuntu/.pyenv/versions/anaconda3-4.3.1/envs/py35/lib/python3.5/site-packages
Requires: theano, pyyaml, six

CFD対象(前回の記事と同じ)

OpenFOAMのチュートリアルからsimpleFOAM/pitzDailyを改造(メッシュを直交化)　末尾コードを参照
乱流モデルは標準k-epsilon

ビッグデータの準備

各cellのdx,dy,ccx,ccyを計算
t=1～10ステップ,最終ステップのp,Ux,Uy,k,epsilon,nutを取得
対象のcellを中心に9x9のcellに関する3,4のデータを取得。tはラベルとして使用

CVEA (Conditional Variational Autoencoder)

810(9x9x10)次元のベクトルから810次元の分散ベクトルと平均ベクトルを作成し、元の810次元のベクトルに戻す
分散・平均のベクトルから元の810次元ベクトルに戻すときにラベルを与えてあげることでラベルに対応した元ベクトルを生成できる
コード含めこちらの記事を大いに参考にさせていただく
こちらの記事も参考にさせていただいた

CVEAソース

必要ライブラリのインポート

import warnings
import numpy as np
from keras.layers import Input, Dense, Lambda, Dropout
from keras.layers.noise import GaussianNoise
from keras.layers.merge import concatenate as concat
from keras.models import Model
from keras import backend as K
from keras.datasets import mnist
from keras.utils import to_categorical
from keras.callbacks import EarlyStopping
from keras.optimizers import Adam
import matplotlib.pyplot as plt

warnings.filterwarnings('ignore')
%pylab inline

学習用のデータと検証用のデータを準備

X = np.load("./X2.npy")
X = X.astype(np.float32)
X = X.reshape(12225,9,9,16,10)
np.random.seed(1)
n_trains = 60000
n_ylabel = 0
n_t = 16
t = 0
c = 0
X_train = np.ones(n_trains * X[0,:,:,0,n_ylabel:].size,dtype=np.float32) * 10
X_train = X_train.reshape(n_trains,
                          X[c,:,:,t,n_ylabel:].shape[0],
                          X[c,:,:,t,n_ylabel:].shape[1],
                          X[c,:,:,t,n_ylabel:].shape[2])

Y_train = np.ones(n_trains * (1 + n_ylabel),dtype=np.float32).reshape(n_trains, (1 + n_ylabel)) * 10

CandT_train = np.ones(2*n_trains, dtype=np.float32).reshape(-1,2) * -1
for i in range(n_trains):
    t = np.random.randint(0,n_t)
    c = np.random.randint(0,12225)
    CandT_train[i,0] = c
    CandT_train[i,1] = t

    X_train[i,:,:,:] = X[c,:,:,t,n_ylabel:]
    Y_train[i,0] = t
n_test = 10000

X_test = np.ones(n_test * X[0,:,:,0,n_ylabel:].size,dtype=np.float32) * 10
X_test = X_test.reshape(n_test,
                          X[c,:,:,t,n_ylabel:].shape[0],
                          X[c,:,:,t,n_ylabel:].shape[1],
                          X[c,:,:,t,n_ylabel:].shape[2])

Y_test = np.ones(n_test * (1 + n_ylabel),dtype=np.float32).reshape(n_test, (1 + n_ylabel)) * 10


CandT_test = np.ones(2*n_test, dtype=np.float32).reshape(-1,2) * -1
for i in range(n_test):
    t = np.random.randint(0,n_t)
    c = np.random.randint(0,12225)
    CandT_test[i,0] = c
    CandT_test[i,1] = t

    X_test[i,:,:,:] = X[c,:,:,t,n_ylabel:]
    Y_test[i,0] = t
n_test2 = 100000

X_test2 = np.ones(n_test2 * X[0,:,:,0,n_ylabel:].size,dtype=np.float32) * 10
X_test2 = X_test2.reshape(n_test2,
                          X[c,:,:,t,n_ylabel:].shape[0],
                          X[c,:,:,t,n_ylabel:].shape[1],
                          X[c,:,:,t,n_ylabel:].shape[2])

Y_test2 = np.ones(n_test2 * (1 + n_ylabel),dtype=np.float32).reshape(n_test2, (1 + n_ylabel)) * 10

CandT_test2 = np.ones(2*n_test2, dtype=np.float32).reshape(-1,2) * -1
for i in range(n_test2):
    t = np.random.randint(0,n_t)
    c = np.random.randint(0,12225)
    CandT_test2[i,0] = c
    CandT_test2[i,1] = t

    X_test2[i,:,:,:] = X[c,:,:,t,n_ylabel:]
    Y_test2[i,0] = t

n_pixels = np.prod(X_train.shape[1:])
X_train = X_train.reshape((len(X_train), n_pixels)) * 0.5 +0.5
X_test = X_test.reshape((len(X_test), n_pixels)) * 0.5 +0.5
X_test2 = X_test2.reshape((len(X_test2), n_pixels)) * 0.5 +0.5

y_train = to_categorical(Y_train)
y_test = to_categorical(Y_test)
y_test2 = to_categorical(Y_test2)

print("X_train:%s" % (X_train.shape,))
print("Y_train:%s" %(Y_train.shape,))
print("X_test:%s" %(X_test.shape,))
print("Y_test:%s" %(Y_test.shape,))
print("X_test2:%s" %(X_test2.shape,))
print("Y_test2:%s" %(Y_test2.shape,))

(60000, 9, 9, 10) (60000, 1)

ハイパーパラメータ

m = 20 # batch size
n_z = 2 # latent space size
encoder_dim1 = 512 # dim of encoder hidden layer
decoder_dim = 512 # dim of decoder hidden layer
decoder_out_dim = n_pixels # dim of decoder output layer
activ = 'relu'
optim = Adam(lr=0.001)
n_x = X_train.shape[1]
n_y = y_train.shape[1]
n_epoch = 100

モデル作成

X = Input(shape=(n_x,))
label = Input(shape=(n_y,))
inputs = concat([X, label])
encoder_h = Dense(encoder_dim1, activation=activ)(inputs)
mu = Dense(n_z, activation='linear')(encoder_h)
l_sigma = Dense(n_z, activation='linear')(encoder_h)
def sample_z(args):
    mu, l_sigma = args
    eps = K.random_normal(shape=(m, n_z), mean=0., stddev=1.)
    return mu + K.exp(l_sigma / 2) * eps

# Sampling latent space
z = Lambda(sample_z, output_shape = (n_z, ))([mu, l_sigma])

z = Lambda(sample_z, output_shape = (n_z, ))([mu, l_sigma])

# merge latent space with label
zc = concat([z, label])

decoder_hidden = Dense(decoder_dim, activation=activ)
decoder_out = Dense(decoder_out_dim, activation='sigmoid')
h_p = decoder_hidden(zc)
outputs = decoder_out(h_p)

def vae_loss(y_true, y_pred):
    recon = K.sum(K.binary_crossentropy(y_true, y_pred), axis=-1)
    kl = 0.5 * K.sum(K.exp(l_sigma) + K.square(mu) - 1. - l_sigma, axis=-1)
    return recon + kl

def KL_loss(y_true, y_pred):
    return(0.5 * K.sum(K.exp(l_sigma) + K.square(mu) - 1. - l_sigma, axis=1))

def recon_loss(y_true, y_pred):
    return K.sum(K.binary_crossentropy(y_true, y_pred), axis=-1)
cvae = Model([X, label], outputs)
encoder = Model([X, label], mu)

d_in = Input(shape=(n_z+n_y,))
d_h = decoder_hidden(d_in)
d_out = decoder_out(d_h)
decoder = Model(d_in, d_out)

cvae.compile(optimizer=optim, loss=vae_loss, metrics = [KL_loss, recon_loss])

cvae.summary()

____________________________________________________________________________________________________
Layer (type)                     Output Shape          Param #     Connected to                     
====================================================================================================
input_1 (InputLayer)             (None, 810)           0                                            
____________________________________________________________________________________________________
input_2 (InputLayer)             (None, 16)            0                                            
____________________________________________________________________________________________________
concatenate_1 (Concatenate)      (None, 826)           0                                            
____________________________________________________________________________________________________
dense_1 (Dense)                  (None, 512)           423424                                       
____________________________________________________________________________________________________
dense_2 (Dense)                  (None, 2)             1026                                         
____________________________________________________________________________________________________
dense_3 (Dense)                  (None, 2)             1026                                         
____________________________________________________________________________________________________
lambda_2 (Lambda)                (None, 2)             0                                            
____________________________________________________________________________________________________
concatenate_2 (Concatenate)      (None, 18)            0                                            
____________________________________________________________________________________________________
dense_4 (Dense)                  (None, 512)           9728                                         
____________________________________________________________________________________________________
dense_5 (Dense)                  (None, 810)           415530                                       
====================================================================================================
Total params: 850,734.0
Trainable params: 850,734.0
Non-trainable params: 0.0
____________________________________________________________________________________________________

モデルをプロットする

from keras.utils import plot_model
fname = "pitzDaily_cvea.png"
plot_model(cvae, to_file='./%s'% fname , show_shapes=True)

下記のようなモデル図が出力される。

学習

GPUじゃないから遅い。。。
GPU買いたいが買えない。

from keras_tqdm import TQDMNotebookCallback
cvae_hist = cvae.fit([X_train, y_train], X_train, verbose = 0, batch_size=m, epochs=n_epoch,
                            validation_data = ([X_test, y_test], X_test),
                            callbacks = [TQDMNotebookCallback(),EarlyStopping(patience = 5)])

TQDMNotebookCallbackを用いてjupyter上でプログレスバーを表示させていたので、学習過程をstdout出力させなかった。
TQDMNotebookCallback()を使う場合も、verbose = 1にした方が良いかも

学習履歴

学習履歴をプロットする。loss関数が下がってきている。
しっかり学習できてる？

import pandas as pd
dfhist = pd.DataFrame(cvae_hist.history)

ax1 = dfhist.loss.plot(legend=True)
ax1 = dfhist.recon_loss.plot(ax=ax1,legend=True)
ax1 = dfhist.val_loss.plot(ax=ax1,legend=True)
ax1 = dfhist.val_recon_loss.plot(ax=ax1,legend=True)
ax1.legend(loc='upper center')
ax2 = ax1.twinx()
ax2 = dfhist.KL_loss.plot(ax=ax2,legend=True, c="b")
ax2 = dfhist.val_KL_loss.plot(ax = ax2,legend=True,  c="c")
ax2.legend(loc='upper right')
plt.savefig("./data2/cvae_hist_plot.png")

encoderの分布

平均・分散のベクトルの分布図を描く

あまりに鮮やかな画像。これは流れ場を特徴を表しているのか？

学習がうまくいっていかいのか？
X_test2をつかって相関図を描く

予測精度(相関図)

ぐは・・・

いまいちな結果に。。。
本来はここからラベルt=iの結果からi+1を予測しようと考えていたが、、、見直してから実施する

まとめ

黒魔術を唱えるためにも最低限のレベルが必要
- 問題設定が悪い？
- 正規分布を仮定して平均分散ベクトル化しているが隠れ層が足りない？
- ビッグデータとしてどうなの？

CVAEで流体解析のあるタイムステップ時の解析結果を生成できなかった件