{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Scalar coupling constant prediction using XGB with scikit-optimizer\n\nIn this kernel, we use features suggested in https://www.kaggle.com/artgor/brute-force-feature-engineering using the code from https://www.kaggle.com/kabure/simple-eda-lightgbm-autotuning-w-hyperopt/notebook\n\n\n* [Additional methods](#addmethods)\n* [Feature generation](#features)\n* [Dataset preprocessing](#dataprep)\n* [Feature selection, algorithm tuning, model training and prediction per molecule type](#preds)\n\n22.07.2019"},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\n\nimport multiprocessing\n\nimport numpy as np \nimport pandas as pd\npd.set_option(\"display.max_columns\", 999)\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm_notebook as tqdm\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.preprocessing import LabelEncoder\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import r2_score\n\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.model_selection import RandomizedSearchCV\nfrom sklearn.model_selection import cross_val_score, cross_validate\nfrom sklearn.model_selection import learning_curve, validation_curve\nfrom sklearn.model_selection import KFold, ShuffleSplit, StratifiedShuffleSplit\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.ensemble import RandomForestRegressor\n\nfrom time import time, ctime\n\nimport xgboost\n\n# HPO\nfrom skopt.space import Integer, Categorical, Real\nfrom skopt.utils import use_named_args\nfrom skopt import gp_minimize, gbrt_minimize\nfrom skopt.plots import plot_convergence\nfrom skopt.callbacks import DeltaXStopper, DeadlineStopper, DeltaYStopper\nfrom skopt.callbacks import EarlyStopper","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Additional methods** <a class=\"anchor\" id=\"addmethods\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_params_SKopt(model, X, Y, space):\n    cv_search = KFold(n_splits=3, shuffle=True, random_state = 0)\n    \n    HPO_PARAMS = {'n_calls':500,\n                  'n_random_starts':10,\n                  'acq_func':'EI',}\n\n    @use_named_args(space)\n    def objective(**params):\n        model.set_params(**params)\n        return -np.mean(cross_val_score(model, \n                                        X, Y, \n                                        cv=cv_search, \n                                        n_jobs = -1, \n                                        scoring='neg_mean_absolute_error'))\n    \n    reg_gp = gbrt_minimize(objective, \n                           space, \n                           verbose = False,\n                           callback = [RepeatedMinStopper(n_best = 30), DeadlineStopper(total_time = 7200)],\n                           **HPO_PARAMS,\n                           random_state = 0)\n\n    model.max_depth = reg_gp.x[0]\n    model.min_child_weight = reg_gp.x[1]\n    model.learning_rate = reg_gp.x[2]\n    model.subsample = reg_gp.x[3]\n    model.colsample_bytree = reg_gp.x[4]\n    model.reg_alpha = reg_gp.x[5]\n    model.reg_lambda = reg_gp.x[6]\n\n    return [model,reg_gp]\n\ndef plotfig (ypred, yactual, strtitle, y_max, reg_gp):\n    fig, ax = plt.subplots(1, 2, figsize=(15,5))\n\n    ax[0].scatter(ypred, yactual.values.ravel())\n    ax[0].set_title(strtitle)\n    ax[0].plot([(0, 0), (y_max, y_max)], [(0, 0), (y_max, y_max)])\n    ax[0].set_xlim(0, y_max)\n    ax[0].set_ylim(0, y_max)\n    ax[0].set_xlabel('Predicted', fontsize=12)\n    ax[0].set_ylabel('Actual', fontsize=12)\n\n    plot_convergence(reg_gp, ax = ax[1]) \n    plt.show()\n\ndef feature_selection(threshold, train, test):\n    corr_matrix = train.corr().abs()\n    upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n    to_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n    print('\\nThere are %d columns to remove.' % (len(to_drop)))\n    train = train.drop(columns = to_drop)\n    test = test.drop(columns = to_drop)  \n    print (f'After dropping {train.shape[1]}' + ' features remain')   \n    return [train, test, to_drop]\n\nclass RepeatedMinStopper(EarlyStopper):\n    \"\"\"Stop the optimization when there is no improvement in the minimum.\n    Stop the optimization when there is no improvement in the minimum\n    achieved function evaluation after `n_best` iterations.\n    \"\"\"\n    def __init__(self, n_best=50):\n        super(EarlyStopper, self).__init__()\n        self.n_best = n_best\n        self.count = 0\n        self.minimum = np.finfo(np.float).max\n\n    def _criterion(self, result):\n        if result.fun < self.minimum:\n            self.minimum = result.fun\n            self.count = 0\n        elif result.fun > self.minimum:\n            self.count = 0\n        else:\n            self.count += 1\n\n        return self.count >= self.n_best\n\ndef get_nan_col(df, N):\n    # Get features with minimum N percentage of null observations \n    n_observ = int(np.round(N*np.size(df, 0)))\n    allnull  = df.isnull().sum(axis=0).reset_index()\n    lst_del  = [allnull.loc[i,'index'] for i in range(len(allnull)) if allnull.loc[i,0] > n_observ]  \n    lst_proc = [allnull.loc[i,'index'] for i in range(len(allnull)) if allnull.loc[i,0] < n_observ and allnull.loc[i,0] > 0]\n    return [lst_del, lst_proc]\n\ndef smart_fillna (common_df, Y, percent , fill_method_all, model_type, cv, scoring):   \n    X = pd.DataFrame()\n    X_test = pd.DataFrame()\n\n    if model_type == 'linear':\n        model = LinearRegression()\n    elif model_type == 'rfr':\n        model = RandomForestRegressor(n_estimators = 50, n_jobs = -1)\n    \n    lst_nan = get_nan_col(common_df, percent)\n    print (f'To delete {len(lst_nan[0])}' + ' features')\n    print (f'To process {len(lst_nan[1])}' + ' features')\n    common_df.drop(lst_nan[0], axis = 1, inplace = True)\n       \n    if len(lst_nan[1]) > 0:\n        print ('Processing features...')\n        for feature in tqdm(lst_nan[1]):\n            mas_score = []\n            best_score = np.inf\n            best_method = ''\n            \n            no_na = common_df.copy()\n            no_na.dropna(axis='columns', inplace=True)\n            \n            for fill_method in fill_method_all:       \n                common_df_copy = common_df.copy()\n\n                if fill_method == 'mean':              \n                    common_df_copy[feature].fillna(common_df_copy.mean()[feature], inplace = True)  \n                elif fill_method == 'median': \n                    common_df_copy[feature].fillna(common_df_copy.median()[feature], inplace = True)  \n                elif fill_method == 'interpolation':\n                    common_df_copy[feature].fillna(common_df_copy.interpolate()[feature], inplace = True)  \n                    \n                X_train_feature = common_df_copy[common_df_copy['train'] == 1][feature]           \n                X_train_feature = pd.DataFrame(np.nan_to_num(X_train_feature), columns = {feature})\n                \n                scaler = StandardScaler()\n                if model_type == 'linear':\n                    scaler.fit(X_train_feature.values.reshape(-1, 1))\n                    X_train = scaler.transform(X_train_feature.values.reshape(-1, 1))  \n                elif model_type == 'rfr':\n                    X_train = X_train_feature.values.reshape(-1, 1)\n                \n                score = -np.mean(cross_val_score(model, \n                                            X_train, Y,      \n                                            cv = cv, \n                                            scoring = scoring))\n                mas_score.append(score)\n\n                if score < best_score:\n                    best_score = score \n                    best_method = fill_method\n                del common_df_copy\n\n            if best_method == 'mean':\n                common_df[feature].fillna(common_df.mean()[feature], inplace = True)\n            elif best_method == 'median': \n                common_df[feature].fillna(common_df.median()[feature], inplace = True)\n            elif best_method == 'interpolation': \n                common_df[feature].fillna(common_df.interpolate()[feature], inplace = True)\n            \n            del no_na\n            \n            print(f'Best score:     {best_score}')\n            \n        X = common_df.loc[common_df['train'] == 1,:]\n        X_test = common_df.loc[common_df['train'] == 0,:] \n        \n    else:\n        print('Zero features with missing values')\n        \n        X = common_df.loc[common_df['train'] == 1,:]\n        X_test = common_df.loc[common_df['train'] == 0,:]\n        \n    X.drop('train', axis = 1, inplace = True)\n    X_test.drop('train', axis = 1, inplace = True)\n    return [X, X_test, Y.reset_index(drop=True), lst_nan[1]]         ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Feature generation** <a class=\"anchor\" id=\"features\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/train.csv')\ndf_structure = pd.read_csv('../input/structures.csv')\ndf_test = pd.read_csv('../input/test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n# Feature generation from https://www.kaggle.com/kabure/simple-eda-lightgbm-autotuning-w-hyperopt/notebook\n\n## Function to reduce the DF size\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)    \n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df\n\ndef create_features(df):\n    df['molecule_couples'] = df.groupby('molecule_name')['id'].transform('count')\n    df['molecule_dist_mean'] = df.groupby('molecule_name')['dist'].transform('mean')\n    df['molecule_dist_min'] = df.groupby('molecule_name')['dist'].transform('min')\n    df['molecule_dist_max'] = df.groupby('molecule_name')['dist'].transform('max')\n    df['atom_0_couples_count'] = df.groupby(['molecule_name', 'atom_index_0'])['id'].transform('count')\n    df['atom_1_couples_count'] = df.groupby(['molecule_name', 'atom_index_1'])['id'].transform('count')\n    df[f'molecule_atom_index_0_x_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['x_1'].transform('std')\n    df[f'molecule_atom_index_0_y_1_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('mean')\n    df[f'molecule_atom_index_0_y_1_mean_diff'] = df[f'molecule_atom_index_0_y_1_mean'] - df['y_1']\n    df[f'molecule_atom_index_0_y_1_mean_div'] = df[f'molecule_atom_index_0_y_1_mean'] / df['y_1']\n    df[f'molecule_atom_index_0_y_1_max'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('max')\n    df[f'molecule_atom_index_0_y_1_max_diff'] = df[f'molecule_atom_index_0_y_1_max'] - df['y_1']\n    df[f'molecule_atom_index_0_y_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['y_1'].transform('std')\n    df[f'molecule_atom_index_0_z_1_std'] = df.groupby(['molecule_name', 'atom_index_0'])['z_1'].transform('std')\n    df[f'molecule_atom_index_0_dist_mean'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('mean')\n    df[f'molecule_atom_index_0_dist_mean_diff'] = df[f'molecule_atom_index_0_dist_mean'] - df['dist']\n    df[f'molecule_atom_index_0_dist_mean_div'] = df[f'molecule_atom_index_0_dist_mean'] / df['dist']\n    df[f'molecule_atom_index_0_dist_max'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('max')\n    df[f'molecule_atom_index_0_dist_max_diff'] = df[f'molecule_atom_index_0_dist_max'] - df['dist']\n    df[f'molecule_atom_index_0_dist_max_div'] = df[f'molecule_atom_index_0_dist_max'] / df['dist']\n    df[f'molecule_atom_index_0_dist_min'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('min')\n    df[f'molecule_atom_index_0_dist_min_diff'] = df[f'molecule_atom_index_0_dist_min'] - df['dist']\n    df[f'molecule_atom_index_0_dist_min_div'] = df[f'molecule_atom_index_0_dist_min'] / df['dist']\n    df[f'molecule_atom_index_0_dist_std'] = df.groupby(['molecule_name', 'atom_index_0'])['dist'].transform('std')\n    df[f'molecule_atom_index_0_dist_std_diff'] = df[f'molecule_atom_index_0_dist_std'] - df['dist']\n    df[f'molecule_atom_index_0_dist_std_div'] = df[f'molecule_atom_index_0_dist_std'] / df['dist']\n    df[f'molecule_atom_index_1_dist_mean'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('mean')\n    df[f'molecule_atom_index_1_dist_mean_diff'] = df[f'molecule_atom_index_1_dist_mean'] - df['dist']\n    df[f'molecule_atom_index_1_dist_mean_div'] = df[f'molecule_atom_index_1_dist_mean'] / df['dist']\n    df[f'molecule_atom_index_1_dist_max'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('max')\n    df[f'molecule_atom_index_1_dist_max_diff'] = df[f'molecule_atom_index_1_dist_max'] - df['dist']\n    df[f'molecule_atom_index_1_dist_max_div'] = df[f'molecule_atom_index_1_dist_max'] / df['dist']\n    df[f'molecule_atom_index_1_dist_min'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('min')\n    df[f'molecule_atom_index_1_dist_min_diff'] = df[f'molecule_atom_index_1_dist_min'] - df['dist']\n    df[f'molecule_atom_index_1_dist_min_div'] = df[f'molecule_atom_index_1_dist_min'] / df['dist']\n    df[f'molecule_atom_index_1_dist_std'] = df.groupby(['molecule_name', 'atom_index_1'])['dist'].transform('std')\n    df[f'molecule_atom_index_1_dist_std_diff'] = df[f'molecule_atom_index_1_dist_std'] - df['dist']\n    df[f'molecule_atom_index_1_dist_std_div'] = df[f'molecule_atom_index_1_dist_std'] / df['dist']\n    df[f'molecule_atom_1_dist_mean'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('mean')\n    df[f'molecule_atom_1_dist_min'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('min')\n    df[f'molecule_atom_1_dist_min_diff'] = df[f'molecule_atom_1_dist_min'] - df['dist']\n    df[f'molecule_atom_1_dist_min_div'] = df[f'molecule_atom_1_dist_min'] / df['dist']\n    df[f'molecule_atom_1_dist_std'] = df.groupby(['molecule_name', 'atom_1'])['dist'].transform('std')\n    df[f'molecule_atom_1_dist_std_diff'] = df[f'molecule_atom_1_dist_std'] - df['dist']\n    df[f'molecule_type_0_dist_std'] = df.groupby(['molecule_name', 'type_0'])['dist'].transform('std')\n    df[f'molecule_type_0_dist_std_diff'] = df[f'molecule_type_0_dist_std'] - df['dist']\n    df[f'molecule_type_dist_mean'] = df.groupby(['molecule_name', 'type'])['dist'].transform('mean')\n    df[f'molecule_type_dist_mean_diff'] = df[f'molecule_type_dist_mean'] - df['dist']\n    df[f'molecule_type_dist_mean_div'] = df[f'molecule_type_dist_mean'] / df['dist']\n    df[f'molecule_type_dist_max'] = df.groupby(['molecule_name', 'type'])['dist'].transform('max')\n    df[f'molecule_type_dist_min'] = df.groupby(['molecule_name', 'type'])['dist'].transform('min')\n    df[f'molecule_type_dist_std'] = df.groupby(['molecule_name', 'type'])['dist'].transform('std')\n    df[f'molecule_type_dist_std_diff'] = df[f'molecule_type_dist_std'] - df['dist']\n\n    df = reduce_mem_usage(df)\n    return df\n\ndef map_atom_info(df, atom_idx):\n    df = pd.merge(df, df_structure, 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\ndf_train = map_atom_info(df_train, 0)\ndf_train = map_atom_info(df_train, 1)\n\ndf_test = map_atom_info(df_test, 0)\ndf_test = map_atom_info(df_test, 1)\n\n## This is a very performative way to compute the distances\ntrain_p_0 = df_train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = df_train[['x_1', 'y_1', 'z_1']].values\ntest_p_0 = df_test[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = df_test[['x_1', 'y_1', 'z_1']].values\n\n## linalg.norm, explanation:\n## This function is able to return one of eight different matrix norms, \n## or one of an infinite number of vector norms (described below),\n## depending on the value of the ord parameter.\ndf_train['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)\ndf_test['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)\n\ndf_train['dist_x'] = (df_train['x_0'] - df_train['x_1']) ** 2\ndf_test['dist_x'] = (df_test['x_0'] - df_test['x_1']) ** 2\ndf_train['dist_y'] = (df_train['y_0'] - df_train['y_1']) ** 2\ndf_test['dist_y'] = (df_test['y_0'] - df_test['y_1']) ** 2\ndf_train['dist_z'] = (df_train['z_0'] - df_train['z_1']) ** 2\ndf_test['dist_z'] = (df_test['z_0'] - df_test['z_1']) ** 2\n\ndf_train['type_0'] = df_train['type'].apply(lambda x: x[0])\ndf_test['type_0'] = df_test['type'].apply(lambda x: x[0])\n\ndf_train = create_features(df_train)\ndf_test = create_features(df_test)\nprint(df_train.shape, df_test.shape)\n\nsubmission = pd.DataFrame({'id': df_test['id'],\n                          'scalar_coupling_constant' : [0] * len(df_test)})\n\ndf_train.drop(['id', 'molecule_name'], axis=1, inplace = True)\ndf_test.drop(['id', 'molecule_name'], axis=1, inplace = True)\nprint(df_train.shape, df_test.shape)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Dataset preprocessing** <a class=\"anchor\" id=\"dataprep\"></a>"},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in ['atom_1', 'atom_0', 'type']:\n        lbl = LabelEncoder()\n        lbl.fit(list(df_train[f].values) + list(df_test[f].values))\n        df_train[f] = lbl.transform(list(df_train[f].values))\n        df_test[f] = lbl.transform(list(df_test[f].values))\n\ndf_train['type_0'] = df_train['type_0'].astype(int)\ndf_test['type_0'] = df_test['type_0'].astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# frac = .1 to take 10% of the dataset to speed up the calculus \nfrac = .1\n\ncol = [c for c in df_train.columns if c not in ['scalar_coupling_constant']]\nX = df_train.sample(frac = frac)[col]\nY = df_train.loc[X.index, 'scalar_coupling_constant']\n\nX_TEST = df_test[col]\n\nX['train'] = 1\nX_TEST['train'] = 0\n\nprint(X.shape, X_TEST.shape, Y.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Feature selection, algorithm tuning, model training and prediction per molecule type** <a class=\"anchor\" id=\"preds\"></a>"},{"metadata":{"scrolled":false,"trusted":true},"cell_type":"code","source":"STATIC_PARAMS = {'metric': 'mae',\n                'n_estimators': 100,\n                'objective' : 'reg:squarederror',\n                'criterion' : 'mae',\n                }\n\nspace_SKopt = [Integer(2, 50, name='max_depth'),\n         Integer(2, 100, name='min_child_weight'),\n         Real(0.01, .1, name='learning_rate'),\n         Real(0.1, 1, name='subsample'),\n         Real(0.1, 1, name='colsample_bytree'),\n         Real(0.1, 1, name='reg_alpha'),\n         Real(0.1, 1, name='reg_lambda')]\n\nY_TEST_preds = pd.DataFrame({'ind': list(X_TEST.index), \n                             'type': X_TEST['type'].values, \n                             'prediction': [0] * len(X_TEST)})\n\nn_fold = 3\ncv = KFold(n_splits=n_fold, shuffle=True, random_state = 0)\n\nfor t in X['type'].unique():\n    moltype = lbl.inverse_transform(X.loc[X['type'] == t]['type'])[0]\n    print('='*45)\n    print(' '*12 + f'Working with type {moltype}')\n    print('='*45)\n    \n    print('Started at', ctime())\n        \n    mae, r2 = [], []\n    X_t = X.loc[X['type'] == t]\n    X_TEST_t = X_TEST.loc[X_TEST['type'] == t]\n    y_t = Y.loc[X['type'] == t]\n    \n    print('\\nFillna was started at', ctime())\n    fill_method_all = ['mean', 'median', 'interpolation']\n    cv_fillna = KFold(n_splits=3, shuffle=True, random_state = 0)\n    \n    # The next method deals with missing values in dataset considering three possible strategies from fill_method_all.\n    # The best strategy is chosen after MAE score evaluation for the linear model (random forest alg can be set using 'rfr' keyword instead of 'linear').\n    # Parameter N defines the percentage of observations with missing values to decide that the feature should be dropped or processed (fillna).  \n    # If you'd like to drop all features with any missing values, N must be set to null.\n    N = 0.0\n    [X_t, X_TEST_t, y_t, lst_proc] = smart_fillna(pd.concat([X_t, X_TEST_t], ignore_index = True), y_t, N, \n                                                                                                    fill_method_all, \n                                                                                                    'linear', \n                                                                                                    cv_fillna,\n                                                                                                    'neg_mean_absolute_error',)\n    print('\\nFeature selection was started at', ctime())\n    [X_t, X_TEST_t, lst_drop] = feature_selection(0.9, X_t, X_TEST_t)\n    print('Common features after fillna and FS:')\n    print(set(lst_proc)&set(lst_drop))\n    \n    oof = np.zeros(len(X_t))\n    X_p = X_t.sample(frac = 0.25)\n    Y_p = y_t.loc[X_p.index]\n    \n    model = xgboost.XGBRegressor(**STATIC_PARAMS)\n    [best_model,reg_gp] = get_params_SKopt(model, X_p, Y_p, space_SKopt)\n      \n    print ('Best score', reg_gp.fun)\n    print ('Best iterations', len(reg_gp.x_iters))\n    \n    best_model.n_jobs = multiprocessing.cpu_count()-1\n    best_model.n_estimators = 3000\n    print(best_model)\n    \n    for fold_n, (train_index, valid_index) in enumerate(cv.split(X_t)):\n        print('\\nFold', fold_n, 'started at', ctime())\n\n        X_train = X_t.iloc[train_index,:]\n        X_valid = X_t.iloc[valid_index,:]\n\n        Y_train = y_t.iloc[train_index]\n        Y_valid = y_t.iloc[valid_index]      \n        \n        best_model.fit(X_train, Y_train, \n               eval_metric = 'mae',    \n               eval_set = [(X_valid, Y_valid)],\n               verbose = False,\n               early_stopping_rounds = 10)\n\n        y_pred = best_model.predict(X_valid, \n                                   ntree_limit = best_model.best_iteration)\n\n        mae.append(mean_absolute_error(Y_valid, y_pred))\n        r2.append(r2_score(Y_valid, y_pred))\n        \n        print('Best score', best_model.best_score) \n        print('Best iteration', best_model.best_iteration)  \n            \n        Y_TEST_preds.loc[Y_TEST_preds['type'] == t, 'prediction'] += best_model.predict(X_TEST_t, \n                                                        ntree_limit = best_model.best_iteration)\n        oof[valid_index] = y_pred\n        \n    Y_TEST_preds.loc[Y_TEST_preds['type'] == t, 'prediction'] /= n_fold\n\n    print('='*45)\n    print('CV mean MAE: {0:.4f}, std: {1:.4f}.'.format(np.mean(mae), np.std(mae)))\n    print('CV mean R2:  {0:.4f}, std: {1:.4f}.'.format(np.mean(r2), np.std(r2)))\n    \n    plotfig(oof, y_t, 'Predicted vs. Actual responses for ' + moltype + ' molecules',max(y_t) + 0.1*max(y_t), reg_gp)   \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['scalar_coupling_constant'] = Y_TEST_preds['prediction'] \nprint(submission.head())\nsubmission.to_csv(\"submission_PMP.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.1"}},"nbformat":4,"nbformat_minor":1}