{"cells":[{"metadata":{},"cell_type":"markdown","source":"## General information\n\nThis kernel is a continuation of my previous kernels:\n- [EDA](https://www.kaggle.com/artgor/molecular-properties-eda-and-models)\n- [feature engineering](https://www.kaggle.com/artgor/brute-force-feature-engineering)\n\nIn this kernel I'll try to use additional data, which is provided only for train dataset.\n\nWe have several files with additional features, which have information only for train and not for test. In this kernel I'll show how to create meta-features both train and test and use them."},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install -U vega_datasets notebook vega","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"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\n\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","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":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\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 v5\nfrom IPython.display import HTML\n\n# 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\n    vega_url = 'https://cdn.jsdelivr.net/npm/vega@' + v5.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 <div/> 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    <div id=\"{id}\"></div><script>\n    require([\"vega-embed\"], function(vg_embed) {{\n        const spec = {chart};     \n        vg_embed(\"#{id}\", spec, {{defaultStyle: true}}).catch(console.warn);\n        console.log(\"anything?\");\n    }});\n    console.log(\"really...anything?\");\n    </script>\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                c_prec = df[col].apply(lambda x: np.finfo(x).precision).max()\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max and c_prec == np.finfo(np.float16).precision:\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 and c_prec == np.finfo(np.float32).precision:\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    \"<script>\",\n    workaround,\n    \"</script>\",\n)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data loading and overview"},{"metadata":{"trusted":true},"cell_type":"code","source":"file_folder = '../input'\ntrain = pd.read_csv(f'{file_folder}/train.csv')\ntest = pd.read_csv(f'{file_folder}/test.csv')\nsub = pd.read_csv(f'{file_folder}/sample_submission.csv')\nstructures = pd.read_csv(f'{file_folder}/structures.csv')\npotential_energy = pd.read_csv(f'{file_folder}/potential_energy.csv')\nmulliken_charges = pd.read_csv(f'{file_folder}/mulliken_charges.csv')\nscalar_coupling_contributions = pd.read_csv(f'{file_folder}/scalar_coupling_contributions.csv')\nmagnetic_shielding_tensors = pd.read_csv(f'{file_folder}/magnetic_shielding_tensors.csv')\ndipole_moments = pd.read_csv(f'{file_folder}/dipole_moments.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's start with `scalar_coupling_contributions`"},{"metadata":{"trusted":true},"cell_type":"code","source":"scalar_coupling_contributions.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we have a direct relationship between this file and train - each line in train has line in this file. We have 4 additional variables here, let's see what is written in competition description for this file:\n```\nThe scalar coupling constants in train.csv are a sum of four terms. scalar_coupling_contributions.csv contain all these terms.\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```\n\nIt seems that these 4 features should be highly correlated with the target variable, so using them should improve our models!\nLet's merge these file and have a look at them!"},{"metadata":{"trusted":true},"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":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.subplot(2, 4, i + 1);\n    plt.scatter(train.loc[train['type'] == t, 'fc'], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'fc vs target \\n for {t} type');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have an amazing correlation!\n\nSo here what I will do:\n\n- generate features from my previous kernels;\n- train a model on them with 1500 estimators for speed;\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":"## Feature generation\n\nIn the hidden cell below I generate features as in my kernel: https://www.kaggle.com/artgor/brute-force-feature-engineering"},{"metadata":{"trusted":true,"_kg_hide-input":true},"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\ntrain = reduce_mem_usage(train)\ntest = reduce_mem_usage(test)\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    \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\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":true},"cell_type":"code","source":"X = train[good_columns].copy()\ny = train['scalar_coupling_constant']\ny_fc = train['fc']\nX_test = test[good_columns].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## On importance of the generated features\n\nAfter the training below you'll see that some of the most important variables are these: `molecule_atom_index_0_dist_max`, `molecule_atom_index_0_dist_mean`, `molecule_atom_index_0_dist_std`. Let's plot them!"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.subplot(2, 4, i + 1);\n    plt.scatter(train.loc[train['type'] == t, 'dist'], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'dist vs target \\n for {t} type');\n    plt.xlim(train['dist'].min(), train['dist'].max());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.subplot(2, 4, i + 1);\n    plt.scatter(train.loc[train['type'] == t, 'molecule_atom_index_0_dist_max'], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'molecule_atom_index_0_dist_max vs target \\n for {t} type');\n    plt.xlim(train['molecule_atom_index_0_dist_max'].min(), train['molecule_atom_index_0_dist_max'].max());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.subplot(2, 4, i + 1);\n    plt.scatter(train.loc[train['type'] == t, 'molecule_atom_index_0_dist_mean'], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'molecule_atom_index_0_dist_mean vs target \\n for {t} type');\n    plt.xlim(train['molecule_atom_index_0_dist_mean'].min(), train['molecule_atom_index_0_dist_mean'].max());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.subplot(2, 4, i + 1);\n    plt.scatter(train.loc[train['type'] == t, 'molecule_atom_index_0_dist_std'], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'molecule_atom_index_0_dist_std vs target \\n for {t} type');\n    plt.xlim(train['molecule_atom_index_0_dist_std'].min(), train['molecule_atom_index_0_dist_std'].max());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train, test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Basic model\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 3\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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 = 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=1500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create out of fold feature"},{"metadata":{"trusted":true},"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=1500)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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":{},"cell_type":"markdown","source":"## Model with oof feature"},{"metadata":{"trusted":true},"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=1500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Deeper model"},{"metadata":{"trusted":true},"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=3000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['scalar_coupling_constant'] = result_dict_lgb2['prediction']\nsub.to_csv('submission.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training separate models for each type"},{"metadata":{"trusted":true},"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    \nsub['scalar_coupling_constant'] = X_short_test['prediction']\nsub.to_csv('submission_t.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot oof predictions vs target"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_data = pd.DataFrame(y)\nplot_data.index.name = 'id'\nplot_data['yhat'] = X_short['oof']\nplot_data['type'] = lbl.inverse_transform(X['type'])\n\ndef plot_oof_preds(ctype, llim, ulim):\n        plt.figure(figsize=(6,6))\n        sns.scatterplot(x='scalar_coupling_constant',y='yhat',\n                        data=plot_data.loc[plot_data['type']==ctype,\n                        ['scalar_coupling_constant', 'yhat']]);\n        plt.xlim((llim, ulim))\n        plt.ylim((llim, ulim))\n        plt.plot([llim, ulim], [llim, ulim])\n        plt.xlabel('scalar_coupling_constant')\n        plt.ylabel('predicted')\n        plt.title(f'{ctype}', fontsize=18)\n        plt.show()\n\nplot_oof_preds('1JHC', 0, 250)\nplot_oof_preds('1JHN', 0, 100)\nplot_oof_preds('2JHC', -50, 50)\nplot_oof_preds('2JHH', -50, 50)\nplot_oof_preds('2JHN', -25, 25)\nplot_oof_preds('3JHC', -25, 100)\nplot_oof_preds('3JHH', -20, 20)\nplot_oof_preds('3JHN', -15, 15)","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}