OpenMDAOでKriging応答曲面

# Kriging応答曲面
理論的にはよく理解してはいないが,散らばりのあるサンプル点から曲面を推定する手法である.
[OpenMDAOでラテン超方格サンプリング](http://qiita.com/TatsuyaKatayama/items/7bcf863bad444c6ad901)に示したラテン超方格サンプリングと相性が良い(散らばりのあるサンプリング)ため,工学的にはCAEを利用して応答曲面を推定する際しばしば利用される.

# OpenMDAOで放物面のKriging応答曲面
取組みの概略は次のとおりである。まず[OpenMDAOでラテン超方格サンプリング](http://qiita.com/TatsuyaKatayama/items/7bcf863bad444c6ad901)で実験したデータ（doe_paraboloidファイル）を読み込む.
次に読み込んだ実験データを用いてKriging近似モデルを鍛え,作成する.
最後に,作成したKriging近似モデルによる内挿値を理論解(ワイヤーフレーム)プロットすることで,モデルの妥当性を確認する.尚,次式が理論解である.

```math
\begin{align}
& f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 \\
    {\rm subject \: to} \: \: \:&  -50.0\leq x \leq 50.0 \\
                                &  -50.0\leq y \leq 50.0
\end{align}
```
---
#Groupの準備
下記のtraining_mm.pyファイルを準備する

```py:training_mm.py
from __future__ import print_function
from openmdao.api import Group,  MetaModel,  FloatKrigingSurrogate

class TrainingMM(Group):
    ''' FloatKriging gives responses as floats '''

    def __init__(self):
        super(TrainingMM, self).__init__()

        # Create meta_model for f_x as the response
        mm = self.add("parabo_mm", MetaModel())
        mm.add_param('x', val=0.)
        mm.add_param('y', val=0.)
        mm.add_output('f_xy', val=0., surrogate=FloatKrigingSurrogate())
```

MetaModelというの近似モデルの訓練と値の評価を行うComponentをを追加している.
f_xyを近似するためにFloatKrigingSurrogate というSurrogateModelをセットしている.
FloatKrigingSurrogate単体でも近似モデルの訓練と値の評価を行う機能を有しているが,
MetaModelを用いれば複数の近似モデルを同時に訓練できると、公式ドキュメントに書かれている.

---
#Problem(問題）の設定

下記kriging_paraboloid.pyを準備する.前半は実験データの読み込み～Kriging近似モデルの訓練までの処理である.

```py:kriging_paraboloid.py
#! /bin/python
import pickle
import numpy as np
import sqlitedict
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from openmdao.api import Problem
from training_mm import TrainingMM

db =sqlitedict.SqliteDict("doe_paraboloid","iterations")
res = np.array([[0.,0.,0.]] * len(db.items()))
for i, v in enumerate(db.items()):
    res[i,0] = v[1]["Unknowns"]["x"]
    res[i,1] = v[1]["Unknowns"]["y"]
    res[i,2] = v[1]["Unknowns"]["f_xy"]

prob = Problem()
prob.root = TrainingMM()
prob.setup()
prob["parabo_mm.train:x"] = res[:,0]
prob["parabo_mm.train:y"] = res[:,1]
prob["parabo_mm.train:f_xy"] = res[:,2]
prob.run()
```

---
後半部を下記に示す。作成したkriging近似モデルと理論解の比較をプロットにて行っている

```py:kriging_paraboloid.py続き
x = np.arange(-50., 50., 4.)
y = np.arange(-50., 50., 4.)

sg = prob.root.parabo_mm._init_unknowns_dict["f_xy"]["surrogate"]
f = open("./kriging_parabo","w")
pickle.dump(sg, f)
f.close()
f = open("./kriging_parabo_mm","w")
pickle.dump(prob.root.parabo_mm, f)
f.close()


xyzkrig = np.array([[xi,yi,sg.predict(np.array([xi,yi]))[0]] \
        for xi in x for yi in y])
#xyzkrig = np.array([[0.,0.,0.]]*(25*25))
#cnt = 0
#for xi in x:
#    for yi in y:
#        xyzkrig[cnt,0] = prob["parabo_mm.x"] = xi
#        xyzkrig[cnt,1] = prob["parabo_mm.y"] = yi
#        prob.run()
#        xyzkrig[cnt,2] = prob["parabo_mm.f_xy"]
#        cnt += 1

fig = plt.figure(figsize=(6,4)); ax = Axes3D(fig)
X, Y = np.meshgrid(x, y)
Z = (X-3.0)**2. + X*Y + (Y+4.0)**2. - 3.0
ax.plot_wireframe(X,Y,Z,label="theoretical surface")
ax.scatter3D(xyzkrig[:,0], xyzkrig[:,1], xyzkrig[:,2], c="r", marker="o",label="Interpolated value")
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('f_xy')
plt.legend()
plt.show()
```

---
kriging_paraboloid.pyの前半部は解説は省略し後半部のみ解説する.
まず冒頭ではKriging近似モデルを評価するサンプル点を作成している.
これはLHS(ラテン超方格サンプリング)した点とは別の点であり, 
$-50 \leq x,y \leq 50$の範囲を25分割した25x25のサンプル点にて評価している. 

```sg = prob.root.parabo_mm._init_unknowns_dict["f_xy"]["surrogate"]```でプライベートクラス変数(_init_unknowns_dict)にアクセスして,MeataModelに設置されているsurrogateモデルを読み込んでいる.
褒められた処置ではないがこれには２点理由がある.
* 鍛えた近似モデルの保存が実装されていない.今回はpickeleでシリアライズした.
* MetaModelで近似モデルの値を評価しようとすると少し煩雑になる(```#xyzkrig = ・・・```以下のコメント行）.

残りはプロット処理である.
理論解をワイヤーフレームにKriging近似モデルによる評価値をドットでプロットしている.

# Kriging近似モデルによる値の内挿
LHSのサンプル数が多すぎたこともあり,かなり精度はよさそう.

---
 
![kriging_paraboloid.png](https://qiita-image-store.s3.amazonaws.com/0/132571/1b5d822a-3e4f-a97d-7799-e8c7e69242a9.png)