PAM_guidedSynthesis-classification.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",
    "\n",
    "  \n",
    "#from sklearn.preprocessing import StandardScaler\n",
    "from sklearn import metrics \n",
    "from sklearn.model_selection import GridSearchCV,StratifiedKFold\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "import time\n",
    "init_time = time.time()\n",
    "\n",
    "from utils import data_handler, plotter\n",
    "\n",
    "\n",
    "#handle warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(clf, X_test,y_test):\n",
    "    y_pred = clf.predict(X_test)\n",
    "    c_mat = metrics.confusion_matrix(y_test,y_pred)\n",
    "    \n",
    "    best_pos_ind = 0\n",
    "    fn = 0\n",
    "    fp = 0\n",
    "    if(len(c_mat)==1):# all correctly predicted as one class\n",
    "        best_pos_ind = -1\n",
    "        if(y_test[0]==0): # all predicted as negative class\n",
    "            tpr_ts = float('nan')\n",
    "            tnr_ts = 1\n",
    "            acc_ts = 1\n",
    "        if(y_test[0]==1): #all predicted as positive class\n",
    "            tpr_ts = 1\n",
    "            tnr_ts = float('nan')   \n",
    "            acc_ts = 1\n",
    "    else:\n",
    "        fn = c_mat[1,0]\n",
    "        fp = c_mat[0,1]\n",
    "        tpr_ts = c_mat[1,1]/(c_mat[1,0]+c_mat[1,1])#TP/(TP+FN)\n",
    "        tnr_ts =  c_mat[0,0]/(c_mat[0,0]+c_mat[0,1])\n",
    "        acc_ts = (c_mat[0,0] + c_mat[1,1])/c_mat.sum()\n",
    "    \n",
    "    pred_prob = clf.predict_proba(X_test)  \n",
    "\n",
    "    best_prob = float('nan')\n",
    "    if(best_pos_ind != -1):\n",
    "        best_pos_ind = np.argmax(pred_prob[:,1])\n",
    "        best_prob = pred_prob[best_pos_ind,1]\n",
    "        \n",
    "    results = [acc_ts,tpr_ts,tnr_ts]\n",
    "    return results, best_pos_ind, best_prob , fp, fn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Handling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_df,Y_df = data_handler.load_XY()\n",
    "X = X_df.as_matrix()\n",
    "Y = Y_df.as_matrix()\n",
    "feature_list = X_df.columns\n",
    "unique, counts = np.unique(Y, return_counts=True) #unique, counts = numpy.unique(a, return_counts=True)\n",
    "tot_can_count = counts[1]\n",
    "tot_cnot_count = counts[0]"
   ]
  },
  {
   "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",
    "totalSamp = X.shape[0]\n",
    "verbose= True\n",
    "\n",
    "\n",
    "def generate_init_sets():\n",
    "    '''\n",
    "         construct initial training/testing dataset\n",
    "         to make sure each class has at least $inner_nsplits$ samples\n",
    "    '''\n",
    "\n",
    "    can_counter = 0\n",
    "    cnot_counter = 0\n",
    "    \n",
    "    # shuffle indexes of data samples\n",
    "    Y_global_max = np.max(Y)\n",
    "    all_ind = np.random.permutation(list(range(0,totalSamp)))\n",
    "    train_ptr = 0\n",
    "\n",
    "    while(can_counter < inner_nsplits or cnot_counter< inner_nsplits):\n",
    "        \n",
    "        next_ind = all_ind[train_ptr]\n",
    "        train_ptr = train_ptr+1\n",
    "        \n",
    "        if(Y[next_ind] ==1):\n",
    "            can_counter = can_counter +1            \n",
    "        else:\n",
    "            cnot_counter = cnot_counter + 1\n",
    "\n",
    "    ret_dict = {}\n",
    "    ret_dict['train_ind'] = list(all_ind[0:train_ptr])\n",
    "    ret_dict['test_ind']  = list(all_ind[train_ptr:len(all_ind)])\n",
    "    return ret_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PAM guided sythesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#critical point \n",
    "Nc = 0\n",
    "\n",
    "# setup initial sets\n",
    "init_sets = generate_init_sets()\n",
    "train_ind = init_sets['train_ind']\n",
    "test_ind = init_sets['test_ind']\n",
    "if(verbose):\n",
    "    print(train_ind)  \n",
    "\n",
    "# Results store \n",
    "init_train_size = len(train_ind)\n",
    "init_cnot_count = list(Y[train_ind]).count(0) \n",
    "init_can_count = init_train_size - init_cnot_count\n",
    "results_mat = np.zeros((totalSamp-init_train_size,8))\n",
    "\n",
    "# setup hyperparameter range to tune\n",
    "tuned_parameters = dict(learning_rate=[0.01],#0.01,0.1,0.2,0.3\n",
    "                      n_estimators=[100,300,500], #100\n",
    "                      gamma=[0,0.2,0.4], #0,0.1,0.2,0.3,0.4\n",
    "                      max_depth =[5,7,9,11], # [4,5,6]\n",
    "                      reg_lambda = [0.1,1,10], \n",
    "                        colsample_bylevel = [0.9],\n",
    "                        subsample=[0.4,0.7,1])\n",
    "\n",
    "\n",
    "# start PAM guided synthesis...\n",
    "for j in range(totalSamp):\n",
    "    inner_cv = StratifiedKFold(n_splits=inner_nsplits, shuffle=True,random_state=j) #StratifiedKFold(n_splits=inner_nsplits, 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",
    "    #count pos/neg of training set\n",
    "    tr_zero_count = list(Y_train).count(0)\n",
    "    tr_total_count = len(train_ind)\n",
    "    pos_tr = tr_total_count - tr_zero_count\n",
    "\n",
    "    # GradientBoost\n",
    "    pipe = xgb.XGBClassifier(objective='binary:logistic',min_child_weight=1,**{'tree_method':'exact'},\n",
    "                             silent=True,n_jobs=4,random_state=3,seed=3, scale_pos_weight=1);\n",
    "\n",
    "\n",
    "    gb_clf = GridSearchCV(pipe,tuned_parameters, cv=inner_cv,scoring='roc_auc',verbose=0,n_jobs=4)\n",
    "    gb_clf.fit(X_train, Y_train)\n",
    "    result_list, next_ind, best_prob,fp_ts, fn_ts = test(gb_clf,X_test,Y_test)\n",
    "\n",
    "\n",
    "    # calculate results\n",
    "    type1_err = (fp_ts + tr_zero_count - init_cnot_count) / (tot_cnot_count - init_cnot_count)\n",
    "    type2_err = fn_ts / (tot_can_count - init_can_count)              \n",
    "    results_mat[j,:] = np.array([tr_total_count] + result_list + [best_prob ,pos_tr, type1_err, type2_err])\n",
    "\n",
    "    next_ind = test_ind[next_ind]\n",
    "    if(verbose):\n",
    "        print(j,'loop, next_ind=',next_ind, ' #tr=',tr_total_count,' pos_tr=',pos_tr,' best_prob=',\"{0:.6f}\".format(best_prob),' type1=',\"{0:.6f}\".format(type1_err),' type2=',\"{0:.6f}\".format(type2_err))\n",
    "\n",
    "    # critical point\n",
    "    if((best_prob <0.5) and (Nc == 0)):            \n",
    "        Nc = tr_total_count            \n",
    "\n",
    "    #stopping condition\n",
    "    if(pos_tr == tot_can_count):\n",
    "        break\n",
    "\n",
    "    #update train/test sets\n",
    "    train_ind = train_ind + [next_ind]      \n",
    "    test_ind.remove(next_ind)\n",
    "\n",
    "print('end at loop ',j, '  Nc = ',Nc)\n",
    "\n",
    "\n",
    "results = pd.DataFrame(data=results_mat[0:j+1],columns=['sample_size','acc_ts','tpr_ts','tnr_ts','best_prob','pos_tr','type1_err','type2_err'])\n",
    "saved_title = data_handler.save_csv(results,title='mos2_PAM_results_Nc_'+str(Nc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samp = results['sample_size']\n",
    "#\n",
    "acc = results['acc_ts']\n",
    "tpr = results['tpr_ts']\n",
    "tnr = results['tnr_ts']\n",
    "#\n",
    "type1_err= results['type1_err']\n",
    "type2_err = results['type2_err']\n",
    "#\n",
    "best_p = results['best_prob']\n",
    "#\n",
    "pos_tr = results['pos_tr']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import medfilt\n",
    "def smooth(in_seq, w = 5):\n",
    "    \n",
    "    out_seq = medfilt(in_seq,w)\n",
    "    \n",
    "    return out_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (12,12/1.618),dpi=1000)\n",
    "\n",
    "w=7\n",
    "pad = int((w-1)/2)\n",
    "n_samp_new = n_samp\n",
    "\n",
    "best_p_smoothed = smooth(best_p, w = w)\n",
    "\n",
    "plt.plot(n_samp_new,best_p_smoothed,lw=2, label='Best Predicted', linestyle='-',color='b')\n",
    "#plt.plot(n_samp_new,best_p,lw=2, label='Highest \"Can grow\"\\nProbability', linestyle='--',color='r')#marker='^', \n",
    "labels = list(range(30,int(Nc),50))\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('\"Can grow\" Probability')\n",
    "plt.show()\n",
    "\n",
    "save_path = plotter.save_fig(fig,'mos2_PAM_a_')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (12,12/1.618))#,dpi=1000\n",
    "plt.plot(n_samp_new,smooth(tnr, w = w),lw=1.5, label='True Negative Rate',  linestyle='-',color='g',alpha=0.8)#marker='^',\n",
    "plt.plot(n_samp_new,smooth(acc, w = w),lw=1.5, label='Accuracy',  linestyle='-',color='y',alpha=0.8)#marker='.',\n",
    "plt.plot(n_samp_new,smooth(tpr, w = w),lw=1.5, label='True Positive Rate', linestyle='-',color='b',alpha=0.8)#marker='^', \n",
    "\n",
    "\n",
    "labels = list(range(30,int(np.max(n_samp)),50))\n",
    "plt.axvline(x=Nc,linestyle='--',color='k')\n",
    "plt.xticks(labels)  \n",
    "\n",
    "plt.grid(False)\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Number of Explored Conditions')\n",
    "plt.ylabel('Performance on Unexplored Conidtions')\n",
    "plt.show()\n",
    "\n",
    "save_path = plotter.save_fig(fig,'mos2_PAM_b_')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (12,12/1.618))#,dpi=1000\n",
    "\n",
    "plt.plot(n_samp,pos_tr,label='Number of \"Can grow\" conditions explored',  lw=2, linestyle='-',color='b')\n",
    "plt.plot([31,210],[21,181],lw=3, linestyle=':',color='r')\n",
    "\n",
    "labels = list(range(30,int(np.max(n_samp)),50))\n",
    "plt.axvline(x=Nc,linestyle='--',color='k')\n",
    "plt.xticks(labels)  \n",
    "#xticklabels(labels) \n",
    "\n",
    "plt.grid(False)\n",
    "#plt.plot([0,1,2,3],[1,1,1,1],lw=2, label='label',marker='.', linestyle=':', color='k')\n",
    "plt.legend(loc='lower right')\n",
    "plt.xlabel('Number of Explored Conditions')\n",
    "plt.ylabel('Count')\n",
    "plt.show()\n",
    "\n",
    "save_path = plotter.save_fig(fig,'mos2_PAM_c_')"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}