{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\n\nThe features are taken from \n* [Giba R + data.table + Simple Features](https://www.kaggle.com/titericz/giba-r-data-table-simple-features-1-17-lb), which was in R and I took the data generated by it.\n* [Quantum Machine 9 - QM9](https://www.kaggle.com/zaharch/quantum-machine-9-qm9), which is now allowed in this competition.\n\nThe training template follows from Andrew's great starter: [Using meta-features to improve model](https://www.kaggle.com/artgor/using-meta-features-to-improve-model), except now the CV is `GroupKFold` using molecule's names, and only essential libraries and functions are loaded.\n\nThis kernel also illustrate how to code a customized loss for LightGBM."},{"metadata":{},"cell_type":"markdown","source":"## Load libs and funcs"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()\n\nfrom tqdm import tqdm_notebook\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn import metrics\npd.options.display.precision = 15\n\nimport lightgbm as lgb\nimport time\nimport datetime\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport gc","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def reduce_mem_usage(df, verbose=False):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)    \n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df\n\n\n\ndef group_mean_log_mae(y_true, y_pred, types, floor=1e-9):\n    \"\"\"\n    Fast metric computation for this competition: https://www.kaggle.com/c/champs-scalar-coupling\n    Code is from this kernel: https://www.kaggle.com/uberkinder/efficient-metric\n    \"\"\"\n    maes = (y_true-y_pred).abs().groupby(types).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()\n\ndef train_lgb_regression_group(X, X_test, y, params, folds, groups,\n                               eval_metric='mae', \n                               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 - Group Kfolds 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                        'sklearn_scoring_function': metrics.mean_absolute_error},\n                    'group_mae': {'lgb_metric_name': 'mae',\n                        'scoring_function': group_mean_log_mae},\n                    'mse': {'lgb_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,groups=groups)):\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        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)], \n                  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        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 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:.6f}, std: {1:.6f}.\\n'.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    \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\", \n                    data=best_features.sort_values(by=\"importance\", ascending=False));\n        plt.title('LGB Features (avg over folds)');\n\n        result_dict['feature_importance'] = feature_importance\n        \n    return result_dict","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"giba_columns = ['inv_dist0', 'inv_dist1', 'inv_distP', 'inv_dist0R', 'inv_dist1R', 'inv_distPR',\n 'inv_dist0E', 'inv_dist1E', 'inv_distPE', 'linkM0', 'linkM1',\n 'min_molecule_atom_0_dist_xyz',\n 'mean_molecule_atom_0_dist_xyz',\n 'max_molecule_atom_0_dist_xyz',\n 'sd_molecule_atom_0_dist_xyz',\n 'min_molecule_atom_1_dist_xyz',\n 'mean_molecule_atom_1_dist_xyz',\n 'max_molecule_atom_1_dist_xyz',\n 'sd_molecule_atom_1_dist_xyz',\n 'coulomb_C.x', 'coulomb_F.x', 'coulomb_H.x', 'coulomb_N.x', 'coulomb_O.x',\n 'yukawa_C.x', 'yukawa_F.x', 'yukawa_H.x', 'yukawa_N.x', 'yukawa_O.x',\n 'vander_C.x', 'vander_F.x', 'vander_H.x', 'vander_N.x', 'vander_O.x',\n 'coulomb_C.y', 'coulomb_F.y', 'coulomb_H.y', 'coulomb_N.y', 'coulomb_O.y',\n 'yukawa_C.y', 'yukawa_F.y', 'yukawa_H.y', 'yukawa_N.y', 'yukawa_O.y',\n 'vander_C.y', 'vander_F.y', 'vander_H.y', 'vander_N.y', 'vander_O.y',\n 'distC0', 'distH0', 'distN0', 'distC1', 'distH1', 'distN1',\n 'adH1', 'adH2', 'adH3', 'adH4', 'adC1', 'adC2', 'adC3', 'adC4',\n 'adN1', 'adN2', 'adN3', 'adN4',\n 'NC', 'NH', 'NN', 'NF', 'NO']\n\nqm9_columns = [\n'rc_A', 'rc_B', 'rc_C', \n'mu', 'alpha', \n'homo','lumo', 'gap', \n'zpve', 'Cv', \n'freqs_min', 'freqs_max', 'freqs_mean',\n'mulliken_min', 'mulliken_max', \n'mulliken_atom_0', 'mulliken_atom_1'\n]\n\nlabel_columns = ['molecule_name',\n'atom_index_0', 'atom_index_1',\n'structure_atom_0','structure_atom_1','type']\n\nindex_columns = ['type','molecule_name','id']\n\ndiff_columns = ['Cv',\n 'alpha', 'freqs_max', 'freqs_mean', 'freqs_min',\n 'gap', 'homo', 'linkM0',\n 'lumo', 'mu', 'mulliken_atom_0', 'mulliken_max', 'mulliken_min',\n 'rc_A', 'rc_B', 'rc_C', 'sd_molecule_atom_1_dist_xyz', 'zpve']","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"%%time\nfile_folder = '../input/champs-scalar-coupling'\n\nprint(\"Load Giba's features...\")\ntrain = pd.read_csv('../input/giba-molecular-features/train_giba.csv/train_giba.csv',\n                   usecols=index_columns+giba_columns+['scalar_coupling_constant'])\ny = train['scalar_coupling_constant']\ntrain = reduce_mem_usage(train,verbose=True)\n\ntest = pd.read_csv('../input/giba-molecular-features/test_giba.csv/test_giba.csv',\n                  usecols=index_columns+giba_columns)\ntest = reduce_mem_usage(test,verbose=True)\n\nprint(\"Load QM9 features...\")\ndata_qm9 = pd.read_pickle('../input/quantum-machine-9-qm9/data.covs.pickle')\ndata_qm9 = data_qm9.drop(columns = ['type', 'linear', 'atom_index_0', 'atom_index_1', \n            'scalar_coupling_constant', 'U', 'G', 'H', \n            'mulliken_mean', 'r2', 'U0'], axis=1)\ndata_qm9 = reduce_mem_usage(data_qm9,verbose=False)\n\ntrain = pd.merge(train, data_qm9, how='left', on=['molecule_name','id'])\ntest = pd.merge(test, data_qm9, how='left', on=['molecule_name','id'])\n\ndel data_qm9\ngc.collect()\n\nprint(\"Encoding label features...\\n\")\nfor f in label_columns: \n    # 'type' has to be the last one\n    # since the this label encoder is used later\n    if f in train.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))\n        \ntrain = train[index_columns+giba_columns+qm9_columns]\ntest = test[index_columns+giba_columns+qm9_columns]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train by type using GroupKFold using a customized objective\n\n* A weighted average of `fair`, `huber` (smoothened version of $L^1$), $L^2$ (regression), and $L^1$ losses.\n* The fold is a GroupKfold by molecule's names."},{"metadata":{"trusted":true},"cell_type":"code","source":"coef = [0.25, 0.5, 0.2, 0.05]\n\ndef custom_objective(y_true, y_pred):\n    \n    # fair\n    c = 0.5\n    residual = y_pred - y_true\n    grad = c * residual /(np.abs(residual) + c)\n    hess = c ** 2 / (np.abs(residual) + c) ** 2\n    \n    # huber\n    h = 1.2  #h is delta in the Huber's formula\n    scale = 1 + (residual / h) ** 2\n    scale_sqrt = np.sqrt(scale)\n    grad_huber = residual / scale_sqrt\n    hess_huber = 1 / scale / scale_sqrt\n\n    # rmse grad and hess\n    grad_rmse = residual\n    hess_rmse = 1.0\n\n    # mae grad and hess\n    grad_mae = np.array(residual)\n    grad_mae[grad_mae > 0] = 1.\n    grad_mae[grad_mae <= 0] = -1.\n    hess_mae = 1.0\n\n    return coef[0] * grad + coef[1] * grad_huber + coef[2] * grad_rmse + coef[3] * grad_mae, \\\n           coef[0] * hess + coef[1] * hess_huber + coef[2] * hess_rmse + coef[3] * hess_mae\n\nparams = {\n'num_leaves': 400,\n'objective': custom_objective,\n'max_depth': 9,\n'learning_rate': 0.1,\n'boosting_type': 'gbdt',\n'metric': 'mae',\n'verbosity': -1,\n'lambda_l1': 2,\n'lambda_l2': 0.2,\n'feature_fraction': 0.6,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_short = pd.DataFrame({'ind': list(train.index), \n                        'type': train['type'].values,\n                        'oof': [0] * len(train), \n                        'target': y.values})\nX_short_test = pd.DataFrame({'ind': list(test.index), \n                             'type': test['type'].values, \n                             'prediction': [0] * len(test)})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nCV_score = 0\nfolds = GroupKFold(n_splits=3)\n####Iters####  [1JHC, 1JHN, 2JHC, 2JHH, 2JHN, 3JHC, 3JHH, 3JHN]\nn_estimators = [6000, 2500, 3500, 3000, 3000, 5000, 3000, 3000]\n\nfor t in train['type'].unique():\n    type_ = lbl.inverse_transform([t])[0]\n    print(f'Training of type {t}: {type_}.\\n')\n    X_t = train.loc[train['type'] == t]\n    X_test_t = test.loc[test['type'] == t]\n    y_t = X_short.loc[X_short['type'] == t, 'target']\n    \n    scaler = StandardScaler()\n    X_t[diff_columns] = scaler.fit_transform(X_t[diff_columns].fillna(-999))\n    X_test_t[diff_columns] = scaler.transform(X_test_t[diff_columns].fillna(-999))\n    \n    molecules_id = X_t.molecule_name\n    \n    result_dict_lgb = train_lgb_regression_group(X=X_t.drop(columns=['molecule_name','id']), \n                                          X_test=X_test_t.drop(columns=['molecule_name','id']), \n                                          y=y_t, params=params, \n                                          folds=folds, groups=molecules_id,\n                                          eval_metric='group_mae', \n                                          plot_feature_importance=False,\n                                          verbose=3000, early_stopping_rounds=200, \n                                          n_estimators=n_estimators[t])\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    CV_score += np.array(result_dict_lgb['scores']).mean()/8 # total 8 types","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(f'{file_folder}/sample_submission.csv')\nsub['scalar_coupling_constant'] = X_short_test['prediction']\ntoday = str(datetime.date.today())\nsub.to_csv(f'LGB_{today}_{CV_score:.4f}.csv', index=False)","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}