{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install --upgrade pip;pip install git+https://github.com/molmod/molmod;pip install pubchempy;pip install auto_ml","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\npd.options.display.max_columns = None\nimport os\nfrom molmod import Molecule\nimport pubchempy as pcp\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor, CatBoostClassifier\nimport keras\n\nfrom auto_ml import Predictor\nfrom auto_ml.utils_models import load_ml_model\n\nimport seaborn as sns\nsns.set()\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport os\nprint(os.listdir(\"../input\"))\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=RuntimeWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\nstructures = pd.read_csv('../input/structures.csv')\nDipoleMovements = pd.read_csv('../input/dipole_moments.csv')\nmagnetic_shielding_tensors = pd.read_csv('../input/magnetic_shielding_tensors.csv')\nmulliken_charges = pd.read_csv('../input/mulliken_charges.csv')\npotential_energy = pd.read_csv('../input/potential_energy.csv')\nscalar_coupling_contributions = pd.read_csv('../input/scalar_coupling_contributions.csv')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Combining tables for training and test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train_comb = train.iloc[:,:-1].append(test, ignore_index=True)\ntrain_comb= train.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape, test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(structures, how='left', left_on=['molecule_name', 'atom_index_0'],\\\n                                        right_on=['molecule_name', 'atom_index'],copy=False, validate=\"m:1\").drop('atom_index',axis=1)\ntrain_comb = train_comb.merge(structures, how='left', left_on=['molecule_name', 'atom_index_1']\\\n                                        ,right_on=['molecule_name', 'atom_index'],  suffixes=('_0', '_1'),copy=False, validate=\"m:1\").drop('atom_index',axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape, test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":false},"cell_type":"code","source":"train_comb = train_comb.merge(DipoleMovements, how='left',on='molecule_name',validate=\"m:1\")\ntrain_comb.rename(columns={'X': 'dipole_moment_X', 'Y': 'dipole_moment_Y','Z': 'dipole_moment_Z'}, inplace=True)\ntrain_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape, test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(magnetic_shielding_tensors, how='left', left_on=['molecule_name', 'atom_index_0'],\\\n                                        right_on=['molecule_name', 'atom_index'],copy=False, validate=\"m:1\").drop('atom_index',axis=1)\ntrain_comb = train_comb.merge(magnetic_shielding_tensors, how='left', left_on=['molecule_name', 'atom_index_1']\\\n                                        ,right_on=['molecule_name', 'atom_index'],  suffixes=('_0', '_1'),copy=False, validate=\"m:1\").drop('atom_index',axis=1)\ntrain_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape , test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(mulliken_charges, how='left', left_on=['molecule_name', 'atom_index_0'],\\\n                                        right_on=['molecule_name', 'atom_index'],copy=False, validate=\"m:1\").drop('atom_index',axis=1)\ntrain_comb = train_comb.merge(mulliken_charges, how='left', left_on=['molecule_name', 'atom_index_1']\\\n                                        ,right_on=['molecule_name', 'atom_index'],  suffixes=('_0', '_1'),copy=False, validate=\"m:1\").drop('atom_index',axis=1)\ntrain_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape , test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(potential_energy, how='left',on='molecule_name',validate=\"m:1\")\ntrain_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape , test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(scalar_coupling_contributions, how='left', left_on=['molecule_name', 'atom_index_0','atom_index_1'],\\\n                                        right_on=['molecule_name', 'atom_index_0','atom_index_1'],copy=False, validate=\"1:1\").drop('type_x',axis=1)\ntrain_comb.rename(columns={'type_y': 'type'},inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape , test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb['#Atoms_in_Molecule'] = train_comb.groupby(['molecule_name'])['id'].transform('count')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Extracting molecular properties as features using molmod (Feature Engineering!)"},{"metadata":{"trusted":true},"cell_type":"code","source":"xyzFiles = os.listdir('../input/structures/')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Molnames = [name[:-4] for name in xyzFiles]\nmolecule_formula_names = []\n\nfor name in xyzFiles:\n    \n    mol = Molecule.from_file(\"../input/structures/\"+name)\n    molecule_formula_names.append(mol.chemical_formula)\n    \nMolnames_df = pd.DataFrame({'Name' : Molnames, 'Chem_formula' : molecule_formula_names})    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Molnames_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb = train_comb.merge(Molnames_df, how='left',left_on='molecule_name',right_on='Name',validate=\"m:1\")\ntrain_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.drop(['id','molecule_name','Name'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.shape, train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.drop(['atom_index_0','atom_index_1'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_comb.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train_comb.iloc[:,1:]\ny = train_comb.iloc[:,0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"list(X.dtypes[X.columns != 'Chem_formula'][X.dtypes == 'object'].index.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_names = list(X.dtypes[X.columns != 'Chem_formula'][X.dtypes == 'object'].index.values)\nX = pd.get_dummies(data=X,columns=cat_names, prefix_sep='_')\nX.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import LabelEncoder\ncats_to_encode = X.dtypes[X.dtypes == 'object'].index.values\nle = LabelEncoder()\nX[cats_to_encode] = X[cats_to_encode].apply(lambda col: le.fit_transform(col))\n\nX.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Basic LGBM"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_test, y_train, y_test = train_test_split(\n    X, y, test_size=0.2, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport time\nimport datetime\nimport json\nimport gc\nfrom numba import jit\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm_notebook\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor, CatBoostClassifier\nfrom sklearn import metrics\n\nfrom itertools import product\n\n\ndef eval_auc(y_true, y_pred):\n    \"\"\"\n    Fast auc eval function for lgb.\n    \"\"\"\n    return 'auc', fast_auc(y_true, y_pred), True\n\n\ndef group_mean_log_mae(y_true, y_pred, types, floor=1e-9):\n    \"\"\"\n    Fast metric computation for this competition: https://www.kaggle.com/c/champs-scalar-coupling\n    Code is from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n    \"\"\"\n    maes = (y_true-y_pred).abs().groupby(types).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()\n\ndef train_model_regression(X, X_test, y, params, folds, model_type='lgb', eval_metric='mae', columns=None, plot_feature_importance=False, model=None,\n                               verbose=10000, early_stopping_rounds=200, n_estimators=50000):\n    \"\"\"\n    A function to train a variety of regression models.\n    Returns dictionary with oof predictions, test predictions, scores and, if necessary, feature importances.\n    \n    :params: X - training data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: X_test - test data, can be pd.DataFrame or np.ndarray (after normalizing)\n    :params: y - target\n    :params: folds - folds to split data\n    :params: model_type - type of model to use\n    :params: eval_metric - metric to use\n    :params: columns - columns to use. If None - use all columns\n    :params: plot_feature_importance - whether to plot feature importance of LGB\n    :params: model - sklearn model, works only for \"sklearn\" model type\n    \n    \"\"\"\n    columns = X.columns if columns is None else columns\n    X_test = X_test[columns]\n    \n    # to set up scoring parameters\n    metrics_dict = {'mae': {'lgb_metric_name': 'mae',\n                        'catboost_metric_name': 'MAE',\n                        'sklearn_scoring_function': metrics.mean_absolute_error},\n                    'group_mae': {'lgb_metric_name': 'mae',\n                        'catboost_metric_name': 'MAE',\n                        'scoring_function': group_mean_log_mae},\n                    'mse': {'lgb_metric_name': 'mse',\n                        'catboost_metric_name': 'MSE',\n                        'sklearn_scoring_function': metrics.mean_squared_error}\n                    }\n\n    \n    result_dict = {}\n    \n    # out-of-fold predictions on train data\n    oof = np.zeros(len(X))\n    \n    # averaged predictions on train data\n    prediction = np.zeros(len(X_test))\n    \n    # list of scores on folds\n    scores = []\n    feature_importance = pd.DataFrame()\n    \n    # split and train on folds\n    for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):\n        print(f'Fold {fold_n + 1} started at {time.ctime()}')\n        if type(X) == np.ndarray:\n            X_train, X_valid = X[columns][train_index], X[columns][valid_index]\n            y_train, y_valid = y[train_index], y[valid_index]\n        else:\n            X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]\n            y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]\n            \n        if model_type == 'lgb':\n            model = lgb.LGBMRegressor(**params, n_estimators = n_estimators, n_jobs = -1)\n            model.fit(X_train, y_train, \n                    eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],\n                    verbose=verbose, early_stopping_rounds=early_stopping_rounds)\n            \n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test, num_iteration=model.best_iteration_)\n            \n        if model_type == 'xgb':\n            train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)\n            valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)\n\n            watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]\n            model = xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, verbose_eval=verbose, params=params)\n            y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n            y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)\n        \n        if model_type == 'sklearn':\n            model = model\n            model.fit(X_train, y_train)\n            \n            y_pred_valid = model.predict(X_valid).reshape(-1,)\n            score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)\n            print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')\n            print('')\n            \n            y_pred = model.predict(X_test).reshape(-1,)\n        \n        if model_type == 'cat':\n            model = CatBoostRegressor(iterations=20000,  eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,\n                                      loss_function=metrics_dict[eval_metric]['catboost_metric_name'])\n            model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)\n\n            y_pred_valid = model.predict(X_valid)\n            y_pred = model.predict(X_test)\n        \n        oof[valid_index] = y_pred_valid.reshape(-1,)\n        if eval_metric != 'group_mae':\n            scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid))\n        else:\n            scores.append(metrics_dict[eval_metric]['scoring_function'](y_valid, y_pred_valid, X_valid['type']))\n\n        prediction += y_pred    \n        \n        if model_type == 'lgb' and plot_feature_importance:\n            # feature importance\n            fold_importance = pd.DataFrame()\n            fold_importance[\"feature\"] = columns\n            fold_importance[\"importance\"] = model.feature_importances_\n            fold_importance[\"fold\"] = fold_n + 1\n            feature_importance = pd.concat([feature_importance, fold_importance], axis=0)\n\n    prediction /= folds.n_splits\n    \n    print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))\n    \n    result_dict['oof'] = oof\n    result_dict['prediction'] = prediction\n    result_dict['scores'] = scores\n    \n    if model_type == 'lgb':\n        if plot_feature_importance:\n            feature_importance[\"importance\"] /= folds.n_splits\n            cols = feature_importance[[\"feature\", \"importance\"]].groupby(\"feature\").mean().sort_values(\n                by=\"importance\", ascending=False)[:50].index\n\n            best_features = feature_importance.loc[feature_importance.feature.isin(cols)]\n\n            plt.figure(figsize=(16, 12));\n            sns.barplot(x=\"importance\", y=\"feature\", data=best_features.sort_values(by=\"importance\", ascending=False));\n            plt.title('LGB Features (avg over folds)');\n            \n            result_dict['feature_importance'] = feature_importance\n        \n    return result_dict","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'objective': 'huber',\n          'max_depth': -1,\n          'learning_rate': 0.01,\n          \"boosting\": \"gbdt\",\n          \"bagging_freq\": 5,\n          \"bagging_seed\": 11,\n          \"metric\": 'mae',\n          \"verbosity\": -1,\n         }\n\n#test our model framework on first 50000 training samo=ples first and first 10000 test samples\n\n\noutput = train_model_regression(X_train[:50000], X_test[:10000], y_train[:50000], params, folds, model_type='lgb', eval_metric='mae', columns=None, plot_feature_importance=False, model=None,\n                               verbose=10000, early_stopping_rounds=200, n_estimators=50000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# print(metrics.mean_absolute_error(y_test[:10000], prediction_lgb[:10000]))\noutput","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Validation MAE for LGBM = {metrics.mean_absolute_error(y_test[:10000], output[\"prediction\"][:10000])}')\nprint(f'Validation MAE for just predicting the mean of test set = {metrics.mean_absolute_error(y_test[:10000], y_test[:10000].mean()*np.ones(y_test[:10000].shape))}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Now we need to find a way to get these features into the test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}