PD06_OPT_ANN_SEL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "###########import packages##########\n",
    "import tensorflow as tf\n",
    "import keras\n",
    "from keras import optimizers\n",
    "from keras import regularizers\n",
    "from keras import backend as K\n",
    "from keras.callbacks import EarlyStopping\n",
    "from keras.callbacks import TensorBoard\n",
    "from keras.constraints import max_norm\n",
    "from keras.models import Sequential \n",
    "from keras.layers import Dense \n",
    "from keras.layers import Dropout \n",
    "from keras.models import Model\n",
    "from keras.layers import BatchNormalization\n",
    "from keras.wrappers.scikit_learn import KerasClassifier \n",
    "from keras.wrappers.scikit_learn import KerasRegressor\n",
    "from keras.constraints import maxnorm \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xgboost as xgb\n",
    "import lightgbm\n",
    "import catboost\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from xgboost import plot_importance\n",
    "from xgboost import XGBClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn import ensemble\n",
    "from sklearn.tree import ExtraTreeRegressor\n",
    "from sklearn import svm\n",
    "from sklearn import neighbors\n",
    "from sklearn import tree\n",
    "from sklearn.impute import SimpleImputer\n",
    "from keras.callbacks import EarlyStopping \n",
    "from sklearn.metrics import r2_score\n",
    "### ljy改5：限制显存\n",
    "gpus = tf.config.experimental.list_physical_devices('GPU')  # 获取GPU列表\n",
    "for gpu in gpus:\n",
    "    tf.config.experimental.set_memory_growth(gpu, True)\n",
    "    # 失效： tf.config.experimental.set_per_process_memory_fraction(0.25)\n",
    "    # 第一个参数为原则哪块GPU，只有一块则是gpu[0],后面的memory_limt是限制的显存大小，单位为M\n",
    "    tf.config.experimental.set_virtual_device_configuration(gpu, [tf.config.experimental.VirtualDeviceConfiguration(memory_limit=1024*8)]) \n",
    "\n",
    "early_stopping=keras.callbacks.EarlyStopping(\n",
    " monitor=\"val_loss\", \n",
    " patience=20, \n",
    " verbose=0, \n",
    " mode=\"auto\"\n",
    ")\n",
    "%matplotlib\n",
    "###########wrapping root mean square error for later calls##########\n",
    "def compute_mae_mse_rmse(target,prediction):\n",
    "    error = []\n",
    "    for i in range(len(target)):\n",
    "        error.append(target[i] - prediction[i])\n",
    "    squaredError = []\n",
    "    absError = []\n",
    "    for val in error:\n",
    "        squaredError.append(val * val)  # target-prediction之差平方\n",
    "        absError.append(abs(val))  # 误差绝对值\n",
    "    mae=sum(absError)/len(absError)  # 平均绝对误差MAE\n",
    "    mse=sum(squaredError)/len(squaredError)  # 均方误差MSE\n",
    "    RMSE=np.sqrt(sum(squaredError)/len(absError))\n",
    "    R2=r2_score(target,prediction)\n",
    "    return mae,mse,RMSE,R2\n",
    "###########loading data##########\n",
    "fdata=pd.read_csv('datatest.csv',encoding=\"gbk\")\n",
    "raw_data=fdata.loc[:,[      \n",
    "                      'Dielectric Constant (C  N-1 M-2)',#0\n",
    "                      'Flash Point(℃)',#1\n",
    "                      'Stirring Method (0 for Mechanical Stirring 1 for Ultrasound)',#2\n",
    "                      'Flow Rate (mL min-1)',#3\n",
    "                      'Water Content (wt%)',#4\n",
    "                      'Anodic Platinum Loading Amount (mgPt cm-2)',#5\n",
    "                      'Cathodic Platinum Loading Amount (mgPt cm-2)',#6\n",
    "                      'I/C',#7\n",
    "                      'Cell Temperature (℃)',#8\n",
    "                      'Anode Flow Rate (sccm)',#9\n",
    "                      'Cathode Flow Rate (sccm)',#10\n",
    "                      'Solid Content (wt%)',#11\n",
    "                      'Backpressure (Mpa)',#12\n",
    "                      'Power Density at 0.6V (mW cm-2)'#13\n",
    "                        ]]\n",
    "\n",
    "###########data standardization##########\n",
    "standardized_data = (raw_data-np.mean(raw_data,axis=0))/np.std(raw_data,axis=0)\n",
    "\n",
    "###########defining a wrapper function for later call from each machine learning algorithms##########\n",
    "raw_input=standardized_data.iloc[:,0:13]\n",
    "raw_output=standardized_data.iloc[:,13]\n",
    "X=raw_input.values.astype(np.float32)\n",
    "y=raw_output.values.astype(np.float32)\n",
    "###########fix random seed for reproducability##########\n",
    "seed=91\n",
    "###########train test splitting##########\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.15,random_state=seed)\n",
    "raw_input_global=raw_data.iloc[:,0:13]\n",
    "raw_output_global=raw_data.iloc[:,13]\n",
    "###########wrap up fuction for later call for OPTIMIZATION##########\n",
    "def evaluate(pre_2,real_2):\n",
    "    pre_2=np.array(pre_2)\n",
    "    real_2=np.array(real_2)\n",
    "    pre_2_series=pd.Series(pre_2)\n",
    "    real_2_series=pd.Series(real_2)\n",
    "    return rmse(pre_2,real_2), round(pre_2_series.corr(real_2_series), 3)\n",
    "def compare(list_name,limit):\n",
    "    judge=1\n",
    "    for a in list_name:\n",
    "        if a < limit:\n",
    "            judge=judge*1\n",
    "        else:\n",
    "            judge=judge*0\n",
    "    return judge\n",
    "def generate_arrays_from_file(path):\n",
    "    while True:\n",
    "        with open(path) as f:\n",
    "            for line in f:\n",
    "                # create numpy arrays of input data\n",
    "                # and labels, from each line in the file\n",
    "                x1, x2, y = process_line(line)\n",
    "                yield ({'input_1': x1, 'input_2': x2}, {'output': y})\n",
    "for batch_size_number in [16,24,32]:\n",
    "    for reg in [0,0.0001,0.001]:\n",
    "        for dropout_rate in [0,0.1,0.2,0.3,0.4,0.5]:\n",
    "            for neurons1 in range(200,1000,100):\n",
    "                for epochs_number in range(150,850,100):\n",
    "                    for learning_rate_search in [0.0005,0.00075,0.001,0.0025,0.005,0.0075,0.01,0.02]:\n",
    "                        for activation1 in ['relu']:\n",
    "                            regularizer=keras.regularizers.l2(reg)\n",
    "                            ###########keras ANN model construction##########\n",
    "                            model = Sequential() \n",
    "                            model.add(Dense(neurons1, input_dim=13, kernel_initializer='random_normal',\n",
    "                                            bias_initializer='random_normal',activation=activation1,kernel_regularizer=regularizer)) \n",
    "                            model.add(Dropout(dropout_rate))\n",
    "                            model.add(Dense(neurons1, input_dim=neurons1, kernel_initializer='random_normal',\n",
    "                                            bias_initializer='random_normal',activation=activation1,kernel_regularizer=regularizer)) \n",
    "                            model.add(Dropout(dropout_rate))\n",
    "                            model.add(Dense(1, input_dim=neurons1, activation='linear'))\n",
    "                            adam=optimizers.Adam(lr=learning_rate_search)\n",
    "                            model.compile(loss='mse', optimizer=adam)\n",
    "                            print('training...')\n",
    "                            model.fit(X_train, y_train,verbose=0, batch_size=batch_size_number,epochs=epochs_number,validation_split=0.2,callbacks=[early_stopping])\n",
    "                            result=model.predict(X_test)\n",
    "                            result_train=model.predict(X_train)\n",
    "                            ###########get RMSE and R2 on the test set##########\n",
    "                            x_prediction_07=result*np.std(raw_output_global,axis=0)+np.mean(raw_output_global,axis=0)\n",
    "                            y_real_07=np.std(raw_output_global,axis=0)*y_test+np.mean(raw_output_global,axis=0)\n",
    "                            x_prediction_07_series=pd.Series(x_prediction_07[:,0])\n",
    "                            y_real_07_series=pd.Series(y_real_07)\n",
    "                            #training set\n",
    "                            x_prediction_07_train=result_train*np.std(raw_output_global,axis=0)+np.mean(raw_output_global,axis=0)\n",
    "                            y_real_07_train=np.std(raw_output_global,axis=0)*y_train+np.mean(raw_output_global,axis=0)\n",
    "                            x_prediction_07_series_train=pd.Series(x_prediction_07_train[:,0])\n",
    "                            y_real_07_series_train=pd.Series(y_real_07_train)\n",
    "                            ###########evaluating the regression quality##########\n",
    "                            corr_ann = round(x_prediction_07_series.corr(y_real_07_series), 5)\n",
    "                            error_val= compute_mae_mse_rmse(x_prediction_07[:,0],y_real_07)\n",
    "                            corr_ann_train = round(x_prediction_07_series_train.corr(y_real_07_series_train), 5)\n",
    "                            error_val_train= compute_mae_mse_rmse(x_prediction_07_train[:,0],y_real_07_train)\n",
    "                            print('TEST SET scatter corr',corr_ann,'scatter error',error_val,'TEST R2',error_val[3])\n",
    "                            print('TRAINING SET scatter corr',corr_ann_train,'scatter error',error_val_train,'R2',error_val_train[3])\n",
    "                            print(neurons1,epochs_number,learning_rate_search,dropout_rate,batch_size_number,reg,activation1)                          \n",
    "                            x_y_x=np.arange(0,1100,0.1)\n",
    "                            x_y_y=np.arange(0,1100,0.1)\n",
    "                            fig = plt.figure()\n",
    "                            ax = fig.add_subplot(111)\n",
    "                            ax.scatter(x_prediction_07[:,0],y_real_07,color='red',label='Artificial Neural Network Test Set',alpha=0.75)\n",
    "                            ax.scatter(x_prediction_07_train[:,0],y_real_07_train,color='blue',label='Artificial Neural Network Training Set',alpha=0.25,marker=\"^\")\n",
    "                            ax.plot(x_y_x,x_y_y)\n",
    "#                                 ax.axis(([-5,40,-5,40]))\n",
    "                            plt.legend()\n",
    "                            plt.xlabel(u\"Predicted_Power_Density@0.6V (mW cm-2)\")\n",
    "                            plt.ylabel(u\"Real_Power_Density@0.6V (mW cm-2)\")\n",
    "                            fig.savefig('sel PD 0.96 %s %s %s %s %s %s %s.png' %(neurons1,epochs_number,learning_rate_search,dropout_rate,batch_size_number,reg,activation1))\n",
    "                            K.clear_session()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import eli5\n",
    "def base_model():\n",
    "    tmodel = Sequential() \n",
    "    tmodel.add(Dense(800, input_dim=13, kernel_initializer='random_normal',\n",
    "                    bias_initializer='random_normal',activation='relu',kernel_regularizer=keras.regularizers.l2(0.0001))) \n",
    "    tmodel.add(Dropout(0))\n",
    "    tmodel.add(Dense(800, input_dim=800, kernel_initializer='random_normal',\n",
    "                bias_initializer='random_normal',activation='relu',kernel_regularizer=keras.regularizers.l2(0.0001))) \n",
    "    tmodel.add(Dropout(0))\n",
    "    tmodel.add(Dense(1, input_dim=800, activation='linear'))\n",
    "    adam=optimizers.Adam(lr=0.02)\n",
    "    tmodel.compile(loss='mse', optimizer=adam)\n",
    "    return tmodel\n",
    "from eli5.sklearn import PermutationImportance\n",
    "my_model = KerasRegressor(build_fn=base_model,nb_epoch=350, batch_size=32, verbose= False)    \n",
    "my_model.fit(X_train, y_train,validation_split=0.2,callbacks=[early_stopping])\n",
    "perm = PermutationImportance(my_model, random_state=1,n_iter=10).fit(X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eli5.show_weights(perm,top=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smodel = Sequential() \n",
    "smodel.add(Dense(800, input_dim=13, kernel_initializer='random_normal',\n",
    "                bias_initializer='random_normal',activation='relu',kernel_regularizer=keras.regularizers.l2(0.0001))) \n",
    "smodel.add(Dropout(0))\n",
    "smodel.add(Dense(800, input_dim=800, kernel_initializer='random_normal',\n",
    "            bias_initializer='random_normal',activation='relu',kernel_regularizer=keras.regularizers.l2(0.0001))) \n",
    "smodel.add(Dropout(0))\n",
    "smodel.add(Dense(1, input_dim=800, activation='linear'))\n",
    "adam=optimizers.Adam(lr=0.02)\n",
    "smodel.compile(loss='mse', optimizer=adam) \n",
    "print('training...')\n",
    "smodel.fit(X_train, y_train,verbose=0, batch_size=32,epochs=350,validation_split=0.2,callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap\n",
    "# %matplotlib\n",
    "standardized_data_SHAP = standardized_data\n",
    "SHAP_INPUT=standardized_data_SHAP.iloc[:,0:13]\n",
    "SHAP_OUTPUT=standardized_data_SHAP.iloc[:,13]\n",
    "X_SHAP=SHAP_INPUT.values.astype(np.float32)\n",
    "y_SHAP=SHAP_OUTPUT.values.astype(np.float32)\n",
    "explainer = shap.DeepExplainer(smodel,X_SHAP)\n",
    "cols=[ \n",
    "            'Dielectric Constant (C  N-1 M-2)',#0\n",
    "                      'Flash Point(℃)',#1\n",
    "                      'Stirring Method (0 for Mechanical Stirring 1 for Ultrasound)',#2\n",
    "                      'Flow Rate (mL min-1)',#3\n",
    "                      'Water Content (wt%)',#4\n",
    "                      'Anodic Platinum Loading Amount (mgPt cm-2)',#5\n",
    "                      'Cathodic Platinum Loading Amount (mgPt cm-2)',#6\n",
    "                      'I/C',#7\n",
    "                      'Cell Temperature (℃)',#8\n",
    "                      'Anode Flow Rate (sccm)',#9\n",
    "                      'Cathode Flow Rate (sccm)',#10\n",
    "                      'Solid Content (wt%)',#11\n",
    "                      'Backpressure (Mpa)'\n",
    "                        ]\n",
    "shap_values = explainer.shap_values(X_SHAP)\n",
    "print(type(shap_values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib\n",
    "shap.summary_plot(shap_values[0], SHAP_INPUT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.abs(shap_values[0]).mean(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}