{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties\n## Can you measure the magnetic interactions between a pair of atoms?\n__CHAMPS (CHemistry And Mathematics in Phase Space)__"},{"metadata":{},"cell_type":"markdown","source":"`url:   `[https://www.kaggle.com/c/champs-scalar-coupling/overview](https://www.kaggle.com/c/champs-scalar-coupling/overview)    \n`Date:     1/Jul/2019`    \n`Updated:  2/Jul/2019`   \n`Author:   Enrique Pérez Herrero`"},{"metadata":{},"cell_type":"markdown","source":"# Predict Mulliken Charges with catboost\nThis netobook predict missing Mulliken charges values for the test data. The input features are available in the competition csv files. The model used is a Catboost regressor. The output file can be used to update structures data."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nfrom matplotlib import rcParams\nimport seaborn as sns\n\nfrom catboost import Pool, CatBoostRegressor, cv\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error as MAE\nfrom sklearn.metrics import mean_squared_error, r2_score","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# User Functions"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# User functions\nDATA_PATH = '../input'\nAUX_PATH = '../input/easy-potential-energy-prediction'\n\n\ndef csv_path(dataset=\"train\", data_path=DATA_PATH):\n    \"\"\"\n    Get csv data path\n    \"\"\"\n    return \"{}/{}.csv\".format(data_path, dataset)\n\n\ndef read_file(data_file, data_path=DATA_PATH):\n    \"\"\"\n    Read csv files from data path\n    \"\"\"\n    assert \"csv\" in data_file\n    file_ext = data_file.split(\".\")[1]\n    dataset = data_file.split(\".\")[0]\n    index_col = None\n    index_type = ['train', 'test']\n    if dataset in index_type:\n        index_col = 'id'\n    data_path = csv_path(dataset, data_path=data_path)\n    return pd.read_csv(data_path, index_col=index_col)\n\n\n# Atom properties\nATOM_ELECTRONEG = {\n    \"H\": 2.2,\n    \"C\": 2.55,\n    \"N\": 3.04,\n    \"O\": 3.44,\n    \"F\": 3.98\n}\n\nATOM_NUMBER = {\n    \"H\": 1,\n    \"C\": 6,\n    \"N\": 7,\n    \"O\": 8,\n    \"F\": 9\n}\n\n# 1 atomic mass unit (amu) corresponds to 1.660539040 × 10−24 gram\nATOM_MASS = {\n    \"H\": 1.00784,\n    \"C\": 12.0107,\n    \"N\": 14.0067,\n    \"O\": 15.9990,\n    \"F\": 18.9984\n}\n\n# Without fudge factor\nATOM_RADIUS = {\n    'H': 0.38,\n    'C': 0.77,\n    'N': 0.75,\n    'O': 0.73,\n    'F': 0.71\n}\n\n\ndef d_0(struc_df):\n    \"\"\"\n    Atom distance to origin\n    struc_df: structure data frame\n    \"\"\"\n    return np.sqrt((struc_df.x ** 2) +\n                   (struc_df.y ** 2) +\n                   (struc_df.z ** 2))\n\n\ndef d_max(struc_df):\n    \"\"\"\n    Atom max coord to origin\n    df: structure data frame\n    \"\"\"\n    return np.max([struc_df.x, struc_df.y, struc_df.z], axis=0)\n\n\ndef d_min(struc_df):\n    \"\"\"\n    Atom min coord to origin\n    df: structure data frame\n    \"\"\"\n    return np.min([struc_df.x, struc_df.y, struc_df.z], axis=0)\n\n\ndef cos_x(struc_df):\n    \"Atom cos axis x\"\n    return struc_df.x / d_0(struc_df)\n\n\ndef cos_y(struc_df):\n    \"Atom cos axis y\"\n    return struc_df.y / d_0(struc_df)\n\n\ndef cos_z(struc_df):\n    \"Atom cos axis z\"\n    return struc_df.z / d_0(struc_df)\n\n\ndef structures_add_features(structures):\n    \"\"\"\n    Add new features to \"structures.csv\"\n    \"\"\"\n    atoms_values = structures['atom'].values\n    structures[\"atom_mass\"] = [ATOM_MASS[x] for x in atoms_values]\n    structures[\"atom_number\"] = [ATOM_NUMBER[x] for x in atoms_values]\n    structures[\"atom_electroneg\"] = [ATOM_ELECTRONEG[x] for x in atoms_values]\n    structures[\"atom_radius\"] = [ATOM_RADIUS[x] for x in atoms_values]\n    structures[\"atom_d0\"] = d_0(structures)\n    structures[\"atom_dmax\"] = d_max(structures)\n    structures[\"atom_dmin\"] = d_min(structures)\n    structures[\"cos_x\"] = cos_x(structures)\n    structures[\"cos_y\"] = cos_y(structures)\n    structures[\"cos_z\"] = cos_z(structures)\n    return structures\n\n\ndef map_atom_info(data_df, structures_df):\n    \"\"\"\n    Merge atoms info in structures\n    \"\"\"\n    index_name = data_df.index.name\n    data_df.reset_index(inplace=True)\n    for atom_idx in range(2):\n        former_cols = list(data_df)\n        data_df = pd.merge(data_df,\n                           structures_df,\n                           how='left',\n                           left_on=['molecule_name', f'atom_index_{atom_idx}'],\n                           right_on=['molecule_name', 'atom_index'],\n                           left_index=False,\n                           right_index=False)\n        data_df = data_df.drop('atom_index', axis=1)\n        actual_cols = list(data_df)\n        new_cols = [x + \"_{}\".format(atom_idx) for x in actual_cols if x not in former_cols]\n        new_cols = former_cols + new_cols\n        data_df.columns = new_cols\n    data_df.set_index(index_name, inplace=True)\n    return data_df\n\n\ndef dist_xyz_2(data_df):\n    \"\"\"\n    Coupling Euclidean distance\n    \"\"\"\n    return ((data_df.x_1 - data_df.x_0) ** 2 +\n            (data_df.y_1 - data_df.y_0) ** 2 +\n            (data_df.z_1 - data_df.z_0) ** 2)\n\n\ndef dist_xyz(data_df):\n    \"\"\"\n    Coupling Euclidean distance\n    \"\"\"\n    return np.sqrt(dist_xyz_2(data_df))\n\n\ndef dist_x(data_df):\n    \"\"\"\n    Coupling X axis distance\n    \"\"\"\n    return np.sqrt((data_df.x_1 - data_df.x_0) ** 2)\n\n\ndef dist_y(data_df):\n    \"\"\"\n    Coupling Y axis distance\n    \"\"\"\n    return np.sqrt((data_df.y_1 - data_df.y_0) ** 2)\n\n\ndef dist_z(data_df):\n    \"\"\"\n    Coupling Z axis distance\n    \"\"\"\n    return np.sqrt((data_df.z_1 - data_df.z_0) ** 2)\n\n\ndef data_add_features(data_df, structures_df):\n    \"\"\"\n    Add new features to train or test data frames\n    \"\"\"\n    data_df = map_atom_info(data_df, structures_df)\n    # Add distance features\n    data_df[\"dist_xyz_2\"] = dist_xyz_2(data_df)\n    data_df[\"dist_xyz\"] = dist_xyz(data_df)\n    data_df[\"dist_x\"] = dist_x(data_df)\n    data_df[\"dist_y\"] = dist_y(data_df)\n    data_df[\"dist_z\"] = dist_z(data_df)\n    # atom_0 is always H\n    data_df.drop(\"atom_0\", axis=1, inplace=True)\n    return data_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Parameters"},{"metadata":{"_kg_hide-input":false,"_kg_hide-output":false,"trusted":true},"cell_type":"code","source":"# kFold Validation Parameters\nRANDOM_STATE = 123\nN_SPLITS = 3\nSHUFFLE = True\nVERBOSE = False\n\n# Script parameters\nFREE_MEMORY = True\nSAVE_SUBMISSION = True\nMODEL_NAME = \"CB_MULLIKEN_001\"\nLINE = 40 * \"-\"\n\n# Ploting Parameters\nFIGSIZE = (10, 6)\nsns.set()\n\n# Data files\nINPUT_FILE_A = \"mulliken_charges.csv\"\nINPUT_FILE_B = \"structures.csv\"\nOUTPUT_FILE = \"mulliken_charges_upd.csv\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read data files"},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_charges = read_file(INPUT_FILE_A)\nstructures = read_file(INPUT_FILE_B)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_charges.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Add new features to structures data frame"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = structures_add_features(structures)\nstructures.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train and test Mulliken data frames"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Basic features\nid_feature = 'molecule_name'\ntarget_feature = 'mulliken_charge'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_train = mulliken_charges.merge(structures)\n# Check merged data frame\nassert mulliken_train.shape[0] == mulliken_charges.shape[0]\nassert mulliken_train.notnull().values.any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_test = mulliken_charges.merge(structures,\n                                       on=[id_feature, 'atom_index'],\n                                       how='right',\n                                       indicator=False)\n\nmulliken_test = mulliken_test[mulliken_test[target_feature].isnull()]\nmulliken_test = mulliken_test.drop(target_feature, axis=1)\n# Check merged data frame\nstructures_shape = structures.shape\nassert mulliken_train.shape[0] + mulliken_test.shape[0] == structures_shape[0]\nassert mulliken_test.notnull().values.any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_train.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"selected_features = list(mulliken_test)\nselected_features.remove(id_feature)\ncategorical_features = ['atom']\nprint(\"Selected Features: \\t{}\".format(selected_features))\nprint(\"Target Feature: \\t{}\".format(target_feature))\nprint(\"Id Feature: \\t\\t{}\".format(id_feature))\nprint(\"Categorical Features: \\t{}\".format(categorical_features))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if FREE_MEMORY:\n    del structures, mulliken_charges","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## k-Fold Validation"},{"metadata":{"trusted":true},"cell_type":"code","source":"X = mulliken_train[selected_features]\ny = mulliken_train[target_feature]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kfold = KFold(n_splits=N_SPLITS,\n              random_state=RANDOM_STATE,\n              shuffle=SHUFFLE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Catboost Regressor Parameters"},{"metadata":{},"cell_type":"markdown","source":"`use_best_model` parameter must be set to `True` during validation,  but in the final fitting, when there is no evaluation pool it must be set to `False`. The maximum `depth` allowed is 16"},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {\n    \"model_name\": MODEL_NAME,\n    \"iterations\": 100,\n    \"depth\": 16,\n    \"learning_rate\": 0.80,\n    \"reg_lambda\": 3.0,\n    \"loss_function\": \"RMSE\",\n    \"verbose\": VERBOSE,\n    \"random_seed\": RANDOM_STATE,\n    \"task_type\": \"GPU\"\n}\n\ncat_reg = CatBoostRegressor(iterations=params[\"iterations\"],\n                            depth=params[\"depth\"],\n                            learning_rate=params[\"learning_rate\"],\n                            loss_function=params[\"loss_function\"],\n                            verbose=params[\"verbose\"],\n                            use_best_model=True,\n                            reg_lambda=params[\"reg_lambda\"],\n                            random_seed=params[\"random_seed\"],\n                            task_type=params[\"task_type\"],\n                            name=params[\"model_name\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## k-Fold Validation Training"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfold = 0\nr2_scores = []\nmse_scores = []\nmae_scores = []\n\n\nfor in_index, oof_index in kfold.split(X, y):\n    fold += 1\n    print(LINE)\n    print(\"- Training Fold: ({}/{})\".format(fold, N_SPLITS))\n    X_in, X_oof = X.loc[in_index], X.loc[oof_index]\n    y_in, y_oof = y.loc[in_index], y.loc[oof_index]\n    # Train and evaluation data pools\n    train_pool = Pool(data=X_in,\n                      label=y_in,\n                      cat_features=categorical_features)\n    eval_pool = Pool(data=X_oof,\n                     label=y_oof,\n                     cat_features=categorical_features)\n    # Fit Regressor\n    hist = cat_reg.fit(train_pool, eval_set=eval_pool)\n    y_pred = cat_reg.predict(X_oof)\n    # Metrics\n    r2 = r2_score(y_oof, y_pred)\n    r2_scores.append(r2)\n    mse_score = mean_squared_error(y_oof, y_pred)\n    mse_scores.append(mse_score)\n    mae_score = MAE(y_oof, y_pred)\n    mae_scores.append(mae_score)\n    print(f'\\t R2:  {r2:.4f}')\n    print(f'\\t MSE: {mse_score:.4f}')\n    print(f'\\t MAE: {mse_score:.4f}')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## k-Fold metrics\nprint('kFold Validation Results:')\nprint(' * Average Variance Score (R2): \\t{:.4f}'.format(np.mean(r2_scores)))\nprint(' * Average Mean squared error (MSE): \\t{:.4f}'.format(np.mean(mse_score)))\nprint(' * Average Mean absolure error (MAE): \\t{:.4f}'.format(np.mean(mae_score)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Plot last fold prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Increase the number of params that matplotlib can handle\nrcParams['agg.path.chunksize'] = 10000\n# Perfect prediction line\nperfect_pred = np.arange(-0.8, 0.9, 0.1)\nplt.figure(figsize=FIGSIZE)\nplt.plot(perfect_pred , perfect_pred , c=\"k\")\nplt.scatter(y_oof, y_pred, s=0.2, alpha=0.1)\nplt.title(\"Fold {} Mulliken Charges Prediction\".format(fold))\nplt.xlabel(\"Validation\")\nplt.ylabel(\"Predicted\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature importances"},{"metadata":{"trusted":true},"cell_type":"code","source":"importances_df = pd.DataFrame({\"features\":list(X),\n                               \"importances\":cat_reg.feature_importances_})\nimportances_df = importances_df.sort_values('importances', ascending=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"importances_df.plot.barh(x='features',\n                         y='importances',\n                         figsize=FIGSIZE,\n                         legend=False,\n                         width=0.9)\nplt.title(\"Catboost Regression Feature Importances\")\nplt.ylabel(\"\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train Full Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_pool = Pool(data=X,\n                  label=y,\n                  cat_features=categorical_features)\n\nfull_cat_reg = CatBoostRegressor(iterations=params[\"iterations\"],\n                                 depth=params[\"depth\"],\n                                 learning_rate=params[\"learning_rate\"],\n                                 loss_function=params[\"loss_function\"],\n                                 verbose=params[\"verbose\"],\n                                 use_best_model=False,\n                                 reg_lambda=params[\"reg_lambda\"],\n                                 random_seed=params[\"random_seed\"],\n                                 task_type=params[\"task_type\"],\n                                 name=params[\"model_name\"])\n# Fit Regressor\nhist = full_cat_reg.fit(train_pool)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Predict test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_test = full_cat_reg.predict(mulliken_test[selected_features])\nmulliken_test[target_feature] = y_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot train and test Mulliken charges density"},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken_train[target_feature].plot.kde(figsize=FIGSIZE, legend=True, label=\"train\")\nmulliken_test[target_feature].plot.kde(figsize=FIGSIZE, legend=True, label=\"test\")\nplt.title(\"Train and test Mulliken Charges Density\")\nplt.ylabel(\"\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Merge train and test in one single data frame"},{"metadata":{"trusted":true},"cell_type":"code","source":"output_features = [id_feature, \"atom_index\", target_feature]\n\nmulliken_charges_upd = pd.concat([\n    mulliken_train[output_features],\n    mulliken_test[output_features]],\n    ignore_index=True)\n    \nif FREE_MEMORY:\n    del mulliken_train, mulliken_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save updated mulliken charges data"},{"metadata":{"trusted":true},"cell_type":"code","source":"assert mulliken_charges_upd.shape[0] == structures_shape[0]\nassert mulliken_charges_upd.notnull().values.any()\nprint(f\"Saving file {OUTPUT_FILE}...\")\nmulliken_charges_upd.to_csv(OUTPUT_FILE, index=False)","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}