{"cells":[{"metadata":{},"cell_type":"markdown","source":"Here, we try to generate two meta-features for the second atom (atom_1), 1) mulliken charge and 2) mean value of tensor orientation contributions that have correlation with the target variable (scalar coupling constant). We use stuctural feature from Artgor's kernel https://www.kaggle.com/artgor/brute-force-feature-engineering and some periodical properties of the chemical elements, such as Pauling electronegativity, atomic, covalent, orbital and ion radiuses, ion charge, 1st ionization energy and electron affinity  "},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import numpy as np \nRANDOM_STATE = 12061985\nnp.random.seed(RANDOM_STATE)\n\nimport gc\nimport math\nimport xgboost\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport pandas as pd\npd.set_option(\"display.max_columns\", 999)\nimport seaborn as sns\nimport multiprocessing\nfrom time import time, ctime\nfrom scipy import linalg,stats\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm_notebook as tqdm\nfrom sklearn.metrics import mean_absolute_error, r2_score\nfrom sklearn.preprocessing import OneHotEncoder, LabelEncoder\nfrom sklearn.model_selection import KFold, ShuffleSplit, StratifiedShuffleSplit, cross_val_score\n\n# HPO\nfrom skopt.utils import use_named_args\nfrom skopt.plots import plot_convergence\nfrom skopt.space import Integer, Categorical, Real\nfrom skopt import gp_minimize, gbrt_minimize, forest_minimize\nfrom skopt.callbacks import EarlyStopper, DeltaXStopper, DeadlineStopper, DeltaYStopper","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plotfig (ypred, yactual, strtitle, y_max,y_min, 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([(y_min, y_min), (y_max, y_max)], [(y_min, y_min), (y_max, y_max)])\n    ax[0].set_xlim(y_min, y_max)\n    ax[0].set_ylim(y_min, 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 plotCorrelationMatrix(df, graphWidth):\n    df = df[[col for col in df if df[col].nunique() > 1]] # keep columns where there are more than 1 unique values\n    if df.shape[1] < 2:\n        print(f'No correlation plots shown: The number of non-NaN or constant columns ({df.shape[1]}) is less than 2')\n        return\n    corr = df.corr()\n    plt.figure(num=None, figsize=(graphWidth, graphWidth), dpi=80, facecolor='w', edgecolor='k')\n    \n    fmt = lambda x,pos: '{:.0%}'.format(x)\n    sns.heatmap(corr, square=True, annot=True, cmap='RdYlGn', annot_kws={\"size\": 10}, fmt='.1f')\n    \n    plt.xticks(range(len(corr.columns)), corr.columns, rotation=90)\n    plt.yticks(range(len(corr.columns)), corr.columns)\n    plt.gca().xaxis.tick_bottom()\n    plt.show()\n\ndef get_params_SKopt(model, X, Y, space, cv_search, opt_method = 'gbrt_minimize', verbose = True,  multi = False, scoring = 'neg_mean_squared_error', n_best = 50, total_time = 7200):\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                                        scoring= scoring))\n    if opt_method == 'gbrt_minimize':\n        \n        HPO_PARAMS = {'n_calls':1000,\n                      'n_random_starts':20,\n                      'acq_func':'EI',}\n        \n        reg_gp = gbrt_minimize(objective, \n                               space, \n                               n_jobs = -1,\n                               verbose = verbose,\n                               callback = [DeltaYStopper(delta = 0.01, n_best = 5), RepeatedMinStopper(n_best = n_best), DeadlineStopper(total_time = total_time)],\n                               **HPO_PARAMS,\n                               random_state = RANDOM_STATE)\n        \n    elif opt_method == 'forest_minimize':\n        \n        HPO_PARAMS = {'n_calls':1000,\n                      'n_random_starts':20,\n                      'acq_func':'EI',}\n        \n        reg_gp = forest_minimize(objective, \n                               space, \n                               n_jobs = -1,\n                               verbose = verbose,\n                               callback = [RepeatedMinStopper(n_best = n_best), DeadlineStopper(total_time = total_time)],\n                               **HPO_PARAMS,\n                               random_state = RANDOM_STATE)\n        \n    elif opt_method == 'gp_minimize':\n        \n        HPO_PARAMS = {'n_calls':1000,\n                      'n_random_starts':20,\n                      'acq_func':'gp_hedge',}        \n        \n        reg_gp = gp_minimize(objective, \n                               space, \n                               n_jobs = -1,\n                               verbose = verbose,\n                               callback = [RepeatedMinStopper(n_best = n_best), DeadlineStopper(total_time = total_time)],\n                               **HPO_PARAMS,\n                               random_state = RANDOM_STATE)\n    \n    TUNED_PARAMS = {} \n    for i, item in enumerate(space):\n        if multi:\n            TUNED_PARAMS[item.name.split('__')[1]] = reg_gp.x[i]\n        else:\n            TUNED_PARAMS[item.name] = reg_gp.x[i]\n    \n    return [TUNED_PARAMS,reg_gp]\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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's used just a fraction of the train and test datasets to speed up the calculus for the demonstration purpose. The processing of the whole datasets can't be commited during session time."},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nnrows = 10000\n\ndf_train = pd.read_csv('/kaggle/input/champs-scalar-coupling/train.csv', nrows = nrows)\ndf_test = pd.read_csv('/kaggle/input/champs-scalar-coupling/test.csv', nrows = nrows)\ndf_tensor = pd.read_csv('/kaggle/input/champs-scalar-coupling/magnetic_shielding_tensors.csv', nrows = nrows)\ndf_structure = pd.read_csv('/kaggle/input/champs-scalar-coupling/structures.csv')\ndf_mulliken = pd.read_csv('/kaggle/input/champs-scalar-coupling/mulliken_charges.csv')","execution_count":null,"outputs":[]},{"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    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    return df\n\ndef map_atom_info(df, atom_idx):\n    df = pd.merge(df, df_structure[['molecule_name', 'atom_index' ,'atom' ,'x' ,'y' ,'z']], 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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"corr_mas = df_train.iloc[:,5:].sample(frac = 0.25, random_state = RANDOM_STATE).corr()['scalar_coupling_constant']\nidx = corr_mas[np.abs(corr_mas) >= 0.1].index\n\ntrain = pd.concat([df_train.iloc[:,:5], df_train.loc[:,idx]], axis = 1)\ndel df_train\n\ntest = pd.concat([df_test.iloc[:,:5], df_test.loc[:,idx]], axis = 1)\ndel df_test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndf_structure = pd.merge(df_structure, df_tensor, how = 'left',\n                  left_on  = ['molecule_name', 'atom_index'],\n                  right_on = ['molecule_name', 'atom_index'])\n\ndf_structure = df_structure[['molecule_name', 'atom_index', 'atom', 'x',\n       'y', 'z', 'XX', 'XY', 'XZ', 'YX', 'YY', 'YZ', 'ZX', 'ZY', 'ZZ']]\nall_means_tensor = []\n\nfor i in tqdm(range(0, len(df_structure))):\n    if math.isnan(df_structure.loc[i, 'XX']):\n\n        all_means_tensor.append(np.nan)\n    else:\n        A = np.array([[df_structure.loc[i, 'XX'],df_structure.loc[i, 'XY'],df_structure.loc[i, 'XZ']],\n                      [df_structure.loc[i, 'YX'],df_structure.loc[i, 'YY'],df_structure.loc[i, 'YZ']],\n                      [df_structure.loc[i, 'ZX'],df_structure.loc[i, 'ZY'],df_structure.loc[i, 'ZZ']]])\n\n        w, v = linalg.eig(A)\n        all_means_tensor.append(np.real((w[0]+w[1]+w[2])/3))\n        \ndf_structure['k_tensor'] = all_means_tensor\ndf_structure.drop(['XX', 'XY', 'XZ', 'YX', 'YY', 'YZ','ZX','ZY', 'ZZ'], axis =1, inplace = True)\nprint(df_structure.shape)\ndf_structure.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_structure = pd.merge(df_structure, df_mulliken, how = 'left',\n                  left_on  = ['molecule_name', 'atom_index'],\n                  right_on = ['molecule_name', 'atom_index'])\ndf_structure.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"e_negativity_pol = {'H': 2.1,'C':2.55,'N':3.04,'O':3.44,'F':3.98}  #  Pauling electronegativity\ne_negativity_ol =  {'H': 2.1,'C':2.5,'N':3.07,'O':3.50,'F':4.1}  #  Allred–Rochow electronegativity\ncov_r     = {'H': 32, 'C': 77, 'N': 75, 'O': 73, 'F': 72} # Covalent radius\natomic_r  = {'H': 46, 'C': 77, 'N': 71, 'O': np.nan, 'F': np.nan} # Atomic radius\norbital_r = {'H': 53, 'C': 62, 'N': 52, 'O': 45, 'F': 40} # Orbital radius \nion_r     = {'H': 136, 'C': 255, 'N': 148, 'O': 136, 'F': 133} # ion radius\noxi_deg     = {'H': 1,   'C': 4,   'N': -3,  'O': -2,  'F' : -1} # ion charge\nfirst_ion_e = {'H': 1312,'C': 1086,'N': 1402,'O': 1314,'F' : 1681} # 1st ionization energy\ne_affinity = {'H': 0.754,'C':1.263,'N':-0.07,'O':1.461,'F':3.399} # Electron affinity\n\nq_e_ext = {'H': 1,'C': 4,'N': 5,'O': 6,'F' : 7}  #  Number of electrons on the external orbital\nq_levels = {'H': 1,'C': 2,'N': 2,'O': 2,'F' : 2}  #  Number of energy levels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_structure['e_negativity_pol'] = df_structure['atom'].apply(lambda x: e_negativity_pol[x]) \n\ndf_structure['cov_r'] = df_structure['atom'].apply(lambda x: cov_r[x]) \ndf_structure['atomic_r'] = df_structure['atom'].apply(lambda x: atomic_r[x]) \n\ndf_structure['orbital_r'] = df_structure['atom'].apply(lambda x: orbital_r[x]) \ndf_structure['ion_r'] = df_structure['atom'].apply(lambda x: ion_r[x]) \n\ndf_structure['oxi_deg'] = df_structure['atom'].apply(lambda x: oxi_deg[x]) \ndf_structure['first_ion_e'] = df_structure['atom'].apply(lambda x: first_ion_e[x]) \n\ndf_structure['e_affinity'] = df_structure['atom'].apply(lambda x: e_affinity[x]) \n\ndf_structure['q_e_ext'] = df_structure['atom'].apply(lambda x: q_e_ext[x]) \ndf_structure['q_levels'] = df_structure['atom'].apply(lambda x: q_levels[x]) \n\ndf_structure.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.merge(train, df_structure[['molecule_name', 'atom_index', 'mulliken_charge', 'k_tensor', 'e_negativity_pol', 'cov_r', 'atomic_r', 'orbital_r', 'ion_r',\n       'oxi_deg', 'first_ion_e', 'e_affinity', 'q_e_ext', 'q_levels']], how = 'left',\n                  left_on  = ['molecule_name','atom_index_0'],\n                  right_on = ['molecule_name','atom_index'])\n\ntrain = pd.merge(train, df_structure[['molecule_name', 'atom_index', 'mulliken_charge', 'k_tensor', 'e_negativity_pol', 'cov_r', 'atomic_r', 'orbital_r', 'ion_r',\n       'oxi_deg', 'first_ion_e', 'e_affinity', 'q_e_ext', 'q_levels']], how = 'left',\n                  left_on  = ['molecule_name','atom_index_1'],\n                  right_on = ['molecule_name','atom_index'])\n\ntest = pd.merge(test, df_structure[['molecule_name', 'atom_index', 'mulliken_charge', 'k_tensor', 'e_negativity_pol', 'cov_r', 'atomic_r', 'orbital_r', 'ion_r',\n       'oxi_deg', 'first_ion_e', 'e_affinity', 'q_e_ext', 'q_levels']], how = 'left',\n                  left_on  = ['molecule_name','atom_index_0'],\n                  right_on = ['molecule_name','atom_index'])\n\ntest = pd.merge(test, df_structure[['molecule_name', 'atom_index', 'mulliken_charge', 'k_tensor','e_negativity_pol', 'cov_r', 'atomic_r', 'orbital_r', 'ion_r',\n       'oxi_deg', 'first_ion_e', 'e_affinity', 'q_e_ext', 'q_levels']], how = 'left',\n                  left_on  = ['molecule_name','atom_index_1'],\n                  right_on = ['molecule_name','atom_index'])\n\ntest.drop(['mulliken_charge_x', 'mulliken_charge_y', 'k_tensor_x','k_tensor_y'], axis = 1, inplace = True)\n\ntest.drop(['id', 'molecule_name'], axis = 1, inplace = True)\ntrain.drop(['id', 'molecule_name'], axis = 1, inplace = True)\n\nprint(train.shape, test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plotCorrelationMatrix(train[['scalar_coupling_constant', 'mulliken_charge_x', 'mulliken_charge_y', 'k_tensor_x', 'k_tensor_y', 'e_negativity_pol_y', 'cov_r_y',\n                            'atomic_r_y','orbital_r_y', 'ion_r_y', 'oxi_deg_y', 'first_ion_e_y', 'e_affinity_y', 'q_e_ext_y', 'q_levels_y']],25)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in ['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))\n        \ntrain.fillna(0, inplace = True)\ntest.fillna(0, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_meta_feature_regression(X, X_TEST, Y):\n    Y_TEST_preds = pd.DataFrame({'ind': list(X_TEST.index), \n                             'prediction': [0] * len(X_TEST)})\n\n    STATIC_PARAMS = {'metric': 'mae',\n                    'n_estimators': 100,\n                    'objective' : 'reg:squarederror',\n                    'random_state' : RANDOM_STATE,\n                    'n_jobs': -1, \n                }\n\n    space_SKopt = [Integer(2, 80, name='max_depth'),\n                   Integer(2, 100, name='min_child_weight'),\n                   Real(0.01, .25, name='learning_rate'),\n                   Real(0.01, 1, name='subsample'),\n                   Real(0.1, 1, name='colsample_bytree'),\n                   Real(0.001, 5, name='reg_alpha'),\n                   Real(0.001, 5, name='reg_lambda')\n                  ]\n    mae, r2 = [],[]\n    oof = np.zeros(len(X))\n\n#     check_features = ['molecule_atom_index_0_y_1_mean_div']\n#     for feature in check_features:\n#         train_test = pd.concat([X[feature], X_TEST[feature]])\n#         limits_IQR = [5, 95]\n#         q99 = np.percentile(train_test, 95)\n#         q01 = np.percentile(train_test, 5)\n\n#         X[feature] = np.clip(X[feature], q01, q99)\n#         X_TEST[feature] = np.clip(X_TEST[feature],q01, q99)\n#         del train_test\n    \n    X_p = X.sample(frac = 0.25, random_state = RANDOM_STATE)\n    Y_p = Y.loc[X_p.index]\n    \n    start_time = time() \n    cv_tune = ShuffleSplit(n_splits=1, test_size = 0.3, random_state = RANDOM_STATE)\n    [TUNED_PARAMS,reg_gp] = get_params_SKopt(xgboost.XGBRegressor(**STATIC_PARAMS), \n                                                             X_p, Y_p, \n                                                             space_SKopt, \n                                                             cv_tune,\n                                                             opt_method = 'forest_minimize',\n                                                             verbose = False,\n                                                             multi = False, \n                                                             scoring = 'neg_mean_absolute_error',\n                                                             n_best = 10,\n                                                             total_time = 7200)\n\n    print('\\nTime for tuning: {0:.2f} minutes'.format((time() - start_time)/60))\n    NEW_PARAMS = {**STATIC_PARAMS, **TUNED_PARAMS}\n    best_model = xgboost.XGBRegressor(**NEW_PARAMS)\n    \n    print ('Best score', reg_gp.fun)\n    print ('Best iterations', len(reg_gp.x_iters))\n    best_model.n_estimators = 2000\n    print(best_model)   \n    \n    n_fold = 3\n    cv = KFold(n_splits=n_fold, shuffle=True, random_state = 0)\n    for fold_n, (train_index, valid_index) in enumerate(cv.split(X)):\n        print('\\nFold', fold_n, 'started at', ctime())\n\n        X_train = X.iloc[train_index,:]\n        X_valid = X.iloc[valid_index,:]\n\n        Y_train = Y.iloc[train_index]\n        Y_valid = Y.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['prediction'] += best_model.predict(X_TEST, \n                                                        ntree_limit = best_model.best_iteration)\n        oof[valid_index] = y_pred\n\n    Y_TEST_preds['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, 'Predicted vs. Actual responses',max(Y) + 0.1*max(Y), \n            min(Y) - 0.1*min(Y), reg_gp)\n    return Y_TEST_preds['prediction']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.drop(['e_negativity_pol_x', 'cov_r_x', 'atomic_r_x', 'orbital_r_x', 'ion_r_x', 'oxi_deg_x', 'first_ion_e_x', 'e_affinity_x', 'q_e_ext_x', 'q_levels_x', 'mulliken_charge_x', 'k_tensor_x'], axis = 1, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.drop(['e_negativity_pol_x', 'cov_r_x', 'atomic_r_x', 'orbital_r_x', 'ion_r_x', 'oxi_deg_x', 'first_ion_e_x', 'e_affinity_x', 'q_e_ext_x', 'q_levels_x'], axis = 1, inplace = True)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"res_test = pd.DataFrame()\nres_train = pd.DataFrame()\n\nall_notin = ['scalar_coupling_constant', 'mulliken_charge_y', 'k_tensor_y']\n\nfor meta_feature in ['mulliken_charge_y', 'k_tensor_y']: \n    print('\\n'+'*'*70 + '\\n' + ' '*12 + f'Processing {meta_feature}' + '\\n' +'*'*70)\n    \n    col = [c for c in train.columns if c not in all_notin]     \n    preds_meta = calc_meta_feature_regression(train[col], test[col], train[meta_feature])\n    \n    test[meta_feature] = preds_meta\n    res_test[meta_feature] = preds_meta\n    res_train[meta_feature] = train[meta_feature]\n    all_notin.remove(meta_feature)\n    \nres_test.to_csv('mulliken_for_test.csv', index = False)\nres_train.to_csv('mulliken_for_train.csv', index = False)","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.3"}},"nbformat":4,"nbformat_minor":1}