PAM_guidedSynthesis-regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import xgboost as xgb\n",
    "import shap\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn import metrics \n",
    "from sklearn.model_selection import KFold, GridSearchCV\n",
    "\n",
    "import time\n",
    "import random\n",
    "\n",
    "from utils import data_handler, plotter\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "\n",
    "#handle warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "np.random.seed(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('font', family='serif')\n",
    "\n",
    "# set font size\n",
    "SMALL_SIZE = 22\n",
    "MEDIUM_SIZE=24\n",
    "BIGGER_SIZE = 26\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=MEDIUM_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metrics: \n",
    "r2,r,MSE - coefﬁcient of determination (R2), Pearson coefﬁcient (r), and mean square error (MSE) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(y_pred,y_true,end_ptr):\n",
    "    '''Evaluate the performance of IAM through 3 metrics; \n",
    "       And find the best predicted condition for the next IAM trial\n",
    "       \n",
    "       Inputs:\n",
    "           X_test      - for testing\n",
    "           y_true      - ground true for X_test\n",
    "           end_ptr     - range of [0,end_ptr] in y_pred and y_true are values observable to IAM\n",
    "           \n",
    "       Uses:(global variables)    \n",
    "           pred_result - of shape (total_samp,), storing all the results predicted by IAM\n",
    "           true_results- of shape (total_samp,), the same order as of pred_results, storing all the ground truth\n",
    "           \n",
    "       Outputs:\n",
    "            \n",
    "    \n",
    "    '''\n",
    "    #All sample metrics    \n",
    "    r2 = metrics.r2_score(y_true,y_pred)   \n",
    "    mse = metrics.mean_squared_error(y_true,y_pred)\n",
    "    [pear,p_value] = pearsonr(y_true,y_pred)\n",
    "    \n",
    "    #metrics based on observable samples\n",
    "    y_true_s = y_true[0:(end_ptr+1)]\n",
    "    y_pred_s = y_pred[0:(end_ptr+1)]\n",
    "    r2_s = metrics.r2_score(y_true_s,y_pred_s)   \n",
    "    mse_s = metrics.mean_squared_error(y_true_s,y_pred_s)\n",
    "    [pear_s,p_value_s] = pearsonr(y_true_s,y_pred_s)\n",
    "    return [r2,pear,p_value,mse,r2_s,pear_s,p_value_s,mse_s] #, best_pos_ind, best_prob \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Handling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df,Y_df = data_handler.load_XY(1)\n",
    "X = X_df.as_matrix()\n",
    "Y = Y_df.as_matrix()/100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up & construct initial dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cross validation settup\n",
    "inner_nsplits = 10\n",
    "init_train_size = 20\n",
    "totalSamp = X.shape[0]\n",
    "\n",
    "\n",
    "Y_global_max = np.max(Y)\n",
    "all_ind = np.random.permutation(list(range(0,totalSamp)))\n",
    "all_ind_wo_max = list(range(0,totalSamp))\n",
    "all_ind_wo_max.remove(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PAM guided sythesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup\n",
    "Nc = 0\n",
    "save_csv= True\n",
    "verbose=True\n",
    "title='cqd_PAM_'\n",
    "batch=1\n",
    "to_break = False # whether to break when the best condition is found\n",
    "\n",
    "## start IAM guided synthesis...\n",
    "init_time = time.time()\n",
    "\n",
    "#construct initial training set\n",
    "results_mat = np.zeros(((totalSamp-init_train_size),12))\n",
    "\n",
    "\n",
    "train_ind = random.sample(all_ind_wo_max, init_train_size)\n",
    "test_ind = [x for x in all_ind if x not in train_ind]\n",
    "\n",
    "# set up result storage to compute eval metrics, in the order of IAM\n",
    "#  ignore the initial training set, as it is not determined by IAM\n",
    "pred_results = np.zeros(totalSamp-init_train_size)\n",
    "true_results = np.zeros(totalSamp-init_train_size)\n",
    "\n",
    "\n",
    "# setup the hyperparameter range for tuning\n",
    "tuned_parameters = dict(learning_rate=[0.01],\n",
    "                        n_estimators=[300,500,700], #100,,300,400,500\n",
    "                        colsample_bylevel = [0.5,0.7,0.9],\n",
    "                      gamma=[0,0.2], #0,0.1,0.2,0.3,0.4\n",
    "                      max_depth =[3,7,11], # [3,7,11]]\n",
    "                      reg_lambda = [0.1,1,10], #[0.1,1,10]\n",
    "                     # reg_alpha = [1],\n",
    "                       subsample=[0.4,0.7,1])\n",
    "\n",
    "j=0\n",
    "loop_count = 0\n",
    "mean_y_only_init = np.mean(Y[train_ind])\n",
    "std_y_only_init = np.std(Y[train_ind])\n",
    "\n",
    "\n",
    "while(j<totalSamp-init_train_size):        \n",
    "    inner_cv = KFold(n_splits=inner_nsplits,shuffle=True, random_state=j) \n",
    "    X_train = X[train_ind]\n",
    "    Y_train = Y[train_ind]\n",
    "    X_test = X[test_ind]\n",
    "    Y_test = Y[test_ind]\n",
    "\n",
    "    last_max = np.max(Y_train)\n",
    "\n",
    "    # GradientBoost\n",
    "    reg = xgb.XGBRegressor(objective=\"reg:linear\",min_child_weight=1,**{'tree_method':'exact'},\n",
    "                             silent=True,n_jobs=4,random_state=3,seed=3);\n",
    "\n",
    "    gb_clf = GridSearchCV(reg,tuned_parameters, cv=inner_cv,scoring='r2',verbose=0,n_jobs=4)\n",
    "    gb_clf.fit(X_train, Y_train)\n",
    "    y_pred = gb_clf.predict(X_test)\n",
    "\n",
    "\n",
    "    # choose the condition with best predicted yield\n",
    "    best_pos_ind = np.argsort(-y_pred)[:batch]       \n",
    "    best_prob = y_pred[best_pos_ind]\n",
    "    next_ind = np.array(test_ind)[best_pos_ind]        \n",
    "\n",
    "    # update results storage\n",
    "    train_size = len(Y_train)\n",
    "    temp = list(range(0,len(y_pred)))\n",
    "    ind_notbest = [x for x in temp if x not in best_pos_ind]\n",
    "\n",
    "    start_ptr = j\n",
    "    end_ptr = np.min([start_ptr + batch, totalSamp-init_train_size]) \n",
    "    pred_results[start_ptr:end_ptr] = best_prob\n",
    "    pred_results[end_ptr:totalSamp-init_train_size]= y_pred[ind_notbest]\n",
    "    true_results[start_ptr:end_ptr] = Y_test[best_pos_ind]\n",
    "    true_results[end_ptr:totalSamp-init_train_size] = Y_test[ind_notbest]\n",
    "\n",
    "    pred_metrics = test(pred_results,true_results,end_ptr-1)    \n",
    "\n",
    "    # calculate results\n",
    "    next_best_true_ind = next_ind[np.argmax(Y_test[best_pos_ind])]\n",
    "    next_best_y_true = np.max(Y_test[best_pos_ind])\n",
    "    result_list = [train_size,next_best_true_ind,next_best_y_true,best_prob[0],] + pred_metrics   \n",
    "    results_mat[loop_count,:] = np.array(result_list)\n",
    "\n",
    "    loop_count = loop_count +1\n",
    "    j = j + batch\n",
    "\n",
    "    if(verbose):\n",
    "        print(loop_count,'->',j,', best_next_ind=',next_best_true_ind, ' best_Y_true=',\"{0:.6f}\".format(next_best_y_true),' train_max=',\"{0:.6f}\".format(last_max),' r2=',pred_metrics[0])\n",
    "\n",
    "    train_ind = [*train_ind , *next_ind]   \n",
    "    test_ind = [x for x in test_ind if x not in next_ind]\n",
    "\n",
    "    ## critical point\n",
    "    if(next_best_y_true==Y_global_max):\n",
    "        Nc = j+init_train_size\n",
    "        print('***Nc=', str(Nc))\n",
    "        if(to_break):\n",
    "            break\n",
    "\n",
    "\n",
    "results = pd.DataFrame(data=results_mat[0:j,:],columns=['sample_size','pred_ind','y_true','best_pred_result','r2','pearson','p_value','mse','r2_s','pearson_s','p_value_s','mse_s'])\n",
    "\n",
    "if(save_csv):\n",
    "    data_handler.save_csv(results,title=title+'_Nc_'+str(Nc))\n",
    "\n",
    "\n",
    "run_time = (time.time() - init_time)/60\n",
    "\n",
    "print('end at loop ',j)\n",
    "print('critical point is when size of training set = ',Nc)\n",
    "print('total time = ',run_time,' mins')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samp = results['sample_size']\n",
    "pred_results = results['best_pred_result'] \n",
    "true_results = results['y_true'] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 3 # window size\n",
    "\n",
    "seq_len = len(pred_results)\n",
    "pad = int((w-1)/2)\n",
    "pred_norm = np.convolve(pred_results, np.ones(w)*1/w , mode='valid')#medfilt(pred_results,w)[pad:seq_len-pad]##max_pooling(pred_results, w=9 , s=1)#\n",
    "#pred_norm = medfilt(pred_norm,3)#[1:len(pred_norm)-1]\n",
    "true_norm = np.convolve(true_results, np.ones(w)*1/w , mode='valid')#max_pooling(true_results, w=w , s=1)#np.convolve(true_results, np.array([0.2, 0.2, 0.2, 0.2, 0.2]) , mode='valid')#\n",
    "n_samp_new = n_samp[pad:seq_len-pad]\n",
    "\n",
    "# save csv\n",
    "to_save_mat = np.zeros((len(n_samp_new), 3))\n",
    "to_save_mat[:,0]  = n_samp_new\n",
    "to_save_mat[:,1] = pred_norm\n",
    "to_save_mat[:,2] = true_norm\n",
    "\n",
    "to_save_df = pd.DataFrame(to_save_mat, columns=['n_sample','pred','true'])\n",
    "save_path = data_handler.save_csv(to_save_df, title=title+'_smoothed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (12,12/1.618),dpi=100)\n",
    "#y = Y[train_ind]\n",
    "plt.plot(n_samp_new,pred_norm,lw=2, label='Best Predicted', linestyle='-',color='b')\n",
    "plt.plot(n_samp_new,true_norm,lw=2, label='Best True', linestyle='--',color='r')\n",
    "#labels = list(range(30,np.max(n_samp),50)) + [n_samp]\n",
    "plt.axvline(x=Nc,linestyle='--',color='k')\n",
    "#plt.xticks(labels)  \n",
    "\n",
    "plt.grid(False)\n",
    "\n",
    "plt.legend(loc='lower left')\n",
    "plt.xlabel('Number of Explored Conditions')\n",
    "plt.ylabel('PLQY')\n",
    "plt.show()\n",
    "\n",
    "save_path = plotter.save_fig(fig, title+'_smoothed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "end_time = time.time()\n",
    "print('total time = ',(end_time - init_time)/60,' mins')"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}