{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel contains feature engineering using the molecule structures file and subsequent training of several lightgbm estimators and prediction.\nThe training set is split because training all at once has led to out of memory.\nThe most recent addition are the openbabel features."},{"metadata":{"trusted":true},"cell_type":"code","source":"!conda install -c openbabel openbabel -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Do the standard imports and have a look at the data \nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\n\nprint(os.listdir(\"../input\"))\n\n\"\"\"\nReduce Mem usage function taken from here:\nhttps://www.kaggle.com/artgor/artgor-utils\n\"\"\"\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)\n            else:\n                c_prec = df[col].apply(lambda x: np.finfo(x).precision).max()\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max and c_prec == np.finfo(np.float16).precision:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max and c_prec == np.finfo(np.float32).precision:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n                    \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\ndf_struct = pd.read_csv('../input/structures.csv')\n\n# Any results you write to the current directory are saved as output.\n\ndf_struct = reduce_mem_usage(df_struct)\n\nprint('The largest molecule has {m} atoms'.format(m=str(df_struct['atom_index'].max())))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_struct['atom'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# df = pd.read_csv('../input/train.csv')\n# df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import openbabel\nobConversion = openbabel.OBConversion()\nobConversion.SetInFormat(\"xyz\")\nxyz_path = '../input/structures'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !ls ../input/structures","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This kernel generates some Features and writes to files for further processing.\n- Distance between the atoms\n- OHE encoded connection type and and atom types\n- Distance to all other atoms in the molecule and the cosine of the angle w.r.t. the vector between the atoms\n\nTraining is done with lightgbm with custom loss and objective. The data split has to be split in three equal parts to prevent memory exhaustion. This drags on the result.\n\nFinally Feature importance and a graphical representation of the results are shown."},{"metadata":{},"cell_type":"markdown","source":"the feature engineering is done with matrix style operations. While it is not as flexible as going line by line it is far more performant.\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"import progressbar as pb\n\n# there are 29 molecules in the largest molecule\nn = 29\ndf_molecules = pd.DataFrame(df_struct['molecule_name'].unique())\ndf_molecules.columns = ['molecule_name']\n# next part is to retrieve some features via the openbabel module\n\ndef enrich_row(row):\n    mol = openbabel.OBMol()\n    obConversion.ReadFile(mol, '{xyz_path}/{molecule_name}.xyz'.format(xyz_path=xyz_path, molecule_name=row['molecule_name']))\n    row['num_atoms'] = mol.NumAtoms()\n    row['num_bonds'] = mol.NumBonds()\n    row['num_hvy'] = mol.NumHvyAtoms()\n    row['num_residues'] = mol.NumResidues()\n    row['num_rotors']  = mol.NumRotors()\n    row['mol_energy'] = mol.GetEnergy()\n    row['mol_wt'] = mol.GetMolWt()\n    row['mol_mass'] = mol.GetExactMass()\n    row['mol_charge'] = mol.GetTotalCharge()\n    row['mol_spin'] = mol.GetTotalSpinMultiplicity()\n    row['mol_dimension'] = mol.GetDimension()\n    return row\n    \ndf_molecules = df_molecules.apply(enrich_row, axis=1)\n\ndf_molecules.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(n):\n    rsuff = '_{idx}'.format(idx=str(i))\n    df_molecules = df_molecules.merge(df_struct[df_struct['atom_index'] == i], how='left',\n                                     on='molecule_name',\n                                     suffixes=('', rsuff))\n    drop_cols = ['atom_index']\n    df_molecules = df_molecules.drop(drop_cols, axis=1)\n\n\n# rename the first columns for ease of use later on\ndf_molecules = df_molecules.rename(index=str, columns={'atom': 'atom_0', 'x': 'x_0', 'y': 'y_0', 'z': 'z_0'})\ndf_molecules.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(df_molecules.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next step is to do the actual feature engineering. The steps are as follows:\n- Join every line with its corresponding molecule structure\n- Rotate all the (x,y,z) cartesian coordinates such that atom 1 is at (0,0,0) and atom 2 is at (x,0,0)\n- OHE the type and all atom types"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"df_train = pd.read_csv('../input/train.csv')\n\nprint(df_train.shape)\natom_list = list(df_struct['atom'].unique())\ntype_list = list(df_train['type'].unique())\n\n# save relative freq for later use as weights. This is meant to improve adherence to the weighted objective\n\nweight_dict = dict()\nfor typ in type_list:\n    weight_dict[typ] = float(df_train.shape[0]) / float(df_train[df_train['type'] == typ].shape[0])\n\nn_train = df_train.shape[0]\ndef get_features(df_train, df_struct, atom_list, type_list, n=29):\n    df_train = df_train.merge(df_struct, how='left',\n                              left_on=['molecule_name', 'atom_index_0'],\n                              right_on=['molecule_name', 'atom_index']).drop('atom_index', axis=1)\n    df_train = df_train.rename(index=str, columns={'atom': 'atom_l', 'x': 'x_l', 'y': 'y_l', 'z': 'z_l'})\n    df_train = df_train.merge(df_struct, how='left',\n                              left_on=['molecule_name', 'atom_index_1'],\n                              right_on=['molecule_name', 'atom_index']).drop('atom_index', axis=1)\n    df_train = df_train.rename(index=str, columns={'atom': 'atom_r', 'x': 'x_r', 'y': 'y_r', 'z': 'z_r'})\n    df_train = df_train.merge(df_molecules, how='left', on='molecule_name')\n    # make OHE\n    # for typ in type_list:\n    #     colname = 'type_' + typ\n    #     df_train[colname] = df_train['type'] == typ\n    # df_train = df_train.drop('type', axis=1)\n    cat_cols = ['atom_l', 'atom_r']\n    for atom in atom_list:\n        df_train['count_{atom}'.format(atom=atom)] = 0\n    for i in range(n):\n        cat_cols.append('atom_{idx}'.format(idx=str(i)))\n    for cat_col in cat_cols:\n        for atom in atom_list:\n            colname = cat_col + '_' + atom\n            df_train[colname] = df_train[cat_col] == atom\n            df_train['count_{atom}'.format(atom=atom)] += df_train[cat_col] == atom\n        df_train = df_train.drop(cat_col, axis=1)\n    # make the atom_index_0 atom the (0,0,0) in the coordinate system\n    for dimcol in ['x', 'y', 'z']:\n        for lridx in ['l', 'r']:\n            tgtcol = dimcol+'_'+lridx\n            srccol = dimcol+'_l'\n            df_train[tgtcol] = df_train[tgtcol] - df_train[srccol]\n    # now we can drop the *_r columns again as they are all 0\n    df_train = df_train.drop(['x_l', 'y_l', 'z_l'], axis=1)\n    # now rotate everything such that atom_index_1 is at (x, 0, 0)\n    # first rotate around (0, 0, 1)\n    df_train['dist'] = np.sqrt(df_train['y_r']*df_train['y_r']+df_train['x_r']*df_train['x_r'])\n    df_train['sintheta'] = -df_train['y_r']/df_train['dist']\n    df_train['costheta'] = df_train['x_r']/df_train['dist']\n    for i in range(n):\n        x_colname = 'x_' + str(i)\n        y_colname = 'y_' + str(i)\n        df_train['x_tmp'] = df_train[x_colname]\n        df_train[x_colname] = df_train[x_colname]*df_train['costheta']-df_train[y_colname]*df_train['sintheta']\n        df_train[y_colname] = df_train['x_tmp']*df_train['sintheta']+df_train[y_colname]*df_train['costheta']\n    x_colname = 'x_r'\n    y_colname = 'y_r'\n    df_train['x_tmp'] = df_train[x_colname]\n    df_train[x_colname] = df_train[x_colname]*df_train['costheta']-df_train[y_colname]*df_train['sintheta']\n    df_train[y_colname] = df_train['x_tmp']*df_train['sintheta']+df_train[y_colname]*df_train['costheta']\n    # now rotate around (0, 1, 0)\n    df_train['dist'] = np.sqrt(df_train['z_r']*df_train['z_r']+df_train['x_r']*df_train['x_r'])\n    df_train['sintheta'] = -df_train['z_r']/df_train['dist']\n    df_train['costheta'] = df_train['x_r']/df_train['dist']\n    for i in range(n):\n        x_colname = 'x_' + str(i)\n        z_colname = 'z_' + str(i)\n        df_train['x_tmp'] = df_train[x_colname]\n        df_train[x_colname] = df_train[x_colname]*df_train['costheta']-df_train[z_colname]*df_train['sintheta']\n        df_train[z_colname] = df_train['x_tmp']*df_train['sintheta']+df_train[z_colname]*df_train['costheta']\n    x_cols = []\n    y_cols = []\n    z_cols = []\n    for i in range(n):\n        x_cols.append('x_{num}'.format(num=str(i)))\n        y_cols.append('y_{num}'.format(num=str(i)))\n        z_cols.append('z_{num}'.format(num=str(i)))        \n    df_train['x_min'] = df_train[x_cols].min(axis=1)\n    df_train['x_max'] = df_train[x_cols].max(axis=1)\n    df_train['x_std'] = df_train[x_cols].std(axis=1)\n    df_train['y_min'] = df_train[y_cols].min(axis=1)\n    df_train['y_max'] = df_train[y_cols].max(axis=1)\n    df_train['y_std'] = df_train[y_cols].std(axis=1)\n    df_train['z_min'] = df_train[z_cols].min(axis=1)\n    df_train['z_max'] = df_train[z_cols].max(axis=1)\n    df_train['z_std'] = df_train[z_cols].std(axis=1)\n    \n    for i in range(n):\n        dist_colname = 'distance_' + str(i)\n        x_colname = 'x_' + str(i)\n        y_colname = 'y_' + str(i)\n        z_colname = 'z_' + str(i)\n        df_train[dist_colname] = df_train[x_colname]*df_train[x_colname]+df_train[y_colname]*df_train[y_colname]+df_train[z_colname]*df_train[z_colname]\n        for spatial_colname in [x_colname, y_colname, z_colname]:\n            df_train[spatial_colname] = df_train[spatial_colname]/df_train[dist_colname]\n    x_colname = 'x_r'\n    z_colname = 'z_r'\n    df_train['x_tmp'] = df_train[x_colname]\n    df_train[x_colname] = df_train[x_colname]*df_train['costheta']-df_train[z_colname]*df_train['sintheta']\n    df_train[z_colname] = df_train['x_tmp']*df_train['sintheta']+df_train[z_colname]*df_train['costheta']\n    # lets drop everything that is useless and lets give it a try\n    df_train = df_train.drop(['x_tmp',\n                              'y_r',\n                              'z_r',\n                              'sintheta',\n                              'costheta',\n                              'dist',\n                              'molecule_name',\n                              'id',\n                              'atom_index_0',\n                              'atom_index_1'], axis=1)\n    return df_train\n\n# df_train = get_features(df_train, df_struct, atom_list, type_list)\n# print(df_train.shape)\n# df_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The ensemble class helps training multiple estimators on different parts of the data. This is a workaround to prevent memory exhaust."},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndel df_train\n\nimport lightgbm as lgb\nfrom sklearn.model_selection import train_test_split\n\n\nclass lgb_ensemble(object):\n    def __init__(self):\n        self.estimators = []\n        \n    def custom_loss(self, y_hat, d_train, type_list, x_valid):\n        \"\"\"\n        Is used to compute the loss in the validation set. To be honest i could not get it to work properly\n        \"\"\"\n        y = d_train.get_label()\n        score = 0\n        n_score = 0\n        for typ in type_list:\n            colname = 'type_' + typ\n            normalization = np.sum(x_valid[colname])\n            if normalization > 0:\n                n_score += 1.0\n                score += np.log(np.sum(np.multiply(x_valid[colname], np.abs(y-y_hat)))/normalization)\n        return 'Grouped Log MAE', score/n_score, False\n    \n    def custom_objective(self, y_hat, d_train, type_list, x_train):\n        \"\"\"\n        Is used to calculate the gradient and diagonal elements of hessian on training set\n        \"\"\"\n        T = y_hat.shape[0]\n        y = d_train.get_label()\n        grad = np.zeros(y.shape[0], )\n        hessian = np.zeros(y.shape[0], )\n        for typ in type_list:\n            colname = 'type_' + typ\n            nt = np.sum(x_train[colname])\n            delta_vec = np.multiply(x_train[colname], y-y_hat)\n            grad += np.multiply(x_train[colname], delta_vec*T/nt)\n            hessian += -np.multiply(x_train[colname], np.sign(delta_vec)*((delta_vec**2)*T/nt))\n        return grad, hessian\n            \n            \n    def train(self, df_train, target, params, typ):\n        x_train, x_valid, y_train, y_valid = train_test_split(df_train, target, test_size=0.15)\n        d_train = lgb.Dataset(x_train, label=y_train)\n        d_valid = lgb.Dataset(x_valid, label=y_valid)\n        # The type_list is defined outside if this class, using it here again can be considered sloppy. However i dont care, call the police if you need to\n        # cust_loss = lambda y_hat, d_train: self.custom_loss(y_hat, d_train, type_list, x_valid)\n        # cust_objective = lambda y_hat, d_train: self.custom_objective(y_hat, d_train, type_list, x_train)\n            \n        n_estimators = 15000\n        watchlist = [d_valid]\n        estim = lgb.train(params,\n                            d_train,\n                            n_estimators,\n                            watchlist,\n                            # fobj=cust_objective,\n                            # feval=cust_loss,\n                            verbose_eval=1000,\n                            early_stopping_rounds=1000)\n        self.estimators.append({'estimator': estim, 'numdata': df_train.shape[0], 'typ': typ})\n        return list(y_valid), list(estim.predict(x_valid))\n    def predict(self, features):\n        y_pred = np.zeros((features.shape[0], ), dtype=np.float32)\n        types = features['type']\n        features = features.drop('type', axis=1)\n        normalization = 0\n        for estim in self.estimators:\n            y_pred[types == estim['typ']] = estim['estimator'].predict(features[types == estim['typ']])\n        return y_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# del df_train\nparams = {\n                    # 'min_data_in_leaf': 20,\n                    'boosting_type': 'gbdt',\n                    'objective': 'regression',\n                    'metric': 'mae', # the actual target is log(mae) weighted per type\n                    'max_depth': 15,\n                    'num_leaves': 120,\n                    'learning_rate': 0.05,\n                    'feature_fraction': 0.9,\n                    'verbose': 1\n                    # 'bagging_fraction': 0.5,\n                    # 'bagging_freq': 10,\n                }\n\n\"\"\"\nThe next step is to train one model for each type there is\n\"\"\"\n\n# chunksize = int(n_train/3)\n# df_train_all = pd.read_csv('../input/train.csv', chunksize=chunksize)\nestimator = lgb_ensemble()\ny_valid_all = []\ny_pred_all = []\nmax_len = 1_200_000\nfor typ in type_list:\n    df_train = pd.read_csv('../input/train.csv')\n    df_train = df_train[df_train['type'] == typ]\n    if df_train.shape[0] > max_len:\n        df_train = df_train.iloc[:max_len, :]\n    df_train = reduce_mem_usage(df_train)\n# for df_train in df_train_all:\n    # weights = [weight_dict[row['type']] for i, row in df_train.iterrows()]\n    df_train = get_features(df_train, df_struct, atom_list, type_list)\n    target = df_train['scalar_coupling_constant']\n    df_train = df_train.drop(['scalar_coupling_constant', 'type'], axis=1)\n    df_train.head()\n    y_valid, y_pred = estimator.train(df_train, target, params, typ)\n    for i in range(len(y_valid)):\n        y_valid_all.append(y_valid[i])\n        y_pred_all.append(y_pred[i])\n\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now apply the model to test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# apply the model to test data\ndel df_train\ndf_test_all = pd.read_csv('../input/test.csv', chunksize=1000)\nisfirst = True\n\nfor df_test in df_test_all:\n    df_test = reduce_mem_usage(df_test)\n    df_test_orig = df_test.copy()\n    df_test = get_features(df_test, df_struct, atom_list, type_list)\n    df_test_orig['scalar_coupling_constant'] = estimator.predict(df_test)\n    if isfirst:\n        df_test_orig.filter(['id', 'scalar_coupling_constant']).to_csv('prediction.csv', index=False, mode='w')\n        isfirst = False\n    else:\n        df_test_orig.filter(['id', 'scalar_coupling_constant']).to_csv('prediction.csv', header=False, index=False, mode='a')\n    \nprint('Written to disk')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Do some visualizations of feature importance and model performance"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport warnings\n%matplotlib inline\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\nfeature_importances = estimator.estimators[0]['estimator'].feature_importance(importance_type='gain')\n\nfor i in range(1, len(estimator.estimators)):\n    feature_importances += estimator.estimators[i]['estimator'].feature_importance(importance_type='gain')\n\nfeature_imp = pd.DataFrame(sorted(zip(feature_importances,df_train.columns)), columns=['Value','Feature'])\nfeature_imp['Value'] = feature_imp['Value']/feature_imp['Value'].sum() # normalize to 1\nsorted_values = feature_imp.sort_values(by=\"Value\", ascending=False)\nsorted_values['cum'] = np.log(sorted_values['Value'].cumsum()/sorted_values['Value'].sum())\nsorted_values = sorted_values[:20]\n\n# save for later use\nimportant_features = sorted_values['Feature'].copy()\n\nplt.figure(figsize=(20, 10))\nsns.barplot(x=\"Value\", y=\"Feature\", data=sorted_values)\nplt.title('LightGBM Feature Importance')\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot predicted vs actual valid y's\n\nprint(mpl.rcParams['agg.path.chunksize'])\nmpl.rcParams['agg.path.chunksize'] = 10000\nplt.figure(figsize=(20, 10))\nplt.plot(y_valid_all, y_pred_all, '.')\nplt.title('Prediction Plot')\nplt.tight_layout()\nplt.show()","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}