{"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":"\n`url:   `[https://www.kaggle.com/c/champs-scalar-coupling/overview](https://www.kaggle.com/c/champs-scalar-coupling/overview)    \n`Date:     23/Jun/2019`    \n`Updated:  23/Jun/2019`   \n`Author:   Enrique Pérez Herrero`"},{"metadata":{},"cell_type":"markdown","source":"# Predict Potential Energy\n\nThis notebook uses the number of dictinct atoms in a molecule to predict, accurately, the molecular potential energy.    The coefficients and the intercept of the linear regression give a closed formula to derive this feature."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import warnings\nimport numpy as np\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error, r2_score\nfrom sklearn.linear_model import LinearRegression","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Ploting Parameters\nFIGSIZE = (10, 6)\nsns.set() # Set seaborn style\n\n# Notebook Parameters\nFREE_MEMORY = True\nSAVE_CSV = True\n\n# kFold Validation Parameters\nRANDOM_STATE = 123\nN_SPLITS = 3\nSHUFFLE = True\nVERBOSE = False\n\n# Data Files\nOUTPUT_FILE_A = \"potential_energy_upd.csv\"\nOUTPUT_FILE_B = \"molecule_train.csv\"\nOUTPUT_FILE_C = \"molecule_test.csv\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## User functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Input data user functions\n\nDATA_PATH = \"../input\"\n\ndef csv_path(dataset=\"train\", data_path=DATA_PATH):\n    return \"{}/{}.csv\".format(data_path, dataset)\n\ndef read_data(dataset='train', data_path=DATA_PATH):\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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Read data and arrange Data Frames"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = read_data(\"train\")\ntest = read_data(\"test\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## List all atoms in the same molecule"},{"metadata":{"trusted":true},"cell_type":"code","source":"# One row per molecule\nmolecule_train = pd.DataFrame({\"molecule_name\" : train[\"molecule_name\"].unique()})\nmolecule_test = pd.DataFrame({\"molecule_name\" : test[\"molecule_name\"].unique()})\nstructures = read_data('structures')\natom_list_df = structures.groupby('molecule_name')['atom'].apply(list)\natom_list_df = atom_list_df.to_frame()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if FREE_MEMORY:\n    del train, test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge atoms list to molecule train and test"},{"metadata":{"trusted":true},"cell_type":"code","source":"molecule_train = pd.merge(molecule_train, atom_list_df, how='left', on='molecule_name')\nmolecule_test = pd.merge(molecule_test, atom_list_df, how='left', on='molecule_name')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Count distinct atoms in a molecule\nThe number of distinct atoms in a molecule are the features to use for `Potential Energy` prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Count atoms by type\natoms_list = structures.atom.unique().tolist()\nprint(\"Distinct atoms in structures data: \\n {}\".format(atoms_list))\n\nfor atom in atoms_list:\n    molecule_train['atom_' + atom] = \\\n        molecule_train['atom'].apply(lambda x: x.count(atom))\n    molecule_test['atom_' + atom] = \\\n        molecule_test['atom'].apply(lambda x: x.count(atom))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge Potential Energy"},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy = read_data(\"potential_energy\")\nmolecule_train = pd.merge(molecule_train, potential_energy)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if FREE_MEMORY:\n    del potential_energy, structures\n    del atom_list_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The column `atom` can be useful to calculate the molecular weight. This feature is not useful to predict `potential_energy` but it can be interesting for other predictions in the competition."},{"metadata":{"trusted":true},"cell_type":"code","source":"# 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\ndef mol_weight(atom_list):\n    \"\"\"\n    Get the molecular weight from a list of atoms\n    \"\"\"\n    return sum(map(lambda x: ATOM_MASS[x], atom_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"molecule_train[\"mol_weight\"] = molecule_train.atom.apply(lambda x: mol_weight(x))\nmolecule_test[\"mol_weight\"] = molecule_train.atom.apply(lambda x: mol_weight(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"molecule_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fit Linear Model and Validate Results"},{"metadata":{},"cell_type":"markdown","source":"## Select Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"id_feature = 'molecule_name'\ntarget_feature = (set(molecule_train) - set(molecule_test)).pop()\nselected_features = list(molecule_test)\nselected_features.remove(id_feature)\nselected_features.remove('atom')\nprint(\"Selected Features: \\t{}\".format(selected_features))\nprint(\"Target Feature: \\t{}\".format(target_feature))\nprint(\"Id Feature: \\t\\t{}\".format(id_feature))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train and test data frames"},{"metadata":{"trusted":true},"cell_type":"code","source":"X = molecule_train[selected_features]\ny = molecule_train[target_feature]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Set k-Fold Validation "},{"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":{"trusted":true},"cell_type":"code","source":"fold = 0\nr2_scores = []\nmse_scores = []\nlin_reg = LinearRegression()\n\nfor in_index, oof_index in kfold.split(X, y):\n    fold += 1\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    \n    lin_reg.fit(X_in, y_in)\n    y_pred = lin_reg.predict(X_oof)\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\nprint('\\nkFold 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)))\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Explained variance score R2 equal to 1 is perfect prediction"},{"metadata":{},"cell_type":"markdown","source":"## Plot last fold validation and predicted Potential Energy"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=FIGSIZE)\nplt.plot(y_oof, y_pred)\nplt.title(\"Fold {} Prediction\".format(fold))\nplt.xlabel(\"Validation Potential Energy\")\nplt.ylabel(\"Predicted Potential Energy\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict Potential Energy for test data\nPrediction using all available train data "},{"metadata":{"trusted":true},"cell_type":"code","source":"lin_reg.fit(X, y)\ny_test = lin_reg.predict(molecule_test[selected_features])\nmolecule_test[target_feature] = y_test\nif FREE_MEMORY:\n    del X, y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Concatenate `molecule_train` and `molecule_test` in the same data frame"},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy_upd = pd.concat([molecule_train[[id_feature, target_feature]],\n                                  molecule_test[[id_feature, target_feature]]],\n                                 ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=FIGSIZE)\nmolecule_train['potential_energy'].plot.kde(figsize=FIGSIZE, legend=True, label=\"train\")\nmolecule_test['potential_energy'].plot.kde(figsize=FIGSIZE, legend=True, label=\"test\")\nplt.title(\"Predicted and Training Potential Energy Density Plot\")\nplt.xlabel(\"Potential Energy\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Arrange updated potencial energy data frame"},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy_upd = potential_energy_upd.sort_values(id_feature)\npotential_energy_upd.reset_index(drop=True, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energy_upd.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Save data as csv\nThis file contains training and predicted testing potential energy, and it can be useful to add this feature in the `scalar_coupling_constant` competition.    \n__Files__\n* \"potential_energy_upd.csv\": Updated `potential_energy.csv` with predicted values for test data.\n* \"molecule_train.csv\": Molecule features, including potential energy, with molecules found in train data.\n* \"molecule_test.csv\": Molecule features, including potential energy, with molecules found in train data."},{"metadata":{"trusted":true},"cell_type":"code","source":"def save_csv(df, file_name):\n    \"\"\"\n    Check and save csv data to kernel main folder \n    \"\"\"\n    assert \"csv\" in file_name, \"Bad file extension\"\n    assert df.notnull().values.any(), \"NaN values in data frame\"\n    print(f\"Saving file {file_name}...\")\n    df.to_csv(file_name, index=False)\n\nif SAVE_CSV:\n    save_csv(potential_energy_upd, OUTPUT_FILE_A)\n    save_csv(molecule_train, OUTPUT_FILE_B)\n    save_csv(molecule_test, OUTPUT_FILE_C)","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}