{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Gaussian process\n",
    "\n",
    "PHYSBO performs Bayesian optimization while running Gaussian process regression.\n",
    "\n",
    "Therefore, it is possible to run Gaussian process regression given training data, and to predict test data using the trained model.\n",
    "\n",
    "In this section, the procedure is introduced.\n",
    "\n",
    "\n",
    "## Preparation of search candidate data\n",
    "\n",
    "In this tutorial, the problem of finding a stable interface structure for Cu is used as an example. The values that have already been evaluated are used, although the evaluation of the objective function, i.e., the structural relaxation calculation, actually takes on the order of several hours per calculation. For more information on the problem setup, please refer to the following references\n",
    "\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",
    "\n",
    "Save the dataset file [s5-210.csv](https://raw.githubusercontent.com/issp-center-dev/PHYSBO/master/examples/grain_bound/data/s5-210.csv) into the subdirectory `data`, and load dataset from this file as the following:"
   ]
  },
  {
   "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": [
    "## Defining training data\n",
    "\n",
    "A randomly selected 10% of the target data will be used as training data, and another randomly selected 10% will be used as test data."
   ]
  },
  {
   "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": [
    "## Learning and Prediction of Gaussian Processes\n",
    "\n",
    "The following process is used to learn the Gaussian process and predict the test data.\n",
    "\n",
    "1. Generate a model of the Gaussian process\n",
    "\n",
    "2. The model is trained using X_train (parameters of the training data) and t_train (objective function value of the training data).\n",
    "\n",
    "3. Run predictions on the test data (X_test) using the trained model.\n",
    "\n",
    "Definition of covariance (Gaussian)"
   ]
  },
  {
   "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": [
    "Definition of mean value"
   ]
  },
  {
   "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": [
    "Definition of likelihood function (Gaussian)"
   ]
  },
  {
   "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": [
    "Generation of a Gaussian Process Model"
   ]
  },
  {
   "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": [
    "Learning a Gaussian process model."
   ]
  },
  {
   "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 17312.31220145003\n",
      "50 -th epoch marginal likelihood 6291.292745798703\n",
      "100 -th epoch marginal likelihood 3269.1167759139516\n",
      "150 -th epoch marginal likelihood 1568.3930580794922\n",
      "200 -th epoch marginal likelihood 664.2847129159145\n",
      "250 -th epoch marginal likelihood -249.28468708456558\n",
      "300 -th epoch marginal likelihood -869.7604930929888\n",
      "350 -th epoch marginal likelihood -1316.6809532065581\n",
      "400 -th epoch marginal likelihood -1546.1623851368954\n",
      "450 -th epoch marginal likelihood -1660.7298135295766\n",
      "500 -th epoch marginal likelihood -1719.5056128528097\n",
      "Done\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gp.fit(X_train, t_train, config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Output the parameters in the learned Gaussian process."
   ]
  },
  {
   "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.81666924]\n",
      "mean parameter in GP prior:  [-1.05939674]\n",
      "covariance parameter in GP prior:  [-0.91578975 -2.45544347]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "gp.print_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Calculating the mean (predicted value) and variance of the test data"
   ]
  },
  {
   "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": [
    "Results of prediction"
   ]
  },
  {
   "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.00420815, -1.10923758, -0.97840623, ..., -1.00323733,\n",
       "       -0.97015759, -1.11076236])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Results of covariance"
   ]
  },
  {
   "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.00056069, 0.00075529, 0.00043006, ..., 0.0016925 , 0.00070103,\n",
       "       0.00073499])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Output mean square error of prediction"
   ]
  },
  {
   "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.008107085662147708"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((fmean-t_test)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "##  Prediction by trained models\n",
    "\n",
    "Read the parameters of the trained model as gp_params and make predictions using them.\n",
    "\n",
    "By storing gp_params and training data (X_train, t_train), prediction by the trained model is possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Prepare the learned parameters (must be done immediately after learning)"
   ]
  },
  {
   "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.81666924, -1.05939674, -0.91578975, -2.45544347])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Prepare the learned parameters as a 1D array\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": [
    "Prepare a model similar to the one used for training as 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": [
    "#Definition of covariance (Gaussian)\n",
    "cov = physbo.gp.cov.gauss( X_train.shape[1],ard = False )\n",
    "\n",
    "#Definition of mean value\n",
    "mean = physbo.gp.mean.const()\n",
    "\n",
    "#Definition of likelihood function (Gaussian)\n",
    "lik = physbo.gp.lik.gauss()\n",
    "\n",
    "#Generation of a Gaussian Process Model\n",
    "gp = physbo.gp.model(lik=lik,mean=mean,cov=cov)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Prepare a model similar to the one used for training as gp"
   ]
  },
  {
   "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": [
    "#Input learned parameters into the Gaussian process.\n",
    "gp.set_params(gp_params)\n",
    "\n",
    "\n",
    "#Calculate the mean (predicted value) and variance of the test data\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": [
    "Results of prediction"
   ]
  },
  {
   "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.00420815, -1.10923758, -0.97840623, ..., -1.00323733,\n",
       "       -0.97015759, -1.11076236])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fmean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Results of covariance"
   ]
  },
  {
   "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.00056069, 0.00075529, 0.00043006, ..., 0.0016925 , 0.00070103,\n",
       "       0.00073499])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Output mean square error of prediction"
   ]
  },
  {
   "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.008107085662147708"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((fmean-t_test)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "Note: In the example above, we used the same pre-registered X to make predictions.\n",
    "If you want to make predictions for parameters X_new that is not included in X using the trained model,\n",
    "you first obtain the mean (X_{mean}) and standard deviation (X_{std}) of the data X and\n",
    "then normalize X_{new} by (X_{new} - X_{mean}) / X_{std}.\n",
    "Also, the data format for X is ndarray format.\n",
    "Therefore, if X_{new} is a single data, it must be transformed to ndarray format.\n",
    "For example, if X_{new} is a real number, you should replace X_new as\n",
    "X_new = np.array(X_new).reshape(1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}