{"cells":[{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\n\n

Can you measure the magnetic interactions between a pair of atoms?

\n\nThis kernel is a combination of multiple kernels. The goal is to organize and explain the code to beginner competitors like me.
\nThis Kernels creates lots of new features and uses lightgbm as model"},{"metadata":{},"cell_type":"markdown","source":"# Table of Contents:\n\n**1. [Problem Definition](#id1)**
\n**2. [Get the Data (Collect / Obtain)](#id2)**
\n**3. [Load the Dataset](#id3)**
\n**4. [Data Pre-processing](#id4)**
\n**5. [Model](#id5)**
\n**6. [Visualization and Analysis of Results](#id6)**
\n**7. [Submittion](#id7)**
\n**8. [References](#ref)**
"},{"metadata":{},"cell_type":"markdown","source":"
\n# **1. Problem Definition:** \n\nThis challenge aims to predict interactions between atoms. The main task is develop an algorithm that can predict the magnetic interaction between two atoms in a molecule (i.e., the scalar coupling constant)
\n\nIn this competition, you will be predicting the scalar_coupling_constant between atom pairs in molecules, given the two atom types (e.g., C and H), the coupling type (e.g., 2JHC), and any features you are able to create from the molecule structure (xyz) files.\n\n**Data**\n* **train.csv** - the training set, where the first column (molecule_name) is the name of the molecule where the coupling constant originates, the second (atom_index_0) and third column (atom_index_1) is the atom indices of the atom-pair creating the coupling and the fourth column (**scalar_coupling_constant**) is the scalar coupling constant that we want to be able to predict\n* **test.csv** - the test set; same info as train, without the target variable\n* **sample_submission.csv** - a sample submission file in the correct format\n* **structures.csv** - this file contains the same information as the individual xyz structure files, but in a single file\n\n**Additional Data**
\n*NOTE: additional data is provided for the molecules in Train only!*\n* **scalar_coupling_contributions.csv** - The scalar coupling constants in train.csv are a sum of four terms. The first column (**molecule_name**) are the name of the molecule, the second (**atom_index_0**) and third column (**atom_index_1**) are the atom indices of the atom-pair, the fourth column indicates the **type** of coupling, the fifth column (**fc**) is the Fermi Contact contribution, the sixth column (**sd**) is the Spin-dipolar contribution, the seventh column (**pso**) is the Paramagnetic spin-orbit contribution and the eighth column (**dso**) is the Diamagnetic spin-orbit contribution.\n\n"},{"metadata":{},"cell_type":"markdown","source":"
\n# **2. Get the Data (Collect / Obtain):** "},{"metadata":{},"cell_type":"markdown","source":"## All imports used in this kernel"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":false},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.metrics import mean_absolute_error\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport xgboost as xgb\nimport time\nimport datetime\nfrom catboost import CatBoostRegressor\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold\nfrom sklearn import metrics\nfrom sklearn import linear_model\nimport gc\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom IPython.display import HTML\nimport json\nimport altair as alt\n\nimport networkx as nx\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport 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\n\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom catboost import CatBoostRegressor, CatBoostClassifier\nfrom sklearn import metrics\n\nfrom itertools import product\n\nimport altair as alt\nfrom altair.vega import v3\nfrom IPython.display import HTML\nalt.renderers.enable('notebook')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## All function used in this kernel"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_kg_hide-input":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"# using ideas from this kernel: https://www.kaggle.com/notslush/altair-visualization-2018-stackoverflow-survey\ndef prepare_altair():\n \"\"\"\n Helper function to prepare altair for working.\n \"\"\"\n vega_url = 'https://cdn.jsdelivr.net/npm/vega@' + v3.SCHEMA_VERSION\n vega_lib_url = 'https://cdn.jsdelivr.net/npm/vega-lib'\n vega_lite_url = 'https://cdn.jsdelivr.net/npm/vega-lite@' + alt.SCHEMA_VERSION\n vega_embed_url = 'https://cdn.jsdelivr.net/npm/vega-embed@3'\n noext = \"?noext\"\n \n paths = {\n 'vega': vega_url + noext,\n 'vega-lib': vega_lib_url + noext,\n 'vega-lite': vega_lite_url + noext,\n 'vega-embed': vega_embed_url + noext\n }\n \n workaround = f\"\"\" requirejs.config({{\n baseUrl: 'https://cdn.jsdelivr.net/npm/',\n paths: {paths}\n }});\n \"\"\"\n \n return workaround\n \n\ndef add_autoincrement(render_func):\n # Keep track of unique
IDs\n cache = {}\n def wrapped(chart, id=\"vega-chart\", autoincrement=True):\n if autoincrement:\n if id in cache:\n counter = 1 + cache[id]\n cache[id] = counter\n else:\n cache[id] = 0\n actual_id = id if cache[id] == 0 else id + '-' + str(cache[id])\n else:\n if id not in cache:\n cache[id] = 0\n actual_id = id\n return render_func(chart, id=actual_id)\n # Cache will stay outside and \n return wrapped\n \n\n@add_autoincrement\ndef render(chart, id=\"vega-chart\"):\n \"\"\"\n Helper function to plot altair visualizations.\n \"\"\"\n chart_str = \"\"\"\n
\n \"\"\"\n return HTML(\n chart_str.format(\n id=id,\n chart=json.dumps(chart) if isinstance(chart, dict) else chart.to_json(indent=None)\n )\n )\n \n\ndef reduce_mem_usage(df, verbose=True):\n numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n start_mem = df.memory_usage().sum() / 1024**2 \n for col in df.columns:\n col_type = df[col].dtypes\n if col_type in numerics:\n c_min = df[col].min()\n c_max = df[col].max()\n if str(col_type)[:3] == 'int':\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n df[col] = df[col].astype(np.int64) \n else:\n if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n df[col] = df[col].astype(np.float16)\n elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float32)\n else:\n df[col] = df[col].astype(np.float64) \n end_mem = df.memory_usage().sum() / 1024**2\n if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n return df\n \n\n@jit\ndef fast_auc(y_true, y_prob):\n \"\"\"\n fast roc_auc computation: https://www.kaggle.com/c/microsoft-malware-prediction/discussion/76013\n \"\"\"\n y_true = np.asarray(y_true)\n y_true = y_true[np.argsort(y_prob)]\n nfalse = 0\n auc = 0\n n = len(y_true)\n for i in range(n):\n y_i = y_true[i]\n nfalse += (1 - y_i)\n auc += y_i * nfalse\n auc /= (nfalse * (n - nfalse))\n return auc\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 \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\n \n\n\ndef train_model_classification(X, X_test, y, params, folds, model_type='lgb', eval_metric='auc', 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 == None else columns\n X_test = X_test[columns]\n \n # to set up scoring parameters\n metrics_dict = {'auc': {'lgb_metric_name': eval_auc,\n 'catboost_metric_name': 'AUC',\n 'sklearn_scoring_function': metrics.roc_auc_score},\n }\n \n result_dict = {}\n \n # out-of-fold predictions on train data\n oof = np.zeros((len(X), len(set(y.values))))\n \n # averaged predictions on train data\n prediction = np.zeros((len(X_test), oof.shape[1]))\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.LGBMClassifier(**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_proba(X_valid)\n y_pred = model.predict_proba(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=n_estimators, evals=watchlist, early_stopping_rounds=early_stopping_rounds, 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_proba(X_test)\n \n if model_type == 'cat':\n model = CatBoostClassifier(iterations=n_estimators, 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\n scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid[:, 1]))\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\n\n# setting up altair\nworkaround = prepare_altair()\nHTML(\"\".join((\n \"\",\n)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n# **3. Load the Dataset** \n\nLet's load all necessary datasets"},{"metadata":{"trusted":false},"cell_type":"code","source":"\ntrain = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\nsub = pd.read_csv('../input/sample_submission.csv')\nstructures = pd.read_csv('../input/structures.csv')\nscalar_coupling_contributions = pd.read_csv('../input/scalar_coupling_contributions.csv')\n\n#Not used by now!\n#magnetic_shielding_tensors = pd.read_csv('../input/magnetic_shielding_tensors.csv')\n#dipole_moments = pd.read_csv('../input/dipole_moments.csv')\n#potential_energy = pd.read_csv('../input/potential_energy.csv')\n#mulliken_charges = pd.read_csv('../input/mulliken_charges.csv')\n\nprint('Train dataset shape is -> rows: {} cols:{}'.format(train.shape[0],train.shape[1]))\nprint('Test dataset shape is -> rows: {} cols:{}'.format(test.shape[0],test.shape[1]))\nprint('Sub dataset shape is -> rows: {} cols:{}'.format(sub.shape[0],sub.shape[1]))\nprint('Structures dataset shape is -> rows: {} cols:{}'.format(structures.shape[0],structures.shape[1]))\nprint('Scalar_coupling_contributions dataset shape is -> rows: {} cols:{}'.format(scalar_coupling_contributions.shape[0],\n scalar_coupling_contributions.shape[1]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For an fast model/feature evaluation, get only 10% of dataset. Final submission must remove/coments this code"},{"metadata":{"trusted":false},"cell_type":"code","source":"'''\nsize = round(0.10*train.shape[0])\ntrain = train[:size]\ntest = test[:size]\nsub = sub[:size]\nstructures = structures[:size]\nscalar_coupling_contributions = scalar_coupling_contributions[:size]\n\nprint('Train dataset shape is now rows: {} cols:{}'.format(train.shape[0],train.shape[1]))\nprint('Test dataset shape is now rows: {} cols:{}'.format(test.shape[0],test.shape[1]))\nprint('Sub dataset shape is now rows: {} cols:{}'.format(sub.shape[0],sub.shape[1]))\nprint('Structures dataset shape is now rows: {} cols:{}'.format(structures.shape[0],structures.shape[1]))\nprint('Scalar_coupling_contributions dataset shape is now rows: {} cols:{}'.format(scalar_coupling_contributions.shape[0],\n scalar_coupling_contributions.shape[1]))\n'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The importante things to know is that the scalar coupling constants in train.csv are a sum of four terms. \n```\n* fc is the Fermi Contact contribution\n* sd is the Spin-dipolar contribution\n* pso is the Paramagnetic spin-orbit contribution\n* dso is the Diamagnetic spin-orbit contribution. \n```\nLet's merge this into train"},{"metadata":{"trusted":false},"cell_type":"code","source":"train = pd.merge(train, scalar_coupling_contributions, how = 'left',\n left_on = ['molecule_name', 'atom_index_0', 'atom_index_1', 'type'],\n right_on = ['molecule_name', 'atom_index_0', 'atom_index_1', 'type'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"train.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"test.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Based in others ideais we can:
\n\n- train a model to predict `fc` feature;\n- add this feature to train and test and train the same model to compare performance;\n- train a better model;"},{"metadata":{},"cell_type":"markdown","source":"
\n# **4. Data Pre-processing** "},{"metadata":{},"cell_type":"markdown","source":"## Feature generation"},{"metadata":{"_kg_hide-input":true,"trusted":false},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n df = pd.merge(df, structures, how = 'left',\n left_on = ['molecule_name', f'atom_index_{atom_idx}'],\n right_on = ['molecule_name', 'atom_index'])\n \n df = df.drop('atom_index', axis=1)\n df = df.rename(columns={'atom': f'atom_{atom_idx}',\n 'x': f'x_{atom_idx}',\n 'y': f'y_{atom_idx}',\n 'z': f'z_{atom_idx}'})\n return df\n\ntrain = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntest = map_atom_info(test, 0)\ntest = map_atom_info(test, 1)\n\ntrain_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ntest['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\ntrain['dist_x'] = (train['x_0'] - train['x_1']) ** 2\ntest['dist_x'] = (test['x_0'] - test['x_1']) ** 2\ntrain['dist_y'] = (train['y_0'] - train['y_1']) ** 2\ntest['dist_y'] = (test['y_0'] - test['y_1']) ** 2\ntrain['dist_z'] = (train['z_0'] - train['z_1']) ** 2\ntest['dist_z'] = (test['z_0'] - test['z_1']) ** 2\n\ntrain['type_0'] = train['type'].apply(lambda x: x[0])\ntest['type_0'] = test['type'].apply(lambda x: x[0])\n\ndef create_features(df):\n df['molecule_couples'] = df.groupby('molecule_name')['id'].transform('count')\n df['molecule_dist_mean'] = df.groupby('molecule_name')['dist'].transform('mean')\n df['molecule_dist_min'] = df.groupby('molecule_name')['dist'].transform('min')\n df['molecule_dist_max'] = df.groupby('molecule_name')['dist'].transform('max')\n df['atom_0_couples_count'] = df.groupby(['molecule_name', 'atom_index_0'])['id'].transform('count')\n df['atom_1_couples_count'] = df.groupby(['molecule_name', 'atom_index_1'])['id'].transform('count')\n df[f'molecule_atom_index_0_x_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['x_1'].transform('std')\n df[f'molecule_atom_index_0_y_1_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('mean')\n df[f'molecule_atom_index_0_y_1_mean_diff'] = df[f'molecule_atom_index_0_y_1_mean'] - df['y_1']\n df[f'molecule_atom_index_0_y_1_mean_div'] = df[f'molecule_atom_index_0_y_1_mean'] / df['y_1']\n df[f'molecule_atom_index_0_y_1_max'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('max')\n df[f'molecule_atom_index_0_y_1_max_diff'] = df[f'molecule_atom_index_0_y_1_max'] - df['y_1']\n df[f'molecule_atom_index_0_y_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('std')\n df[f'molecule_atom_index_0_z_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['z_1'].transform('std')\n df[f'molecule_atom_index_0_dist_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('mean')\n df[f'molecule_atom_index_0_dist_mean_diff'] = df[f'molecule_atom_index_0_dist_mean'] - df['dist']\n df[f'molecule_atom_index_0_dist_mean_div'] = df[f'molecule_atom_index_0_dist_mean'] / df['dist']\n df[f'molecule_atom_index_0_dist_max'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('max')\n df[f'molecule_atom_index_0_dist_max_diff'] = df[f'molecule_atom_index_0_dist_max'] - df['dist']\n df[f'molecule_atom_index_0_dist_max_div'] = df[f'molecule_atom_index_0_dist_max'] / df['dist']\n df[f'molecule_atom_index_0_dist_min'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('min')\n df[f'molecule_atom_index_0_dist_min_diff'] = df[f'molecule_atom_index_0_dist_min'] - df['dist']\n df[f'molecule_atom_index_0_dist_min_div'] = df[f'molecule_atom_index_0_dist_min'] / df['dist']\n df[f'molecule_atom_index_0_dist_std'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('std')\n df[f'molecule_atom_index_0_dist_std_diff'] = df[f'molecule_atom_index_0_dist_std'] - df['dist']\n df[f'molecule_atom_index_0_dist_std_div'] = df[f'molecule_atom_index_0_dist_std'] / df['dist']\n df[f'molecule_atom_index_1_dist_mean'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('mean')\n df[f'molecule_atom_index_1_dist_mean_diff'] = df[f'molecule_atom_index_1_dist_mean'] - df['dist']\n df[f'molecule_atom_index_1_dist_mean_div'] = df[f'molecule_atom_index_1_dist_mean'] / df['dist']\n df[f'molecule_atom_index_1_dist_max'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('max')\n df[f'molecule_atom_index_1_dist_max_diff'] = df[f'molecule_atom_index_1_dist_max'] - df['dist']\n df[f'molecule_atom_index_1_dist_max_div'] = df[f'molecule_atom_index_1_dist_max'] / df['dist']\n df[f'molecule_atom_index_1_dist_min'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('min')\n df[f'molecule_atom_index_1_dist_min_diff'] = df[f'molecule_atom_index_1_dist_min'] - df['dist']\n df[f'molecule_atom_index_1_dist_min_div'] = df[f'molecule_atom_index_1_dist_min'] / df['dist']\n df[f'molecule_atom_index_1_dist_std'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('std')\n df[f'molecule_atom_index_1_dist_std_diff'] = df[f'molecule_atom_index_1_dist_std'] - df['dist']\n df[f'molecule_atom_index_1_dist_std_div'] = df[f'molecule_atom_index_1_dist_std'] / df['dist']\n df[f'molecule_atom_1_dist_mean'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('mean')\n df[f'molecule_atom_1_dist_min'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('min')\n df[f'molecule_atom_1_dist_min_diff'] = df[f'molecule_atom_1_dist_min'] - df['dist']\n df[f'molecule_atom_1_dist_min_div'] = df[f'molecule_atom_1_dist_min'] / df['dist']\n df[f'molecule_atom_1_dist_std'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('std')\n df[f'molecule_atom_1_dist_std_diff'] = df[f'molecule_atom_1_dist_std'] - df['dist']\n df[f'molecule_type_0_dist_std'] = df.groupby(['molecule_name', 'type_0'])['dist'].transform('std')\n df[f'molecule_type_0_dist_std_diff'] = df[f'molecule_type_0_dist_std'] - df['dist']\n df[f'molecule_type_dist_mean'] = df.groupby(['molecule_name', 'type'])['dist'].transform('mean')\n df[f'molecule_type_dist_mean_diff'] = df[f'molecule_type_dist_mean'] - df['dist']\n df[f'molecule_type_dist_mean_div'] = df[f'molecule_type_dist_mean'] / df['dist']\n df[f'molecule_type_dist_max'] = df.groupby(['molecule_name', 'type'])['dist'].transform('max')\n df[f'molecule_type_dist_min'] = df.groupby(['molecule_name', 'type'])['dist'].transform('min')\n df[f'molecule_type_dist_std'] = df.groupby(['molecule_name', 'type'])['dist'].transform('std')\n df[f'molecule_type_dist_std_diff'] = df[f'molecule_type_dist_std'] - df['dist']\n df = reduce_mem_usage(df)\n return df\n\ntrain = create_features(train)\ntest = create_features(test)\n\ngood_columns = [\n'molecule_atom_index_0_dist_min',\n'molecule_atom_index_0_dist_max',\n'molecule_atom_index_1_dist_min',\n'molecule_atom_index_0_dist_mean',\n'molecule_atom_index_0_dist_std',\n'dist',\n'molecule_atom_index_1_dist_std',\n'molecule_atom_index_1_dist_max',\n'molecule_atom_index_1_dist_mean',\n'molecule_atom_index_0_dist_max_diff',\n'molecule_atom_index_0_dist_max_div',\n'molecule_atom_index_0_dist_std_diff',\n'molecule_atom_index_0_dist_std_div',\n'atom_0_couples_count',\n'molecule_atom_index_0_dist_min_div',\n'molecule_atom_index_1_dist_std_diff',\n'molecule_atom_index_0_dist_mean_div',\n'atom_1_couples_count',\n'molecule_atom_index_0_dist_mean_diff',\n'molecule_couples',\n'atom_index_1',\n'molecule_dist_mean',\n'molecule_atom_index_1_dist_max_diff',\n'molecule_atom_index_0_y_1_std',\n'molecule_atom_index_1_dist_mean_diff',\n'molecule_atom_index_1_dist_std_div',\n'molecule_atom_index_1_dist_mean_div',\n'molecule_atom_index_1_dist_min_diff',\n'molecule_atom_index_1_dist_min_div',\n'molecule_atom_index_1_dist_max_div',\n'molecule_atom_index_0_z_1_std',\n'y_0',\n'molecule_type_dist_std_diff',\n'molecule_atom_1_dist_min_diff',\n'molecule_atom_index_0_x_1_std',\n'molecule_dist_min',\n'molecule_atom_index_0_dist_min_diff',\n'molecule_atom_index_0_y_1_mean_diff',\n'molecule_type_dist_min',\n'molecule_atom_1_dist_min_div',\n'atom_index_0',\n'molecule_dist_max',\n'molecule_atom_1_dist_std_diff',\n'molecule_type_dist_max',\n'molecule_atom_index_0_y_1_max_diff',\n'molecule_type_0_dist_std_diff',\n'molecule_type_dist_mean_diff',\n'molecule_atom_1_dist_mean',\n'molecule_atom_index_0_y_1_mean_div',\n'molecule_type_dist_mean_div',\n'type']\n\nfor f in ['atom_1', 'type_0', 'type']:\n if f in good_columns:\n lbl = LabelEncoder()\n lbl.fit(list(train[f].values) + list(test[f].values))\n train[f] = lbl.transform(list(train[f].values))\n test[f] = lbl.transform(list(test[f].values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"X = train[good_columns].copy()\ny = train['scalar_coupling_constant']\ny_fc = train['fc']\nX_test = test[good_columns].copy()\n\n#del train, test\n#gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n# **5. Model** \n"},{"metadata":{"trusted":false},"cell_type":"code","source":"n_fold = 3\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create out of fold feature"},{"metadata":{"trusted":false},"cell_type":"code","source":"params = {'num_leaves': 128,\n 'min_child_samples': 79,\n 'objective': 'regression',\n 'max_depth': 9,\n 'learning_rate': 0.2,\n \"boosting_type\": \"gbdt\",\n \"subsample_freq\": 1,\n \"subsample\": 0.9,\n \"bagging_seed\": 11,\n \"metric\": 'mae',\n \"verbosity\": -1,\n 'reg_alpha': 0.1,\n 'reg_lambda': 0.3,\n 'colsample_bytree': 1.0\n }\nresult_dict_lgb_oof = train_model_regression(X=X, X_test=X_test, y=y_fc, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=False,\n verbose=500, early_stopping_rounds=200, n_estimators=1000)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"X['oof_fc'] = result_dict_lgb_oof['oof']\nX_test['oof_fc'] = result_dict_lgb_oof['prediction']","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def create_bunch_of_features(dtrain,dtest,cat_features):\n n_new_features = 0\n train_objs_num = len(dtrain)\n df_merge = pd.concat(objs=[dtrain, dtest], axis=0)\n \n for feature in cat_features:\n #Log Transform\n df_merge[feature+'_log'] = np.log (df_merge[feature])\n n_new_features = n_new_features +1\n\n dtrain = df_merge[:train_objs_num]\n dtest = df_merge[train_objs_num:]\n del df_merge\n gc.collect()\n print('Features Created: {} \\nTotal Features {}'.format(n_new_features,len(dtrain.columns)))\n return dtrain,dtest\n\nfeatures = list(X.columns)\nX, X_test = create_bunch_of_features(X,X_test,features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"good_columns = ['oof_fc_log',\n 'oof_fc',\n 'molecule_atom_index_0_dist_min',\n 'molecule_dist_mean',\n 'molecule_couples',\n 'molecule_atom_index_0_dist_std',\n 'molecule_atom_index_0_y_1_mean_div_log',\n 'molecule_atom_index_0_dist_mean',\n 'molecule_atom_index_0_dist_max',\n 'molecule_atom_1_dist_mean',\n 'molecule_dist_max',\n 'y_0',\n 'molecule_atom_index_0_z_1_std',\n 'molecule_atom_1_dist_std_diff',\n 'molecule_atom_index_0_x_1_std',\n 'molecule_atom_index_0_y_1_std',\n 'molecule_atom_index_1_dist_min',\n 'molecule_atom_index_1_dist_std',\n 'molecule_atom_index_0_y_1_mean_diff_log',\n 'molecule_atom_index_1_dist_mean',\n 'molecule_dist_min',\n 'molecule_atom_index_1_dist_max',\n 'molecule_type_dist_std_diff',\n 'molecule_type_dist_max',\n 'molecule_type_dist_mean_diff_log',\n 'molecule_type_dist_min',\n 'molecule_atom_index_0_y_1_max_diff',\n 'molecule_atom_index_0_y_1_mean_div',\n 'molecule_atom_index_0_y_1_mean_diff',\n 'molecule_atom_index_0_dist_mean_diff_log',\n 'molecule_type_0_dist_std_diff',\n 'molecule_atom_index_1_dist_std_diff',\n 'y_0_log',\n 'molecule_atom_index_1_dist_mean_diff_log',\n 'molecule_atom_index_0_dist_std_diff',\n 'molecule_atom_1_dist_min_diff',\n 'molecule_type_dist_mean_diff',\n 'molecule_atom_index_0_dist_max_div',\n 'molecule_atom_index_1_dist_min_diff',\n 'atom_index_1',\n 'molecule_atom_index_0_dist_std_div',\n 'atom_1_couples_count',\n 'molecule_atom_index_0_dist_mean_div',\n 'atom_index_0',\n 'molecule_atom_index_1_dist_std_div',\n 'molecule_atom_index_0_dist_max_diff',\n 'molecule_type_dist_mean_div',\n 'molecule_atom_index_1_dist_max_div',\n 'molecule_atom_index_1_dist_max_diff',\n 'molecule_atom_index_0_x_1_std_log',\n 'type']\n\nX = X[good_columns].copy()\nX_test = X_test[good_columns].copy()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model with oof feature"},{"metadata":{"trusted":false},"cell_type":"code","source":"params = {'num_leaves': 128,\n 'min_child_samples': 79,\n 'objective': 'regression',\n 'max_depth': 9,\n 'learning_rate': 0.2,\n \"boosting_type\": \"gbdt\",\n \"subsample_freq\": 1,\n \"subsample\": 0.9,\n \"bagging_seed\": 11,\n \"metric\": 'mae',\n \"verbosity\": -1,\n 'reg_alpha': 0.1,\n 'reg_lambda': 0.3,\n 'colsample_bytree': 1.0\n }\nresult_dict_lgb2 = train_model_regression(X=X, X_test=X_test, y=y, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=True,\n verbose=500, early_stopping_rounds=200, n_estimators=3000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#Best Features?\nfeature_importance = result_dict_lgb2['feature_importance']\nbest_features = feature_importance[['feature','importance']].groupby(['feature']).mean().sort_values(\n by='importance',ascending=False).iloc[:50,0:0].index.tolist()\nbest_features","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n# **6. Visualization and Analysis of Results** "},{"metadata":{"trusted":false},"cell_type":"code","source":"params = {'num_leaves': 128,\n 'min_child_samples': 79,\n 'objective': 'regression',\n 'max_depth': 9,\n 'learning_rate': 0.25,\n \"boosting_type\": \"gbdt\",\n \"subsample_freq\": 1,\n \"subsample\": 0.9,\n \"bagging_seed\": 11,\n \"metric\": 'mae',\n \"verbosity\": -1,\n 'reg_alpha': 0.1,\n 'reg_lambda': 0.3,\n 'colsample_bytree': 1.0\n }\nresult_dict_lgb2 = train_model_regression(X=X, X_test=X_test, y=y, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=True,\n verbose=500, early_stopping_rounds=200, n_estimators=1000)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training models for each type"},{"metadata":{"trusted":false},"cell_type":"code","source":"X_short = pd.DataFrame({'ind': list(X.index), 'type': X['type'].values, 'oof': [0] * len(X), 'target': y.values})\nX_short_test = pd.DataFrame({'ind': list(X_test.index), 'type': X_test['type'].values, 'prediction': [0] * len(X_test)})\nfor t in X['type'].unique():\n print(f'Training of type {t}')\n X_t = X.loc[X['type'] == t]\n X_test_t = X_test.loc[X_test['type'] == t]\n y_t = X_short.loc[X_short['type'] == t, 'target']\n result_dict_lgb3 = train_model_regression(X=X_t, X_test=X_test_t, y=y_t, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=False,\n verbose=500, early_stopping_rounds=200, n_estimators=3000)\n X_short.loc[X_short['type'] == t, 'oof'] = result_dict_lgb3['oof']\n X_short_test.loc[X_short_test['type'] == t, 'prediction'] = result_dict_lgb3['prediction']\n \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n"},{"metadata":{},"cell_type":"markdown","source":"
\n# **7. Submittion** "},{"metadata":{"trusted":false},"cell_type":"code","source":"#Deep model\nsub['scalar_coupling_constant'] = result_dict_lgb2['prediction']\nsub.to_csv('submission_deepmodel.csv', index=False)\nsub.head()\n\n#Training models for type\nsub['scalar_coupling_constant'] = X_short_test['prediction']\nsub.to_csv('submission_type.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n# **8. References** \n\n[1] OOF Model: https://www.kaggle.com/adarshchavakula/out-of-fold-oof-model-cross-validation
\n[2] Using Meta Features: https://www.kaggle.com/artgor/using-meta-features-to-improve-model
\n[3] Lot of Features: https://towardsdatascience.com/understanding-feature-engineering-part-1-continuous-numeric-data-da4e47099a7b"}],"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.4"}},"nbformat":4,"nbformat_minor":1}