{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# ガウス過程\n",
    "\n",
    "PHYSBOではガウス過程回帰を実行しながらベイズ最適化を行なっています。\n",
    "\n",
    "そのため、学習データが与えられた際にガウス過程回帰を実行することもでき、学習済みモデルを利用したテストデータの予測も行うことができます。\n",
    "\n",
    "ここでは、その手順について紹介します。\n",
    "\n",
    "\n",
    "## 探索候補データの準備\n",
    "\n",
    "本チュートリアルでは例として、Cuの安定した界面構造の探索問題を扱います。 目的関数の評価にあたる構造緩和計算には、実際には1回あたり数時間といったオーダーの時間を要しますが、本チュートリアルでは既に評価済みの値を使用します。問題設定については、以下の文献を参照してください。\n",
    "S. Kiyohara, H. Oda, K. Tsuda and T. Mizoguchi, “Acceleration of stable interface structure searching using a kriging approach”, Jpn. J. Appl. Phys. 55, 045502 (2016).\n",
    "\n",
    "データセットファイル [s5-210.csv](https://raw.githubusercontent.com/issp-center-dev/PHYSBO/master/examples/grain_bound/data/s5-210.csv) を `data` ディレクトリ以下に保存し、次のように読み出します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:41.987250Z",
     "start_time": "2020-12-04T06:11:41.537168Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import physbo\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "def load_data():\n",
    "    A =  np.asarray(np.loadtxt('data/s5-210.csv',skiprows=1, delimiter=',') )\n",
    "    X = A[:,0:3]\n",
    "    t  = -A[:,3]\n",
    "    return X, t\n",
    "\n",
    "X, t = load_data()\n",
    "X = physbo.misc.centering( X )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 学習データの定義\n",
    "\n",
    "対象データのうち、ランダムに選んだ1割をトレーニングデータとして利用し、別のランダムに選んだ1割をテストデータとして利用します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:51.077070Z",
     "start_time": "2020-12-04T06:11:51.072211Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ntrain = 1798\n",
      "Ntest = 1798\n"
     ]
    }
   ],
   "source": [
    "N = len(t)\n",
    "Ntrain = int(N*0.1)\n",
    "Ntest = min(int(N*0.1), N-Ntrain)\n",
    "\n",
    "id_all   = np.random.choice(N, N, replace=False)\n",
    "id_train  = id_all[0:Ntrain]\n",
    "id_test = id_all[Ntrain:Ntrain+Ntest]\n",
    "\n",
    "X_train = X[id_train]\n",
    "X_test = X[id_test]\n",
    "\n",
    "t_train = t[id_train]\n",
    "t_test = t[id_test]\n",
    "\n",
    "print(\"Ntrain =\", Ntrain)\n",
    "print(\"Ntest =\", Ntest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## ガウス過程の学習と予測\n",
    "\n",
    "以下のプロセスでガウス過程を学習し、テストデータの予測を行います。\n",
    "\n",
    "1. ガウス過程のモデルを生成します。\n",
    "\n",
    "2. X_train（学習データのパラメータ）, t_train（学習データの目的関数値）を用いてモデルを学習します。\n",
    "\n",
    "3. 学習されたモデルを用いてテストデータ(X_test)に対する予測を実行します。\n",
    "\n",
    "共分散の定義(ガウシアン)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:55.403677Z",
     "start_time": "2020-12-04T06:11:55.399915Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "cov = physbo.gp.cov.gauss( X_train.shape[1],ard = False )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "平均の定義"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:56.279543Z",
     "start_time": "2020-12-04T06:11:56.277082Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mean = physbo.gp.mean.const()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "尤度関数の定義(ガウシアン)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:57.077507Z",
     "start_time": "2020-12-04T06:11:57.075581Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "lik = physbo.gp.lik.gauss()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "ガウス過程モデルの生成"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:11:57.832602Z",
     "start_time": "2020-12-04T06:11:57.828902Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "gp = physbo.gp.model(lik=lik,mean=mean,cov=cov)\n",
    "config = physbo.misc.set_config()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "ガウス過程モデルを学習"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.218792Z",
     "start_time": "2020-12-04T06:11:58.261609Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start the initial hyper parameter searching ...\n",
      "Done\n",
      "\n",
      "Start the hyper parameter learning ...\n",
      "0 -th epoch marginal likelihood 16051.69452976001\n",
      "50 -th epoch marginal likelihood 4551.39626443153\n",
      "100 -th epoch marginal likelihood 2141.377872729846\n",
      "150 -th epoch marginal likelihood 595.8361411907399\n",
      "200 -th epoch marginal likelihood -373.23922309413774\n",
      "250 -th epoch marginal likelihood -929.2472009209496\n",
      "300 -th epoch marginal likelihood -1273.8727959019732\n",
      "350 -th epoch marginal likelihood -1413.2553901394206\n",
      "400 -th epoch marginal likelihood -1477.3889625983586\n",
      "450 -th epoch marginal likelihood -1525.339082571432\n",
      "500 -th epoch marginal likelihood -1539.787541261617\n",
      "Done\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gp.fit(X_train, t_train, config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "学習されたガウス過程におけるパラメタを出力"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.227479Z",
     "start_time": "2020-12-04T06:12:58.221821Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "likelihood parameter =   [-2.57036368]\n",
      "mean parameter in GP prior:  [-1.0654197]\n",
      "covariance parameter in GP prior:  [-0.59460765 -2.4232173 ]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gp.print_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "テストデータの平均値（予測値）および分散を計算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.605713Z",
     "start_time": "2020-12-04T06:12:58.244883Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "gp.prepare(X_train, t_train)\n",
    "fmean = gp.get_post_fmean(X_train, X_test)\n",
    "fcov = gp.get_post_fcov(X_train, X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "予測の結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.618218Z",
     "start_time": "2020-12-04T06:12:58.607794Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.01865696, -0.98356729, -0.97797072, ..., -0.99978278,\n",
       "       -0.98145533, -0.9956255 ])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "分散の結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.628483Z",
     "start_time": "2020-12-04T06:12:58.622345Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00046688, 0.0010622 , 0.0006136 , ..., 0.00043492, 0.0005969 ,\n",
       "       0.00053435])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "予測値の平均二乗誤差の出力"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.636081Z",
     "start_time": "2020-12-04T06:12:58.631461Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.004179032574484333"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((fmean-t_test)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 訓練済みモデルによる予測\n",
    "\n",
    "学習済みモデルのパラメタをgp_paramsとして読み出し、これを用いた予測を行います。\n",
    "\n",
    "gp_paramsおよび学習データ(X_train, t_train)を記憶しておくことで、訓練済みモデルによる予測が可能となります。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "学習されたパラメタを準備(学習の直後に行う必要あり)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.645968Z",
     "start_time": "2020-12-04T06:12:58.639012Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.57036368, -1.0654197 , -0.59460765, -2.4232173 ])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#学習したパラメタを１次元配列として準備\n",
    "gp_params =  np.append(np.append(gp.lik.params, gp.prior.mean.params), gp.prior.cov.params)\n",
    "\n",
    "gp_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "学習に利用したモデルと同様のモデルをgpとして準備"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:58.666019Z",
     "start_time": "2020-12-04T06:12:58.653259Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "#共分散の定義 (ガウシアン)\n",
    "cov = physbo.gp.cov.gauss( X_train.shape[1],ard = False )\n",
    "\n",
    "#平均の定義\n",
    "mean = physbo.gp.mean.const()\n",
    "\n",
    "#尤度関数の定義 (ガウシアン)\n",
    "lik = physbo.gp.lik.gauss()\n",
    "\n",
    "#ガウス過程モデルの生成\n",
    "gp = physbo.gp.model(lik=lik,mean=mean,cov=cov)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "学習済みのパラメタをモデルに入力し予測を実行"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:59.016429Z",
     "start_time": "2020-12-04T06:12:58.673034Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "#学習済みのパラメタをガウス過程に入力\n",
    "gp.set_params(gp_params)\n",
    "\n",
    "\n",
    "#テストデータの平均値（予測値）および分散を計算\n",
    "gp.prepare(X_train, t_train)\n",
    "fmean = gp.get_post_fmean(X_train, X_test)\n",
    "fcov = gp.get_post_fcov(X_train, X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "予測の結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:59.020795Z",
     "start_time": "2020-12-04T06:12:59.017606Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.01865696, -0.98356729, -0.97797072, ..., -0.99978278,\n",
       "       -0.98145533, -0.9956255 ])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "分散の結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:59.026523Z",
     "start_time": "2020-12-04T06:12:59.023035Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00046688, 0.0010622 , 0.0006136 , ..., 0.00043492, 0.0005969 ,\n",
       "       0.00053435])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "予測値の平均二乗誤差の出力"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-04T06:12:59.033497Z",
     "start_time": "2020-12-04T06:12:59.027871Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.004179032574484333"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((fmean-t_test)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "(注) 上の例では事前に登録されているXと同じものを利用して予測を行いました。\n",
    "学習済みのモデルを利用してXに含まれていないパラメータ X_new に対して予測をしたい場合には、\n",
    "学習モデルで使用したデータXの平均(X_{mean})と標準偏差(X_{std})を求めていただいたうえで、\n",
    "X_{new} = (X_{new} - X_{mean}) / X_{std}\n",
    "の変形を行うことで予測を行うことができます。\n",
    "また、渡す際のデータ形式はndarray形式になっています。\n",
    "そのため、X_{new}が一つのデータの場合には事前に変換する必要があります。\n",
    "例えば、X_{new}が実数である場合には、\n",
    "X_new = np.array(X_new).reshape(1)\n",
    "などとして変換する必要があります。"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}