{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport xgboost\nimport time\nimport gc\nimport shap\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dtypes = {'atom_index_0':'uint8', \n          'atom_index_1':'uint8', \n          'scalar_coupling_constant':'float32', \n          'num_C':'uint8', \n          'num_H':'uint8', \n          'num_N':'uint8', \n          'num_O':'uint8', \n          'num_F':'uint8',\n          'total_atoms':'uint8',\n          'num_bonds':'uint8', \n          'num_mol_bonds':'uint8', \n          'min_d':'float32', \n          'mean_d':'float32', \n          'max_d':'float32', \n          'space_dr':'float32', \n          'bond_dr':'float32',\n          'bond_1':'uint8', \n          'bond_2':'uint8', \n          'bond_3':'uint8', \n          'atom_0_pc':'float32', \n          'atom_end_pc':'float32',\n          'atom_2_hyb':'uint8', \n          'atom_3_hyb':'uint8', \n          'atom_end_hyb':'uint8', \n          'path_count':'uint8', \n          'atom_0_min':'float32',\n          'atom_0_mean':'float32', \n          'atom_0_max':'float32', \n          'atom_0_Cmin':'float32', \n          'atom_0_Cmean':'float32',\n          'atom_0_Cmax':'float32', \n          'atom_0_Omin':'float32', \n          'atom_0_Omean':'float32',\n          'atom_0_Omax':'float32', \n          'atom_0_Nmin':'float32', \n          'atom_0_Nmean':'float32', \n          'atom_0_Nmax':'float32',\n          'atom_0_Fmin':'float32', \n          'atom_0_Fmean':'float32', \n          'atom_0_Fmax':'float32', \n          'atom_end_min':'float32',\n          'atom_end_mean':'float32', \n          'atom_end_max':'float32', \n          'atom_end_Cmin':'float32', \n          'atom_end_Cmean':'float32',\n          'atom_end_Cmax':'float32', \n          'atom_end_Omin':'float32', \n          'atom_end_Omean':'float32',\n          'atom_end_Omax':'float32', \n          'atom_end_Nmin':'float32', \n          'atom_end_Nmean':'float32', \n          'atom_end_Nmax':'float32',\n          'atom_end_Fmin':'float32', \n          'atom_end_Fmean':'float32', \n          'atom_end_Fmax':'float32',\n          'Dmin_COM':'float32', \n          'Dmean_COM':'float32', \n          'Dmax_COM':'float32',\n          'COM_dr_0': 'float32',\n          'COM_dr_1': 'float32',\n          'bond2_angle': 'float32',\n          'bond3_angle': 'float32'\n         }","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/predmolprop-featureengineering-final/train_extend.csv\",dtype=dtypes)\ntest = pd.read_csv(\"../input/predmolprop-featureengineering-finaltest/test_extend.csv\",dtype=dtypes)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"PopList = ['molecule_name', 'atom_index_0','atom_index_1','num_bonds','atom_end_type',\n           'bond_1','is_linear','atom_2_hyb','atom_3_hyb','atom_end_hyb',\n           'atom_0_fc', 'atom_end_fc','atom_0_val','atom_end_val','atom_0_sm', 'atom_end_sm']\n\nfor col in PopList:\n    if col in train:\n        train.pop(col)\n        test.pop(col)\n    \ntrain.fillna(value ='',inplace= True)\ntest.fillna(value='',inplace=True)\n\ncoupling_types = sorted(list(train.type.unique()))\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Encode categorical features\nfrom sklearn.preprocessing import LabelEncoder\n\ncols = ['atom_0_type2','atom_2_type','atom_3_type','atom_end_type2']\nfor col in cols:\n    enc = LabelEncoder()\n    train[col]=enc.fit_transform(train[col]).astype(np.uint8)\n    test[col]=enc.transform(test[col]).astype(np.uint8)\ndel cols","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[train.type=='1JHC'].loc[:,train.columns.map(lambda x: x in set(['space_dr','bond_dr']))].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# feature lists\nprefix_train= ['id', 'type', 'scalar_coupling_constant']\nprefix_test= ['id', 'type']\n\nfc_distance = ['space_dr','min_d','mean_d', 'max_d']\n\nfc_COM = ['Dmin_COM', 'Dmean_COM', 'Dmax_COM']\n\nfc_size = ['num_mol_bonds', 'total_atoms','num_C', 'num_H', 'num_N', 'num_O', 'num_F']\n\nfc_atom_0 = ['atom_0_pc','atom_0_type2','COM_dr_0',\n             'atom_0_min','atom_0_mean', 'atom_0_max',\n             'atom_0_Cmin', 'atom_0_Cmean','atom_0_Cmax',\n             'atom_0_Nmin', 'atom_0_Nmean', 'atom_0_Nmax',\n             'atom_0_Omin', 'atom_0_Omean','atom_0_Omax',\n             'atom_0_Fmin', 'atom_0_Fmean', 'atom_0_Fmax']\n\nfc_atom_end = ['atom_end_pc', 'atom_end_type2', 'COM_dr_1',\n               'atom_end_min','atom_end_mean', 'atom_end_max',\n               'atom_end_Cmin', 'atom_end_Cmean','atom_end_Cmax',\n               'atom_end_Nmin', 'atom_end_Nmean', 'atom_end_Nmax',\n               'atom_end_Omin', 'atom_end_Omean','atom_end_Omax',\n               'atom_end_Fmin', 'atom_end_Fmean', 'atom_end_Fmax']\n\nfc_2 = ['bond_dr','path_count','atom_2_type','bond_2','bond2_angle']\nfc_3 = ['atom_3_type','bond_3','bond3_angle']\n\nfc1 = set(fc_distance+fc_COM+fc_size+fc_atom_0+fc_atom_end)\n\n# reorder\ntrain = train[prefix_train+fc_distance+fc_COM+fc_size+fc_atom_0+fc_atom_end+fc_2+fc_3]\ntest = test[prefix_test+fc_distance+fc_COM+fc_size+fc_atom_0+fc_atom_end+fc_2+fc_3]\nlen(train.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Process Data\nfrom sklearn.model_selection import train_test_split\ndef ProcessData(df,features,test_size=0.25):\n    if test_size == 0:\n        train_Y = df.pop('scalar_coupling_constant')\n        train_type = df.pop('type')\n        df.pop('id')\n        return df.loc[:,df.columns.map(lambda x: x in features)], train_Y, train_type\n    \n    train_X, val_X, train_Y, val_Y = train_test_split(df.loc[:,df.columns.map(lambda x: x in features or x=='type')], \n                                                      df.scalar_coupling_constant,test_size=test_size,random_state=42)\n    \n    train_type = train_X.pop('type')\n    val_type = val_X.pop('type')\n    return train_X, train_Y, train_type, val_X, val_Y, val_type","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"[https://www.kaggle.com/uberkinder/efficient-metric](http://)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def CalcLMAE(y_true, y_pred, groups, floor=1e-9):\n    maes = (y_true-y_pred).abs().groupby(groups).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from xgboost import XGBRegressor\nfrom xgboost import XGBRFRegressor\nfrom math import log\n\ndef SingleRun(df,features, test_size=0.25, model_fn=XGBRegressor, includeType=False, early_stopping_rounds=None, do_SHAP=False, **kwargs):\n    data = ProcessData(df,features,test_size)\n    if(test_size==0):\n        train_X,train_Y,train_type = data\n    else:\n        train_X,train_Y,train_type,val_X,val_Y,val_type = data\n    if includeType:\n        train_X=train_X.join(train_type)\n        val_X=val_X.join(val_type)\n        enc = LabelEncoder()\n        train_X.type=enc.fit_transform(train_X.type).astype(np.uint8)\n        val_X.type=enc.transform(val_X.type).astype(np.uint8)\n    \n    model = model_fn(**kwargs)\n    t1=time.time()\n    if early_stopping_rounds is None:\n        model.fit(train_X,train_Y, verbose=False)\n    else:\n        if test_size==0:\n            print('ERROR: need test data for early_stopping_rounds')\n            return -1,-1,-1\n        model.fit(train_X,train_Y, early_stopping_rounds=5, eval_set=[(val_X, val_Y)], verbose=False)\n        print('best LMAE:',log(model.best_score))\n        print('best ntree:', model.best_ntree_limit)\n    t2=time.time()\n    print('training time:',t2-t1)\n    train_predict = pd.Series(model.predict(train_X),index=train_X.index)\n    train_LMAE = CalcLMAE(train_Y,train_predict,train_type)\n    print('\\ttrain LMAE:',train_LMAE)\n    if(test_size==0):\n        val_LMAE = None\n#        g = sns.FacetGrid(pd.DataFrame({'type':train_type,'scalar_coupling_constant': train_Y,'predictions':train_predict}), \n#                          col=\"type\", col_order = coupling_types,sharex=False,sharey=False)\n#        g.map(sns.scatterplot, \"scalar_coupling_constant\",\"predictions\")\n        plt.figure(figsize=[6,6])\n        sns.scatterplot(x=train_Y, y=train_predict)\n        plt.plot(train_Y,train_Y,color='black')\n        plt.title(train_type.values[0])\n        plt.xlabel('scalar coupling constant')\n        plt.ylabel('predicted value')\n        plt.show()\n        \n        if do_SHAP:\n            samp = train_X.sample(200)\n            shap.summary_plot(model.get_booster().predict(xgboost.DMatrix(samp),pred_contribs=True)[:,:-1], samp)#, max_display=len(features))\n        #explainer = shap.TreeExplainer(model)\n        #shap.summary_plot(explainer.shap_values(samp), samp, max_display=len(features))\n    \n    else:\n        val_predict = pd.Series(model.predict(val_X),index=val_X.index)\n        val_LMAE = CalcLMAE(val_Y,val_predict,val_type)\n        print('\\tval LMAE:',val_LMAE)\n#        g = sns.FacetGrid(pd.DataFrame({'type':val_type,'scalar_coupling_constant': val_Y,'predictions':val_predict}), \n#                          col=\"type\", col_order = coupling_types,sharex=False,sharey=False)\n#        g.map(sns.scatterplot, \"scalar_coupling_constant\",\"predictions\")\n        plt.figure(figsize=[6,6])\n        sns.scatterplot(x=val_Y, y=val_predict)\n        plt.plot(val_Y,val_Y,color='black')\n        plt.title(train_type.values[0])\n        plt.xlabel('scalar coupling constant')\n        plt.ylabel('predicted value')\n        plt.show()\n        \n        if do_SHAP:\n            samp = val_X.sample(200)\n            shap.summary_plot(model.get_booster().predict(xgboost.DMatrix(samp),pred_contribs=True)[:,:-1], samp)#, max_display=len(features))\n        #explainer = shap.TreeExplainer(model)\n        #shap.summary_plot(explainer.shap_values(samp), samp, max_display=len(features))\n    \n    gc.collect()\n    return model, train_LMAE, val_LMAE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coupling_type = '3JHN'\ntrain_sample = train[train.type==coupling_type]\n\nfeatures=fc1.copy()\nif(coupling_type[0]=='2'):\n    features.update(fc_2)\nelif(coupling_type[0]=='3'):\n    features.update(fc_2[:-1]+fc_3)\n\nmodel_3JHN,_,_=SingleRun(train_sample,features,test_size=0.2,model_fn=XGBRegressor,includeType=False,early_stopping_rounds=5, \n                         do_SHAP = True,\n                         max_depth=11, learning_rate=0.1, n_estimators=3000, \n                         verbosity=1, \n                         objective='reg:squarederror', booster='gbtree',tree_method= 'gpu_hist',\n                         n_jobs=4, \n                         gamma=0, min_child_weight=1, max_delta_step=0, \n                         subsample=1,colsample_bytree=1, colsample_bylevel=1, colsample_bynode=1, \n                         reg_alpha=0, reg_lambda=1, scale_pos_weight=1, base_score=0.5, \n                         random_state=0, seed=None, missing=None, importance_type='gain', eval_metric='mae')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = xgboost.plot_tree(model_3JHN,num_trees=0)\nfig = ax.figure\nfig.set_size_inches(100,100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax=xgboost.plot_importance(model_3JHN,importance_type='gain',show_values=False,log=True)\nfig = ax.figure\nfig.set_size_inches(5, 10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax=xgboost.plot_importance(model_3JHN,importance_type='weight',show_values=False,log=True)\nfig = ax.figure\nfig.set_size_inches(5, 10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax=xgboost.plot_importance(model_3JHN,importance_type='cover',show_values=False,log=True)\nfig = ax.figure\nfig.set_size_inches(5, 10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def RunByType(df,test_size=0.25,model_fn=XGBRegressor,includeType=False,early_stopping_rounds=None,**kwargs):\n    model_dict={}\n    train_LMAE_dict={}\n    val_LMAE_dict={}\n    for coupling_type in coupling_types:\n        print('Now training type:',str(coupling_type))\n\n        df_type = df[df['type']==coupling_type]\n\n        features=fc1.copy()\n        if(coupling_type[0]=='2'):\n            features.update(fc_2)\n        elif(coupling_type[0]=='3'):\n            features.update(fc_2[:-1]+fc_3)\n        \n        model,train_LMAE,val_LMAE = SingleRun(df_type,features,test_size,model_fn=model_fn,includeType=False,early_stopping_rounds=early_stopping_rounds,**kwargs)\n\n        model_dict[coupling_type]=model\n        train_LMAE_dict[coupling_type]=train_LMAE\n        val_LMAE_dict[coupling_type]=val_LMAE\n        \n    train_LMAE = sum(train_LMAE_dict.values())/len(coupling_types)\n    if test_size>0:\n        val_LMAE = sum(val_LMAE_dict.values())/len(coupling_types)\n    else:\n        val_LMAE=None\n    print('total train LMAE:',train_LMAE)\n    print('total val LMAE', val_LMAE)\n    gc.collect()\n    return model_dict#,train_LMAE, val_LMAE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def PredictByType(df_X,model_dict,df_Y=None):\n    predictions = pd.DataFrame()\n\n    for coupling_type in coupling_types:\n        print('predicting type:',str(coupling_type))\n        model = model_dict[coupling_type]\n        df_type = df_X[df_X['type']==coupling_type]\n        \n        features=fc1.copy()\n        if(coupling_type[0]=='2'):\n            features.update(fc_2)\n        elif(coupling_type[0]=='3'):\n            features.update(fc_2[:-1]+fc_3)\n        \n        predict = model.predict(df_type.loc[:,df_X.columns.map(lambda col: col in features)])\n        predict = pd.Series(predict,index=df_type.index)\n        predict = pd.DataFrame({'id':df_type.id,'scalar_coupling_constant': predict})\n        predict.set_index('id',inplace=True)\n        predictions = pd.concat([predictions,predict])\n    \n    if df_Y is not None:\n        predictions.columns = ['predictions']\n        merged = pd.merge(train,predictions,how='left',left_on='id',right_index=True)\n        LMAE = CalcLMAE(df_Y, merged.predictions,merged.type)\n        print('LMAE:', LMAE)\n        g = sns.FacetGrid(merged,col=\"type\", col_order = coupling_types,sharex=False,sharey=False)\n        g.map(sns.scatterplot, \"scalar_coupling_constant\",\"predictions\")\n        del merged\n    \n    gc.collect()\n    return predictions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_dict=RunByType(train,test_size=0.2,includeType=False,early_stopping_rounds=5,\n                max_depth=11, learning_rate=0.1, n_estimators=10000, \n                verbosity=1, \n                objective='reg:squarederror', booster='gbtree',tree_method= 'gpu_hist',\n                n_jobs=4, \n                gamma=0, min_child_weight=1, max_delta_step=0, \n                subsample=1,colsample_bytree=1, colsample_bylevel=1, colsample_bynode=1, \n                reg_alpha=0, reg_lambda=1, scale_pos_weight=1, base_score=0.5, \n                random_state=0, seed=None, missing=None, importance_type='gain', eval_metric='mae')\n\n#train_predict = PredictByType(val_sample,model_dict,val_sample.scalar_coupling_constant)\ntest_predict = PredictByType(test,model_dict)\ntest_predict.to_csv('test_predict_TypeModel.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('FINISHED!')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}