{"cells":[{"metadata":{},"cell_type":"markdown","source":"## General information\n\nThis kernel is created using data of `Predicting Molecular Properties` competition.\nWe have information about atom couples in molecules and need to predict `scalar_coupling_constant` between these atoms.\n\n![](http://www.et.byu.edu/~rowley/VLEfinal/methane_dimer.gif)\n\nIn this kernel I'll do EDA and will try some approaches to modelling.\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":{},"cell_type":"markdown","source":"importing libraries"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install -U vega_datasets notebook vega","execution_count":null,"outputs":[]},{"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 the fuctions"},{"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.float64).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    :params: verbose - parameters for gradient boosting models\n    :params: early_stopping_rounds - parameters for gradient boosting models\n    :params: n_estimators - parameters for gradient boosting models\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 classification 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    :params: verbose - parameters for gradient boosting models\n    :params: early_stopping_rounds - parameters for gradient boosting models\n    :params: n_estimators - parameters for gradient boosting models\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":{"_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 {structures['atom'].nunique()} unique atoms.\")\nprint(f\"There are {train['type'].nunique()} unique types.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So in out main data files we have information about moleculas and pairs of atoms\n- test set in ~2 times smaller that train set;\n- we have 28 unique atoms types and 8 coupling types;"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"atom_count = train['atom_index_0'].value_counts().reset_index().rename(columns={'atom_index_0': 'count', 'index': 'atom_index_0'})\nchart1 = alt.Chart(atom_count).mark_bar().encode(\n    x=alt.X(\"atom_index_0:N\", axis=alt.Axis(title='atom_index_0')),\n    y=alt.Y('count:Q', axis=alt.Axis(title='Count')),\n    tooltip=['atom_index_0', 'count']\n).properties(title=\"Counts of atom_index_0\", width=350).interactive()\n\natom_count = train['atom_index_1'].value_counts().reset_index().rename(columns={'atom_index_1': 'count', 'index': 'atom_index_1'})\nchart2 = alt.Chart(atom_count).mark_bar().encode(\n    x=alt.X(\"atom_index_1:N\", axis=alt.Axis(title='atom_index_1')),\n    y=alt.Y('count:Q', axis=alt.Axis(title='Count')),\n    tooltip=['atom_index_1', 'count']\n).properties(title=\"Counts of atom_index_1\", width=350).interactive()\n\ntype_count = train['type'].value_counts().reset_index().rename(columns={'type': 'count', 'index': 'type'})\nchart3 = alt.Chart(type_count).mark_bar().encode(\n    x=alt.X(\"type:N\", axis=alt.Axis(title='type')),\n    y=alt.Y('count:Q', axis=alt.Axis(title='Count')),\n    tooltip=['type', 'count']\n).properties(title=\"Counts of type\", width=350).interactive()\n\nhist_df = pd.cut(train['scalar_coupling_constant'], 20).value_counts().sort_index().reset_index().rename(columns={'index': 'bins'})\nhist_df['bins'] = hist_df['bins'].astype(str)\nchart4 = alt.Chart(hist_df).mark_bar().encode(\n    x=alt.X(\"bins:O\", axis=alt.Axis(title='Target bins')),\n    y=alt.Y('scalar_coupling_constant:Q', axis=alt.Axis(title='Count')),\n    tooltip=['scalar_coupling_constant', 'bins']\n).properties(title=\"scalar_coupling_constant histogram\", width=400).interactive()\n\n\nrender((chart1 | chart2) & (chart3 | chart4))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (18, 6))\nplt.subplot(1, 2, 1);\nplt.hist(train['scalar_coupling_constant'], bins=20);\nplt.title('Basic scalar_coupling_constant histogram');\nplt.subplot(1, 2, 2);\nsns.violinplot(x='type', y='scalar_coupling_constant', data=train);\nplt.title('Violinplot of scalar_coupling_constant by type');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are many interesting things here:\n- among first atoms there is a little number of atoms with index lower than 7 or higher than 24;\n- among second atoms there is a little number of atoms with index higher than 24. Also index with atom with index 9 in quite rare;\n- coupling types are unevenly distributed. There are 3 very popular, 3 quite rare and 2 with medium frequency;\n- target variable has a bimodal distribution;\n- different coupling types have really different values of target variable. Maybe it would make sense to build separate models for each of them;"},{"metadata":{},"cell_type":"markdown","source":"## Plotting network graphs by type\n\nWe have molecules, atom pairs, so this means data, which is interconnected. Network graphs should be useful to visualize such data!"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 12))\nfor i, t in enumerate(train['type'].unique()):\n    train_type = train.loc[train['type'] == t]\n    G = nx.from_pandas_edgelist(train_type, 'atom_index_0', 'atom_index_1', ['scalar_coupling_constant'])\n    plt.subplot(2, 4, i + 1);\n    nx.draw(G, with_labels=True);\n    plt.title(f'Graph for type {t}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Cool! We can see that atom connections have different shapes for different types. Type 2JHH has an expecially unique scheme.\nAlso we can see that some atoms are connected only to several other atoms."},{"metadata":{},"cell_type":"markdown","source":"## Better network graphs\nBut there is a little problem: as we saw earlier, there are atoms which are very rare, as a result graphs will be skewed due to them. Now I'll drop atoms for each type which are present in less then 1% of connections"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (20, 12))\nfor i, t in enumerate(train['type'].unique()):\n    train_type = train.loc[train['type'] == t]\n    bad_atoms_0 = list(train_type['atom_index_0'].value_counts(normalize=True)[train_type['atom_index_0'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms_1 = list(train_type['atom_index_1'].value_counts(normalize=True)[train_type['atom_index_1'].value_counts(normalize=True) < 0.01].index)\n    bad_atoms = list(set(bad_atoms_0 + bad_atoms_1))\n    train_type = train_type.loc[(train_type['atom_index_0'].isin(bad_atoms_0) == False) & (train_type['atom_index_1'].isin(bad_atoms_1) == False)]\n    G = nx.from_pandas_edgelist(train_type, 'atom_index_0', 'atom_index_1', ['scalar_coupling_constant'])\n    plt.subplot(2, 4, i + 1);\n    nx.draw(G, with_labels=True);\n    plt.title(f'Graph for type {t}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now the graphs are much more clear!"},{"metadata":{},"cell_type":"markdown","source":"## Feature engineering\n\nFor now I'll use a basic approach to feature engineering.\n\nI'll use this useful 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":{"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])\ntrain['type_1'] = train['type'].apply(lambda x: x[1:])\ntest['type_1'] = test['type'].apply(lambda x: x[1:])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So, everyone uses this distance feature, let's have a look at it!"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (18, 8))\nplt.subplot(1, 2, 1);\nplt.hist(train['dist'], bins=20);\nplt.title('Basic dist_speedup histogram');\nplt.subplot(1, 2, 2);\nsns.violinplot(x='type', y='dist', data=train);\nplt.title('Violinplot of dist_speedup by type');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems that the values are quite different for different types!"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['dist_to_type_mean'] = train['dist'] / train.groupby('type')['dist'].transform('mean')\ntest['dist_to_type_mean'] = test['dist'] / test.groupby('type')['dist'].transform('mean')\n\ntrain['dist_to_type_0_mean'] = train['dist'] / train.groupby('type_0')['dist'].transform('mean')\ntest['dist_to_type_0_mean'] = test['dist'] / test.groupby('type_0')['dist'].transform('mean')\n\ntrain['dist_to_type_1_mean'] = train['dist'] / train.groupby('type_1')['dist'].transform('mean')\ntest['dist_to_type_1_mean'] = test['dist'] / test.groupby('type_1')['dist'].transform('mean')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[f'molecule_type_dist_mean'] = train.groupby(['molecule_name', 'type'])['dist'].transform('mean')\ntest[f'molecule_type_dist_mean'] = test.groupby(['molecule_name', 'type'])['dist'].transform('mean')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Basic 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`."},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in ['atom_0', 'atom_1', 'type_0', 'type_1', 'type']:\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":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train.drop(['id', 'molecule_name', 'scalar_coupling_constant'], axis=1)\ny = train['scalar_coupling_constant']\nX_test = test.drop(['id', 'molecule_name'], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\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': 13,\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=1000, early_stopping_rounds=200, n_estimators=10000)\n","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":{"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}