{"cells":[{"metadata":{},"cell_type":"markdown","source":"## General information\n\nWe all know that feature engineering is the most important things in competitions. If you have deep domain knowledge, you can create better features!\n\nSadly, I have little domain knowledge (last time I paid attention to chemical formulas was in school)...\n\nSo I'll to create a lot of features and hope that some of them work!\n\n**Important!**\nI have realized the main problem of features in my previous [kernel](https://www.kaggle.com/artgor/molecular-properties-eda-and-models): I created features based on the whole train/test set, but it makes much more sense to generate features based on molecules (for example mean distance per molecules). Let's try this!\n\n![](https://imgs.xkcd.com/comics/machine_learning.png)\n\n~Thanks to the new kaggle update we can write code in kernels and import it. This is much more convenient and useful. I'm moving all the functions I can into this script: https://www.kaggle.com/artgor/artgor-utils So if you see somewhere code like artgot_utils.function_name(parameters) - it is from this script~\n**Important**:\nIt seems that after forking the utility script isn't copied, so I have decided to move all the functions back into the kernel. But utility kernel can't be removed from the kernel... I hope I was able to deal with it."},{"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":"importing 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 functions"},{"metadata":{"_kg_hide-input":true,"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\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/champs-scalar-coupling' if 'champs-scalar-coupling' in os.listdir('../input/') else '../input'\nos.listdir(file_folder)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a lot of files, let's focus on the main ones for now."},{"metadata":{"trusted":true},"cell_type":"code","source":"np.set_printoptions(precision=15)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = 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')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'There are {train.shape[0]} rows in train data.')\nprint(f'There are {test.shape[0]} rows in test data.')\n\nprint(f\"There are {train['molecule_name'].nunique()} distinct molecules in train data.\")\nprint(f\"There are {test['molecule_name'].nunique()} distinct molecules in test data.\")\nprint(f\"There are {train['atom_index_0'].nunique()} unique atoms.\")\nprint(f\"There are {train['type'].nunique()} unique types.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Basic feature engineering\n\nCode for joining data is from this kernel:\nhttps://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark"},{"metadata":{"trusted":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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"distances"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I want to create two new features: one will show the first character of the `type`, the second one will show other characters."},{"metadata":{"trusted":true},"cell_type":"code","source":"train['type_0'] = train['type'].apply(lambda x: x[0])\ntest['type_0'] = test['type'].apply(lambda x: x[0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I calculated these three features in my previous kernel... But they are hardly useful - they contain statistics on all the dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = reduce_mem_usage(train)\ntest = reduce_mem_usage(test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Much more features\n\nI have decided to create a lot of features:\n- these features calculate statistics per molecules;\n- several features are calculated for molecules in general;\n- most features are calculated for molecules and one of categorical variables;\n- I wanted to use all combinations of features and more aggregates, but kernel RAM isn't enough :( \n- also I wanted to calculated differences between aggregated and original values\n\nYou can see the code to generate all features (400+) below in function `create_features_full`. I generated all these features, trained an lgb and selected top-50 features. I'll generate them now.\n\nI'd upload all the generated features on Kaggle dataset, but the total size of train and test data exceeds the allowed limits :("},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_features_full(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['molecule_dist_std'] = df.groupby('molecule_name')['dist'].transform('std')\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    num_cols = ['x_1', 'y_1', 'z_1', 'dist', 'dist_x', 'dist_y', 'dist_z']\n    cat_cols = ['atom_index_0', 'atom_index_1', 'type', 'atom_1', 'type_0']\n    aggs = ['mean', 'max', 'std', 'min']\n    for col in cat_cols:\n        df[f'molecule_{col}_count'] = df.groupby('molecule_name')[col].transform('count')\n\n    for cat_col in tqdm_notebook(cat_cols):\n        for num_col in num_cols:\n            for agg in aggs:\n                df[f'molecule_{cat_col}_{num_col}_{agg}'] = df.groupby(['molecule_name', cat_col])[num_col].transform(agg)\n                df[f'molecule_{cat_col}_{num_col}_{agg}_diff'] = df[f'molecule_{cat_col}_{num_col}_{agg}'] - df[num_col]\n                df[f'molecule_{cat_col}_{num_col}_{agg}_div'] = df[f'molecule_{cat_col}_{num_col}_{agg}'] / df[num_col]\n\n    df = reduce_mem_usage(df)\n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Generating subsert of features for this kernel."},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = create_features(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = create_features(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"good_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']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot new features\n\nFor now let's have a look at one of the new features."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 10))\nfor i, t in enumerate(train['type'].unique()):\n    plt.scatter(train.loc[train['type'] == t, good_columns[0]], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\nplt.title(f'{good_columns[0]} vs target for 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, good_columns[0]], train.loc[train['type'] == t, 'scalar_coupling_constant'], label=t);\n    plt.title(f'{good_columns[0]} vs target \\n for {t} 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.hist(train.loc[train['type'] == t, good_columns[0]], label='train');\n    plt.hist(test.loc[test['type'] == t, good_columns[0]], label='test');\n    plt.title(f'{good_columns[0]} distribution \\n for {t} type');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Prepare data for model\n\nI'll use the function for metric calculation from this kernel: https://www.kaggle.com/abhishek/competition-metric\n\nUPD: use faster metric from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n\nYou can see the code in my utility script. Please note that to use this metric for calculation you need to pass value `group_mae` to parameter `eval_metric`.\n\nFor now I'll simply train a single model on folds"},{"metadata":{"trusted":true},"cell_type":"code","source":"for 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']\nX_test = test[good_columns].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train, test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training model on selected features\n\nFor now I'll train model with simply 1500 estimators to check that it works."},{"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=5000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['scalar_coupling_constant'] = result_dict_lgb['prediction']\nsub.to_csv('submission.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(result_dict_lgb['prediction'], bins=40);\nplt.title('Distribution of predictions');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot oof predictions vs target\n\nI use the code from benchmark kernel: https://www.kaggle.com/inversion/atomic-distance-benchmark\n\nNotice that while using `LabelEncoder` on categorical variables, `type` was the last feature, so `lbl` contains encodings for it and we can inverse_transform `type`."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_data = pd.DataFrame(y)\nplot_data.index.name = 'id'\nplot_data['yhat'] = result_dict_lgb['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":{},"cell_type":"markdown","source":"We can see that some types has better predictions than others. Maybe we should focus on improving models for these types?"},{"metadata":{},"cell_type":"markdown","source":"## Training separate models for each type\n\nWe can see that different types have really different distributions of target variable, so it seems to be reasonable to build 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_lgb = 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=2000)\n    X_short.loc[X_short['type'] == t, 'oof'] = result_dict_lgb['oof']\n    X_short_test.loc[X_short_test['type'] == t, 'prediction'] = result_dict_lgb['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":{"trusted":true},"cell_type":"code","source":"plt.hist(result_dict_lgb['prediction'], bins=40);\nplt.title('Distribution of predictions');","execution_count":null,"outputs":[]},{"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_short['type'])\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}