{"cells":[{"metadata":{},"cell_type":"markdown","source":"I decided to be part of this competition because proteins' structure and conformation intrigues me. Although I imagine that many wonderful kernels were generated with little knowledge on the topic, I decided to take the opportunity to learn about it.\n\nI tried to study it using books, but couldn't begin to understand. This amazing video course available at Khan Academy have helped me better understand how the data is generated. I totally suggest it as an information source: https://www.khanacademy.org/science/organic-chemistry/spectroscopy-jay. From now on I'll call it Khan to reference purposes."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nfrom numpy.random import permutation\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport lightgbm as lgb\n\nfrom sklearn import metrics\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import LabelEncoder\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nos.listdir(\"../input\")\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv(f'../input/train.csv')\ntest = pd.read_csv(f'../input/test.csv')\nstructures = pd.read_csv(f'../input/structures.csv')\n\ntrain[train['molecule_name'] == \"dsgdb9nsd_000001\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### **Generating the Euclidean distances¶**\n\nAs the test set is limited to spatial coordinates, I started with the classical Euclidean distance used in almost all kernels.\n\nTo do so, I merged each atom info from the structures_df to the train and test dataframes.\nAndrew Lukyanenko's kernel was my reference in this early stage of the process."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Adapted from the Andrew Lukyanenko's kernel. \n# https://www.kaggle.com/artgor/molecular-properties-eda-and-models\ndef map_atom_info(df, df2, atom_idx, col_names):\n    df = pd.merge(df, df2, 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    for i in range(len(col_names)):\n        df = df.rename(columns={f'{col_names[i]}': f'{col_names[i]}_{atom_idx}'})\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_change = ['x','y','z','atom']\ntrain = map_atom_info(train, structures, 0, cols_change)\ntrain = map_atom_info(train, structures, 1, cols_change)\n\ntest = map_atom_info(test, structures, 0, cols_change)\ntest = map_atom_info(test, structures, 1, cols_change)\n\nlist(test.columns)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The merging proccess worked well. Let's calculate Euclidean Distances between the atoms as Andrew did."},{"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":"We decided to labelEncode our data rigth now to able to compare data before and after predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# LabelEncoding the the character values, splitting and reducing the dataframe sizes\nfor f in ['atom_0', 'atom_1', 'type']:\n    lbl = LabelEncoder()\n    lbl.fit(list(train[f].values) + list(test[f].values))\n    lbl_name_mapping = dict(zip(lbl.classes_, lbl.transform(lbl.classes_)))\n    train[f] = lbl.transform(list(train[f].values))\n    test[f] = lbl.transform(list(test[f].values))\n    print(lbl_name_mapping)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Khan Academy course doesn't explain how to understand the types of J-coupling presented in both datasets, but from the begging (and the Organic Chemistry books I found), they seemed to be a strong influence on the scalar coupling values.\n\nWhat I understood is that they give us a clue as to which atoms the hydrogens (protons) are closest to and, as O, N and F have a higher electronegativity, they de-shield the proton and expose it to different levels of magnetic fielding (If anyone who knows this better hits this kernel, please correct me if I explained it wrong).\n\nAnyways, I decided to use it somehow.\n\nThese are the available types and some basic statistics about scalar coupling:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.groupby('type')['scalar_coupling_constant'].agg(\n    ['count', 'std', 'mean', 'median', 'var', 'skew']).reset_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Type 1JHC (code 0), which the following site (http://sopnmr.ucsd.edu/coupling.htm) describes as \"Peptide one bond\" had the highest variance and standard deviation. If I understand correctly, these are direct C-H bonds and they can vary greatly due to other molecular bonds.\n\nThe second biggest variance was 1JHN (code 1). The second biggest variance was 1JHN. I used this awsome kernel to understand a bit more the molecule possible structures: https://www.kaggle.com/cdeotte/molecule-animations\nSince N requires 3 bonds, I played a little with this kernel and tried to validate some of them. I have found that N is sometimes bonded to a C or an O and, to make 3 bonds, sometimes a double bond with some of these atoms is required. Khan explained that this could also be a reason for proton de-shielding.\n\nAs our distances can describe the different patterns that cause the large range of scalar coupling in these cases, I decided to run the basic model of distances before generating any other feature."},{"metadata":{},"cell_type":"markdown","source":"### **Auxiliary functions for automating model generation**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# This memory reduction function will make the execution faster\n# https://www.kaggle.com/artgor/artgor-utils\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# evaluation metric for validation\n# https://www.kaggle.com/abhishek/competition-metric\ndef metric(df, preds):\n    df[\"prediction\"] = preds\n    maes = []\n    for t in df.type.unique():\n        y_true = df[df.type==t].scalar_coupling_constant.values \n        y_pred = df[df.type==t].prediction.values\n        mae = np.log(metrics.mean_absolute_error(y_true, y_pred))\n        maes.append(mae)\n    return np.mean(maes)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Splitting the \"train set\" in two to allow internal model validation (adapted from)\n# https://www.kaggle.com/robertburbidge/using-estimated-mulliken-charges\ndef split_data(df):\n    \n    molecule_names = pd.DataFrame(permutation(df['molecule_name'].unique()),columns=['molecule_name'])\n    nm = molecule_names.shape[0]\n    ntrn = int(0.9*nm)\n    nval = int(0.1*nm)\n\n    tmp_train = pd.merge(df, molecule_names[0:ntrn], how='right', on='molecule_name')\n    tmp_train = reduce_mem_usage(tmp_train)\n    tmp_val = pd.merge(df, molecule_names[ntrn:nm], how='right', on='molecule_name')\n    tmp_val = reduce_mem_usage(tmp_val)\n\n    X_train = tmp_train.drop(['id', 'molecule_name', 'scalar_coupling_constant'], axis=1)\n    y_train = tmp_train['scalar_coupling_constant']\n\n    X_val = tmp_val.drop(['id', 'molecule_name', 'scalar_coupling_constant'], axis=1)\n    y_val = tmp_val['scalar_coupling_constant']\n\n    return X_train, y_train, X_val, y_val","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I'd love to do some tuning and cross-validation, but I'll leave it for the next competition."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, y_train, X_val, y_val = split_data(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'num_leaves': 128,\n          'min_child_samples': 79,\n          'objective': 'regression_l1',\n          'max_depth': 13,\n          'learning_rate': 0.2,\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'\n          : 1.0\n         }\n\nmodel_basic = lgb.LGBMRegressor(**params, n_estimators = 5000, n_jobs = -1)\nmodel_basic.fit(X_train, y_train, \n                    eval_set=[(X_train, y_train), (X_val, y_val)],\n                    verbose=500, early_stopping_rounds=100)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the automatic results of the model we can see overfitting. This may be addressed with some parameters tuning.\n\nI'll predict over the train data to so I'll be able to generate the same statistics by type to see if the distances had enough pattern to identify the different scenarios in 1JHC and 1JHN."},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_train = model_basic.predict(X_train)\npred_train_median = np.full(y_train.shape, np.median(y_train))\nprint(metrics.mean_absolute_error(y_train, pred_train) / metrics.mean_absolute_error(y_train, pred_train_median))\nmetric(pd.concat([X_train, y_train], axis=1), pred_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_train_df = pd.DataFrame(pred_train, columns=['pred_sc'])\nX_train_types = pd.concat([X_train.reset_index(drop=True), pred_train_df], axis=1)\nX_train_types.groupby('type')['pred_sc'].agg(\n    ['count', 'std', 'mean', 'median', 'var', 'skew']).reset_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we couldn't see any significative variance change in the train dataset for the real and the predicted scalar coupling values, we can assume the distances themselves probably justify most of the scalar coupling range differences in the first tow types. The third one (2JHC - code 2), presented a representative variation. \n\nThey surelly don't explain all the influence, so I decided to check the regression plot to see if some types are more out of phase than others."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_types = pd.concat([X_train_types.reset_index(drop=True), y_train.rename('original_sc')], axis=1)\nsns.lmplot(data=X_train_types, x=\"pred_sc\", y=\"original_sc\", hue=\"type\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And the 1JHC seem more disperse in our predictions, as the type 3JHC (code 5), which may confirm the Khan course discussions of electronegativity and de-sheilding. \nThese may be a topic for another Kernel."},{"metadata":{},"cell_type":"markdown","source":"### **Applying the model for the validation and test data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_val = model_basic.predict(X_val)\npred_val_median = np.full(y_val.shape, np.median(y_val))\nprint(metrics.mean_absolute_error(y_val, pred_val) / metrics.mean_absolute_error(y_val, pred_val_median))\nmetric(pd.concat([X_val, y_val], axis=1), pred_val)\n\npred_val_df = pd.DataFrame(pred_val, columns=['pred_sc'])\nX_val_types = pd.concat([X_val.reset_index(drop=True), pred_val_df], axis=1)\nX_val_types = pd.concat([X_val_types.reset_index(drop=True), y_val.rename('original_sc')], axis=1)\nsns.lmplot(data=X_val_types, x=\"pred_sc\", y=\"original_sc\", hue=\"type\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The validation data showed some JHC bond dispersion as the train did, which indicate that we don't have bias in the data division. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# We can't forget the memory usage reduction\ntest = reduce_mem_usage(test)\nX_test = test.drop(['id', 'molecule_name'], axis=1)\n\npred_test = model_basic.predict(X_test)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"sub = pd.read_csv(f'../input/sample_submission.csv')\nsub['scalar_coupling_constant'] = pred_test\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Conclusions: ** The distances between the atoms already explain a lot of their scalar coupling values. This aspect is in agreement with what the Khan course discusses on bondings and the impact caused by atoms in the proton's vicinity. The 0.25 for mean absolute error isn't awful, but a lot of transformations are possible given the provided data. \n\n**Additional features:** With the 2 weeks I took to understand the topic, study the available kernels and try to implement my own, I finish this work with a lot of transformation ideas. Some of them are:\n- Chris Deotte's Molecule Visualization kernel can be used to calculate which atoms are near each H and we can apply some degrees of electronegativity to each of them (https://www.kaggle.com/cdeotte/molecule-animations);\n- Binding types (single, double or triple) can be estimated using the expected amounts of bonds for each atom and degrees of de-shielding can be applied to the estimated potential energy of the molecule. (calculated in the Robert Burbidge kernel - https://www.kaggle.com/robertburbidge/imputing-molecular-features).\n\nMy goal was to understand a little more about molecular structures and to improve my knowledge of machine learning with python. I certainly accomplished it, and will try to see if my feature ideas make sense and are capable of producing better prediction results."},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}