{"cells":[{"metadata":{},"cell_type":"markdown","source":"\n\n# Reach Top-54 using only LGBM\n\n\n\n\n\n
\n\nWe had problems with the GPUs since Kaggle only lets you use one. MPNN was running on GCP and we used kernels for SchNet + NN but the performance wasn't that good, so we tried to improve our **lgbm baseline** ```-1.956 public LB``` with the aim of being able to run 8 lgbm models at the same time.\n\nAt my post [42nd Solution, explanations and apologies](https://www.kaggle.com/c/champs-scalar-coupling/discussion/106263#latest-611323) you can find more information.\nAlso I realized that the **43rd place team** has a similar approach, and similar LGBM, check their post [43 place solution of 2 Experts and the farmer](https://www.kaggle.com/c/champs-scalar-coupling/discussion/106280#latest-611320)\n\n**References**\n\n- LGBM based on [Distance - is all you need. LB -1.481](https://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481) by @criskiev \n\nother important kernels:\n- [Using RDKit for Atomic Feature and Visualization](https://www.kaggle.com/sunhwan/using-rdkit-for-atomic-feature-and-visualization) by @sunhwan \n- [Molecule with OpenBabel](https://www.kaggle.com/jmtest/molecule-with-openbabel) by @jmtest \n"},{"metadata":{},"cell_type":"markdown","source":"# Pipeline\n\n0. Data: we decided to split **1JHC** into 2 (you will see *1JHC_0* and *1JHC_1*) due to the discussions.\n1. Features.\n2. Manual finetuning\n3. Bayessian Optimization (1 morning)\n4. Run 8 models (1 per type) using 8 kernels and at the end stack the solutions\n"},{"metadata":{},"cell_type":"markdown","source":"# Cross Validation\n\n| type | CV | \n|---|---|---|\n|1JHC | -0.91 |\n| 1JHN | -1.52 |\n| 2JHH | -2.49 |\n| 2JHC | -1.772 | \n| 2JHN | -2.319 |\n| 3JHH | -2.46 |\n| 3JHC | -1.701 |\n| 3JHN | -2.64 |\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Features\n\n![\"inbox-2779868-c7a27bb40ebfeb48f2f68ddd76e57f97-Screenshot-from-2019-08-25-09-42-58\"](\"https://i.ibb.co/7kdn8s7/inbox-2779868-c7a27bb40ebfeb48f2f68ddd76e57f97-Screenshot-from-2019-08-25-09-42-58.png\")
\n\nI started with **distance features** and then started to add based on my knowledge, papers, discussions etc\n\n### 1. Distance Features. \n\nDistance between C-C bonds is important. My baseline was obviously the kernel [Distance - is all you need. LB -1.481](https://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481) by @criskiev . Distances between atom helps to know more about the geometry and strenghts, bond type, electonegativity... remember the atoms have charge and attract and repel each other.\n\n### 2. Angles: Bond Angles (2J) and Dihedral angels (3J)\n\nWe used this kernel in order to get this features: [Molecule with OpenBabel](https://www.kaggle.com/jmtest/molecule-with-openbabel) by @jmtest \n\n**Karplus Ecuation**\n\nBothner-By equation\nJHH = 7 -cos Θ + 5 cos 2Θ\n\nwhere Θ is the torsion angle ... the problem wasn't obtain those angles, the problem was that I checked openbabel and RDKit and were not well calculated! For example we have **CH4** (1st molecule), the angles should be 60 but I obtained random numbers like: 83, 74,57 etc ([check here](http://www.ochempal.org/index.php/alphabetical/c-d/dihedral-angle/)).\nEven the bond angles, I tried with water **H20** and insted of 104.5 I obtained random results like: 122, 97...\nI read about it [here](https://www.rdkit.org/docs/GettingStartedInPython.html) and in the case of RDKit, it uses an algorithm based on distance geometry for conorming molecules from 3D (xyz), probably that's the reason :(\n\n> Note that the conformations that result from this procedure tend to be fairly ugly. They should be cleaned up using a force field. This can be done within the RDKit using its implementation of the Universal Force Field (UFF).\n\n\n### 3. Bond type\n\nIf the bond is simple, double, triple ... we didn't count the number of each type, but I could be a molecule feature. As I said, in the theory all these features are connected, based on the distances and atom types you can guess the bond **type**. However, I think is way better to add this features the explicit form when training/feeding the model.\n\n### 4. Atom features\n\n\n- Atom type\n- Hybridization\n- Aromatization\n- Electronegativity\n- Valences \n- Charges\n\nWe used RDKit and I learned how to get them in the kernel - [Using RDKit for Atomic Feature and Visualization](https://www.kaggle.com/sunhwan/using-rdkit-for-atomic-feature-and-visualization) by @sunhwan.\n\n### 5. Molecule Features\n\nHere you could use molecular properties like molecular polarity, potential energy etc. We added number of components and substitutens like numer of benzenes, number of aromatic nitrogens, number of hydroxil groups etc.\nWe used RDKit, see: [rdkit.Chem.Fragments](http://rdkit.org/docs_temp/source/rdkit.Chem.Fragments.html) (Number of aromatic nitrogens, Number of carboxylic acids etc)\n\nImportant post: **[Is 1JHC really one class](https://www.kaggle.com/c/champs-scalar-coupling/discussion/104241#latest-606224)**\n\n> Giba: The two groups are easily splited setting a threshold in 1J coupling distance to 1.065.\n\nThe reason is this:\n\n![\"inbox-2779868-4cef58f29e1062c9271e1b17e9bc80f9-Screenshot-from-2019-08-25-09-43-22\"](\"https://i.ibb.co/2j9DG4r/inbox-2779868-4cef58f29e1062c9271e1b17e9bc80f9-Screenshot-from-2019-08-25-09-43-22.png\")
\n\n**References**\n\n- https://www.ucl.ac.uk/nmr/NMR_lecture_notes/L3_3_97_web.pdf\n- https://www.chem.wisc.edu/areas/reich/nmr/notes-5-hmr-5-vic-coupling.pdf\n- https://www.chem.wisc.edu/areas/reich/nmr/05-hmr-05-3j.htm\n\n"},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T09:58:27.940391Z","start_time":"2019-08-19T09:58:26.307970Z"},"scrolled":false,"id":"bNC7bQBVZOfu","colab_type":"code","colab":{},"_kg_hide-input":false,"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nfrom sklearn.metrics import mean_absolute_error\nimport pandas as pd\nimport numpy as np \nfrom sklearn.model_selection import KFold, StratifiedKFold\nimport lightgbm as lgb\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Basic config"},{"metadata":{"trusted":true},"cell_type":"code","source":"PATH = '../input/mol-features'","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T14:28:57.364457Z","start_time":"2019-08-19T14:28:57.360950Z"},"id":"n_D9Tnf0ZOf5","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"type_name = {\n 0: '1_1JHC_0',\n 1: '1_1JHC_1',\n 2: '2_2JHH',\n 3: '3_1JHN',\n 4: '4_2JHN',\n 5: '5_2JHC',\n 6: '6_3JHH',\n 7: '7_3JHC',\n 8: '8_3JHN', \n }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/mol-features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/mol-features/1_1jhc-20190826t180727z-001 #compressed file","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"folder = {\n '1_1JHC_0': '1_1jhc_0-20190825t153133z-001',\n '1_1JHC_1': '1_1jhc_1-20190826t180741z-001',\n '2_2JHH': '2_2jhh-20190825t170952z-001',\n '3_1JHN': '3_1jhn-20190825t171137z-001',\n '4_2JHN': '4_2jhn-20190825t170517z-001',\n '5_2JHC':'5_2jhc-20190825t175318z-001',\n '6_3JHH':'6_3jhh-20190825t170153z-001',\n '7_3JHC':'7_3jhc-20190825t153045z-001',\n '8_3JHN':'8_3jhn-20190825t170413z-001'\n \n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Parameters\n\nWe did manual finetuning + bayesian optimization."},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T14:29:01.107228Z","start_time":"2019-08-19T14:29:01.101940Z"},"id":"4v0Y9c0IZOgD","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"type_params = {\n \n 0: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'num_leaves': 511,\n 'sub_feature': 0.50,\n 'sub_row': 0.5,\n 'bagging_freq': 1,\n 'metric': 'mae'},\n\n 1: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'num_leaves': 100,\n 'sub_feature': 0.50,\n 'sub_row': 0.5,\n 'bagging_freq': 1,\n 'metric': 'mae'},\n \n 2: { \n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.01,\n 'bagging_freq': 1,\n 'metric': 'mae',\n 'min_data_in_leaf': 130, \n 'num_leaves': 150, \n 'reg_alpha': 0.5, \n 'reg_lambda': 0.6000000000000001, \n 'sub_feature': 0.30000000000000004, \n 'sub_row': 0.4},\n \n \n 3: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.01,\n 'bagging_freq': 1,\n 'metric': 'mae',\n 'min_data_in_leaf': 96, \n 'num_leaves': 30, \n 'reg_alpha': 0.2, \n 'reg_lambda': 0.4, \n 'sub_feature': 0.4, \n 'sub_row': 0.5},\n \n \n 4: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'bagging_freq': 1,\n 'metric': 'mae',\n 'min_data_in_leaf': 21, \n 'num_leaves': 200, \n 'reg_alpha': 0.30000000000000004, \n 'reg_lambda': 0.2, \n 'sub_feature': 0.4, \n 'sub_row': 1.0},\n \n 5: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'bagging_freq': 1,\n 'metric': 'mae',\n 'num_leaves': 1023, \n 'sub_feature': 0.5, \n 'sub_row': 0.5},\n \n 6:{\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'min_data_in_leaf': 50, \n 'num_leaves': 700, \n 'reg_alpha': 0.30000000000000004, \n 'reg_lambda': 0.8, \n 'sub_feature': 0.5, \n 'sub_row': 0.5},\n \n 7: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.1,\n 'bagging_freq': 1,\n 'metric': 'mae',\n 'num_leaves': 1023, \n 'sub_feature': 0.5, \n 'sub_row': 0.5},\n \n\n 8: {\n 'boosting_type': \"gbdt\",\n 'objective': \"huber\",\n 'learning_rate': 0.01,\n 'min_data_in_leaf': 38,\n 'num_leaves': 350,\n 'reg_alpha': 0.30000000000000004,\n 'reg_lambda': 0.6000000000000001,\n 'sub_feature': 0.6000000000000001,\n 'sub_row': 0.5,\n 'metric': 'mae'}\n\n\n}","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T14:29:06.557881Z","start_time":"2019-08-19T14:29:05.586101Z"},"id":"ky1-RLbaZOgK","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(f'../input/champs-scalar-coupling/sample_submission.csv', low_memory=False)\nsub ['typei'] = pd.read_csv(f'../input/champs-scalar-coupling/test.csv', low_memory=False)['type']\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Run\n\n**Important** we used 8 kaggle kernels, that means run 8 different models at the same time. In this part you should indicate a specific *bond type* by modifying ```type_name``` (a dictionary at the beginning)"},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T16:55:45.804343Z","start_time":"2019-08-19T14:29:10.522665Z"},"scrolled":true,"id":"NoNVlf4mZOgQ","colab_type":"code","outputId":"da600548-cd69-4352-ba69-f3ffa4b7a69a","colab":{},"trusted":true},"cell_type":"code","source":"\nscore = []\n\n# fit 8(type) model\nfor idx in type_name:\n\n ntype = type_name[idx]\n print(ntype)\n direct = folder[ntype]\n x_train_val = pd.read_csv(PATH+f'/{direct}/{ntype}/x_train.csv', index_col=0, low_memory=False)\n print('x_train:', x_train_val.shape)\n x_test_val =pd.read_csv(PATH+f'/{direct}/{ntype}/x_test.csv', index_col=0, low_memory=False)\n print('x_test:', x_test_val.shape)\n ID = pd.read_pickle(PATH+f'/{direct}/{ntype}/ID.csv')\n y_train_val = pd.read_pickle(PATH+f'/{direct}/{ntype}/y_train.csv')\n\n break\n \n print(f'------------type{ntype}------------')\n\n maes = []\n predictions = np.zeros(len(x_test_val))\n preds_train = np.zeros(len(x_train_val))\n\n n_fold = 5\n folds = StratifiedKFold(n_splits=n_fold, shuffle=False, random_state=42)\n\n for fold_, (trn_idx, val_idx) in enumerate(folds.split(x_train_val, ID)):\n strLog = \"fold {}\".format(fold_)\n print(strLog)\n\n x_tr, x_val = x_train_val.iloc[trn_idx], x_train_val.iloc[val_idx]\n y_tr, y_val = y_train_val.iloc[trn_idx], y_train_val.iloc[val_idx]\n\n model = lgb.LGBMRegressor(**type_params[idx], n_estimators=30000,random_state=1)\n model.fit(x_tr,\n y_tr,\n eval_set=[(x_tr, y_tr), (x_val, y_val)],\n eval_metric='mae',\n verbose=1000,\n early_stopping_rounds=200\n )\n\n # predictions\n preds = model.predict(\n x_test_val) # , num_iteration=model.best_iteration_)\n predictions += preds / folds.n_splits\n preds = model.predict(\n x_train_val) # , num_iteration=model.best_iteration_)\n preds_train += preds / folds.n_splits\n\n preds = model.predict(x_val)\n\n # mean absolute error\n mae = mean_absolute_error(y_val, preds)\n print('MAE: %.6f' % mae)\n print('Score: %.6f' % np.log(mae))\n maes.append(mae)\n print('')\n\n sub.loc[sub['typei'] == idx, 'scalar_coupling_constant'] = predictions\n score.append(np.mean(maes))\n print(f'{ntype} MAE:', np.mean(maes))\n print(f'{ntype} Score:', np.log(np.mean(maes)))\n\n print('')\n\nprint('')\nprint('----------------------')\n# print('train score:', sum(np.log(score)) / 8)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T16:59:45.638663Z","start_time":"2019-08-19T16:59:36.384090Z"},"id":"GOhmJGBcZOgY","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"#sub.to_csv(f'{config.DATA_DIR}/ano_20000.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"ExecuteTime":{"end_time":"2019-08-19T16:55:46.454876Z","start_time":"2019-08-19T16:55:46.430700Z"},"id":"NhrI7r13ZOge","colab_type":"code","colab":{},"trusted":true},"cell_type":"code","source":"final = sub.drop(['typei'], axis=1)\nfinal.to_csv(f'example.csv', index=False)\nfinal.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now I'm going to use the submission we obtained"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/molsubs/lgbm_final.csv')\nsubmission.to_csv('lgbm_final.csv', index=False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion"},{"metadata":{},"cell_type":"markdown","source":"At the post: [1JHC best CV score?](https://www.kaggle.com/c/champs-scalar-coupling/discussion/98444#latest-602104) CPMP posted this:\n\n\n\n> I don't think lgb is competitive here, will move to graph NNs now.\n\nWe used that as reference and definitely he was right (not a surprise), but LGBM performs really good with the right **features**, and actually our LGBM single model (well, 1 model per type) could reach top-54.\n\n> **NOTE:** we didn't submit LGBM alone (during the last week), so we realized about this yesterday :)\n"}],"metadata":{"hide_input":false,"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"},"toc":{"base_numbering":1,"nav_menu":{},"number_sections":true,"sideBar":true,"skip_h1_title":false,"title_cell":"Table of Contents","title_sidebar":"Contents","toc_cell":false,"toc_position":{},"toc_section_display":true,"toc_window_display":false},"varInspector":{"cols":{"lenName":16,"lenType":16,"lenVar":40},"kernels_config":{"python":{"delete_cmd_postfix":"","delete_cmd_prefix":"del ","library":"var_list.py","varRefreshCmd":"print(var_dic_list())"},"r":{"delete_cmd_postfix":") ","delete_cmd_prefix":"rm(","library":"var_list.r","varRefreshCmd":"cat(var_dic_list()) "}},"types_to_exclude":["module","function","builtin_function_or_method","instance","_Feature"],"window_display":false},"colab":{"name":"lgbm_model.ipynb","version":"0.3.2","provenance":[]}},"nbformat":4,"nbformat_minor":1}