{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f7b6d7b6-698e-4ea2-b8e9-bad6f279fcfe",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jwchen25/Modeling_Lab_Class/blob/main/Simple_Tutorial_ML4Chem.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd7ee09d-524b-4f90-ac5c-cc37187ba4c4",
   "metadata": {},
   "source": [
    "# Tutorial - Machine Learning in Chemistry\n",
    "\n",
    "## Table of content\n",
    "\n",
    "0. Relevant packages\n",
    "1. Basic handling of molecular data\n",
    "2. Supervised learning\n",
    "3. A simple example of regression task in chemistry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a15d3c35",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 0. Relevant packages\n",
    "\n",
    "### RDKit\n",
    "`RDKit` is an open-source software toolkit for cheminformatics, designed to assist in the analysis and design of small molecules and chemical compounds. It provides a set of libraries and tools for the manipulation and analysis of molecular structures, molecular descriptors, molecular fingerprints, molecular similarity, molecular visualization, and more. The toolkit is widely used in academia, as well as in the pharmaceutical, biotech, and chemical industries for a variety of tasks such as virtual screening, lead optimization, and chemical database management.\n",
    "\n",
    "### Scikit-learn\n",
    "Scikit-learn is an open source machine learning library that supports supervised and unsupervised learning. It also provides various tools for model fitting, data preprocessing, model selection, model evaluation, and many other utilities. We will learn to use scikit-learn to do machine learning work. You can also browse the scikit-learn [user guide](https://scikit-learn.org/stable/user_guide.html) and [tutorials](https://scikit-learn.org/stable/tutorial/index.html) for additional details.\n",
    "### Essential Libraries and Tools \n",
    "Scikit-learn depends on two other Python packages, NumPy and SciPy. For plotting and interactive development, you should also install matplotlib, IPython, and the Jupyter Notebook.\n",
    "- **NumPy** is one of the fundamental packages for scientific computing in Python. In scikit-learn, the NumPy array is the fundamental data structure. Any data you’re using will have to be converted to a NumPy array.\n",
    "- **SciPy** is a collection of functions for scientific computing in Python. Scikit-learn draws from SciPy’s collection of functions for implementing its algorithms.\n",
    "- **Matplotlib** is the primary scientific plotting library in Python. It provides functions for making publication-quality visualizations such as line charts, histograms, scatter plots, and so on.\n",
    "- **Pandas** Python library for data wrangling and analysis. It can ingest from a great variety of file formats and databases, like SQL, Excel files, and comma-separated values (CSV) files.\n",
    "\n",
    "### XGBoost\n",
    "XGBoost (eXtreme Gradient Boosting) is an optimized distributed gradient boosting library designed to be highly efficient, flexible and portable. It implements machine learning algorithms under the Gradient Boosting framework. XGBoost provides a parallel tree boosting (also known as GBDT, GBM) that solve many data science problems in a fast and accurate way. You can also browse the [XGBoost Documentation](https://xgboost.readthedocs.io/en/stable/) for additional details.\n",
    "\n",
    "### DeepChem\n",
    "DeepChem is a high quality open-source toolchain that democratizes the use of deep-learning in chemistry, biology and materials science. It also provides various tools for dataset loader, splitters, molecular featurization, model construction and hyperparameter tuning. You can also browse the [DeepChem Ducumentation](https://deepchem.readthedocs.io/en/latest/) for additional details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44fbf541",
   "metadata": {},
   "source": [
    "We will first install the required libraries. We also need `RDKit` library to process and analyze molecules, like calculating molecular descriptors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee96d46-2b8a-4043-8289-a5944f5311b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install all libraries\n",
    "! pip install numpy scipy matplotlib scikit-learn pandas rdkit xgboost deepchem mordred pycm\n",
    "\n",
    "# Download all data\n",
    "! mkdir data\n",
    "! wget https://raw.githubusercontent.com/schwallergroup/ai4chem_course/main/notebooks/02%20-%20Supervised%20Learning/data/esol.csv -O data/esol.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e2999e5",
   "metadata": {},
   "source": [
    "# 1. Basic handling of molecular data\n",
    "\n",
    "Let's start looking at a specific molecule to play around with RDKit functionalities, and look at caffeine. \n",
    "\n",
    "The name `caffeine` does not contain any structural information on the molecule. Just by having this name a computer does not know how many (heavy) atoms `caffeine` contains and what bonds they form. We need machine-accessible representations. One of the most commonly used ones is `SMILES`.\n",
    "\n",
    "Firstly, import relevant modules of `RDKit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a3a672",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "from rdkit.Chem.Draw import IPythonConsole"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74a740ba",
   "metadata": {},
   "source": [
    "## 1.1 SMILES\n",
    "\n",
    "`SMILES` (Simplified Molecular Input Line Entry System) is a line notation representation of molecular structures. It is a way of representing chemical compounds as strings of characters, which can be easily processed and analyzed by computer algorithms.\n",
    "\n",
    "Each SMILES string consists of symbols that represent the elements in the molecule, as well as brackets and other characters that describe the bonding between the atoms. For example, the SMILES string for ethanol (C2H5OH) would be `CCO`. In SMILES, each carbon atom is represented by the letter \"C\", each hydrogen atom by the letter \"H\", and each oxygen atom by the letter \"O\". The bonding between the atoms is indicated by the arrangement of the characters in the string.\n",
    "\n",
    "SMILES is widely used in cheminformatics and computational chemistry, as it provides a compact and standardized way of representing molecular structures in a machine-readable form. This makes it possible to compare and analyze large numbers of chemical compounds, as well as to generate predictions about their properties and behavior.\n",
    "\n",
    "\n",
    " Look up the SMILES string of `caffeine` on Wikipedia/PubChem, or use this [link](https://opsin.ch.cam.ac.uk/) to get SMILES from molecule name. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7c487b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "caffeine_smiles = '' # fill in the SMILES "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f2e4cce",
   "metadata": {},
   "source": [
    "## 1.2 - Creating and visualizing molecules\n",
    "\n",
    "Let's start with the most basic rdkit action: creating a `Mol` object (or variable, as you prefer). `Mol` objects represent molecules, and can be created from different molecular representations (SMILES, .sdf files, etc.). We will use the basic `MolFromSmiles` function to create a variable `caffeine` representing our caffeine molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d8bc855",
   "metadata": {},
   "outputs": [],
   "source": [
    "caffeine = Chem.MolFromSmiles(caffeine_smiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e53f3c26",
   "metadata": {},
   "source": [
    "We can display the value of a variable in the notebook by typing the name and then running the cell. In this case, we can visualize the molecule this way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda88f8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "caffeine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbb775d4",
   "metadata": {},
   "source": [
    "Another interesting option is saving the mol object as an image. This way, you can download it or save it in your working directory. We can create an image file using the function `MolToImage` and our mol object as the argument. Check the directory to see the image after running the following code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67f59804",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Draw\n",
    "\n",
    "#Create image file\n",
    "im = Draw.MolToImage(caffeine)\n",
    "\n",
    "#Save image as a PNG file in our current directory\n",
    "im.save('caffeine.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6dcd51e",
   "metadata": {},
   "source": [
    "But what is actually a Mol object? Nothing simpler than a graph representing the molecule! In this graph, vertices represents the atoms and edges the bonds in the molecule. Therefore, we can retrieve the atoms and bonds from the graph and play with them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cb24ce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Note: Hydrogen atoms are ignored!\")\n",
    "\n",
    "# Get total number of atoms\n",
    "n_atoms = caffeine.GetNumAtoms()\n",
    "print(f'Total number of atoms: {n_atoms}')\n",
    "\n",
    "# Get total number of bonds\n",
    "n_bonds = caffeine.GetNumBonds()\n",
    "print(f'Total number of bonds: {n_bonds}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba4951d3",
   "metadata": {},
   "source": [
    "We can also get the properties of atoms and bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e3fdadd",
   "metadata": {},
   "outputs": [],
   "source": [
    "IPythonConsole.drawOptions.addAtomIndices = True  # adding indices for atoms\n",
    "IPythonConsole.drawOptions.addBondIndices = False  # not adding indices for bonds\n",
    "IPythonConsole.molSize = 300, 300  # the image size of the molecule\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cd4c53e",
   "metadata": {},
   "outputs": [],
   "source": [
    "atom_symbols = []  # a list of element symbols for all atoms\n",
    "atom_numbers = []  # a list of the corresponding element numbers of all atoms\n",
    "for atom in caffeine.GetAtoms():\n",
    "    atom_symbols.append(atom.GetSymbol())\n",
    "    atom_numbers.append(atom.GetAtomicNum())\n",
    "\n",
    "print(\"The list of element symbols for all atoms:\")\n",
    "print(atom_symbols)\n",
    "print(\"\\nThe list of the corresponding element numbers of all atoms:\")\n",
    "print(atom_numbers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "037fb358",
   "metadata": {},
   "outputs": [],
   "source": [
    "IPythonConsole.drawOptions.addAtomIndices = False  # not adding indices for atoms\n",
    "IPythonConsole.drawOptions.addBondIndices = True  # adding indices for bonds\n",
    "IPythonConsole.molSize = 300, 300  # the image size of the molecule\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8e71bba",
   "metadata": {},
   "outputs": [],
   "source": [
    "bond_types = []  # a list of bond type for all bonds\n",
    "for bond in caffeine.GetBonds():\n",
    "    bond_types.append(str(bond.GetBondType()))\n",
    "\n",
    "print(\"The list of bond type for all bonds:\")\n",
    "print(bond_types)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dba2fc2",
   "metadata": {},
   "source": [
    "# 2. Supervised learning\n",
    "Two major types of supervised machine learning problems:\n",
    "- **Classification**, where the task task is to predict a class label (e.g. what color, what smell, state of aggregation, etc).\n",
    "- **Regression**, where the task is to predict a real number (e.g. solubility in water, yield, selectivity, etc).\n",
    "\n",
    "## 2.1 Algorithms\n",
    "- k-Nearest Neighbors (k-NN)\n",
    "- Linear Models\n",
    "- Support Vector Machines (SVM)\n",
    "- Decision Trees\n",
    "- Ensembles of Decision Trees\n",
    "  - Random forests\n",
    "  - Gradient boosting machines\n",
    "\n",
    "The package [scikit-learn](https://scikit-learn.org/stable/index.html) lets us create a large number ML models in a very conveinent way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86793e89",
   "metadata": {},
   "outputs": [],
   "source": [
    "# linear regressor\n",
    "from sklearn.linear_model import LinearRegression\n",
    "lin_reg = LinearRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18fd6e5e-2960-4446-8a57-920f0e324129",
   "metadata": {},
   "source": [
    "### Exercise: Similarly to the linear regression model above, create a [k-NN regression model](https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsRegressor.html) using scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26f1caf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "### YOUR CODE\n",
    "\n",
    "# import\n",
    "' your code '\n",
    "\n",
    "# create knn model\n",
    "knn_clf = ' your code '\n",
    "\n",
    "### END"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3441adba",
   "metadata": {},
   "source": [
    "Next time, you can browse the scikit-learn [user guide](https://scikit-learn.org/stable/user_guide.html) to learn about supported algorithms and how to create the model you want."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7bd577f-19c7-4121-b0c7-753f1e500356",
   "metadata": {},
   "source": [
    "## 2.2 Model evaluation and data splitting\n",
    "\n",
    "### Why do we need to split the dataset?\n",
    "\n",
    "We want models to learn from data so we can use them in the future, but we also need to know how good our models perform when they see new examples, and so we reserve some part of our dataset for `testing`: we keep it to evaluate how the model would do the real world, with data it hasn't seen before.\n",
    "\n",
    "In addition, typically we implement multiple models and select the best one, but how to assess which one is the best, without revealing the `test set`? Well, take another subset of the data, and this one we call `validation` set. \n",
    "\n",
    "In the end, we end up splitting our data into `training` (used for training the models), `validation` (for selecting the models), and `test` (for testing the resulting models). If you have more time, you can read [this article](https://towardsdatascience.com/how-to-split-data-into-three-sets-train-validation-and-test-and-why-e50d22d3e54c) for more details.\n",
    "\n",
    "### Evaluation metrics\n",
    "The metrics used to evaluate the ML models are very important. The choice of metrics to use influences how model performance is measured and compared. The main evaluation metrics for regression and classification tasks are illustrated below. If you have more time, you can read this [article](https://blog.knoldus.com/model-evaluation-metrics-for-machine-learning-algorithms/) for more details.\n",
    "\n",
    "<div align=\"center\">\n",
    "<img src=\"https://www.oreilly.com/api/v2/epubs/9781492073048/files/assets/mlbf_0407.png\" width=\"500\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a633f2c",
   "metadata": {},
   "source": [
    "## 2.3 The path to a ML model.\n",
    "0. Define the task\n",
    "1. Prepare data & split data\n",
    "2. Choose the model\n",
    "3. Train the model\n",
    "4. Evaluate the model\n",
    "5. Use the model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "048efaf2",
   "metadata": {},
   "source": [
    "# 3. A simple example of regression task in chemistry\n",
    "\n",
    "**Aqueous solubility is one of the key physical properties of interest to a medicinal or agrochemical chemist**. Solubility affects the uptake/distribution of biologically active compounds in living material and the environment, thus affecting their potential efficacy and marketability. However, solubility determination is a time-consuming experiment, and it is useful to be able to assess solubility in the absence of a physical sample. \n",
    "\n",
    "Our goal here will be to build a ML model that can predict **aqueous solubility** of organic molecules.\n",
    "\n",
    "### We will use the [ESOL dataset](https://pubs.acs.org/doi/10.1021/ci034243x) for this task.\n",
    "This dataset contains structures and water solubility data for 1128 compounds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf98acc5",
   "metadata": {},
   "source": [
    "### Load dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf79eb55",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# load dataset from the CSV file (we have downloaded it in the Section 0)\n",
    "esol_df = pd.read_csv('data/esol.csv')\n",
    "esol_df.sample(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "794e236d",
   "metadata": {},
   "source": [
    "The original dataset contains 2 columns, the `smiles` column, which encodes the molecular structure, and the `log solubility (mol/L)` columns represents aqueous solubility of molecules, which is the our target for this task."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2950f29b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get NumPy arrays from DataFrame for the input and target\n",
    "smiles = esol_df['smiles'].values\n",
    "y = esol_df['log solubility (mol/L)'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "343779ed",
   "metadata": {},
   "source": [
    "### Molecular descriptor calculation\n",
    "We need to convert the SMILES strings of molecules into numerical values that can be used as input to the ML models. Many molecular descriptors can be calculated from SMILES strings using software like `RDKit`, `DeepChem` and [Mordred](https://github.com/mordred-descriptor/mordred).\\\n",
    "For this task we will use [DeepChem featurizers](https://deepchem.readthedocs.io/en/latest/api_reference/featurizers.html) to compute molecular descriptors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5961b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we use molecular descriptors from RDKit, like molecular weight, number of valence electrons, maximum and minimum partial charge, etc.\n",
    "from deepchem.feat import RDKitDescriptors\n",
    "featurizer = RDKitDescriptors()\n",
    "features = featurizer.featurize(smiles)\n",
    "print(f\"Number of generated molecular descriptors: {features.shape[1]}\")\n",
    "\n",
    "# Drop the features containing invalid values\n",
    "import numpy as np\n",
    "features = features[:, ~np.isnan(features).any(axis=0)]\n",
    "print(f\"Number of molecular descriptors without invalid values: {features.shape[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d9c63ba",
   "metadata": {},
   "source": [
    "**Exercise**: You can try to use [MACCSKeysFingerprint](https://deepchem.readthedocs.io/en/latest/api_reference/featurizers.html#maccskeysfingerprint) in [DeepChem featurizers](https://deepchem.readthedocs.io/en/latest/api_reference/featurizers.html) to compute molecular descriptors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09f941ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "### YOUR CODE\n",
    "\n",
    "# import the featurizer\n",
    "' your code '\n",
    "\n",
    "# create featurizer\n",
    "mf_featurizer = ' your code '\n",
    "\n",
    "# compute molecular descriptors\n",
    "mf_features = ' your code '"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97180e44",
   "metadata": {},
   "source": [
    "### Feature Selection\n",
    "\n",
    "Feature selection is a crucial step in machine learning that involves selecting the most relevant features (or variables) from a dataset. This process helps to improve the accuracy and efficiency of a model by reducing the amount of noise and redundancy in the data. Essentially, feature selection allows us to focus on the most important information in the dataset, while ignoring the irrelevant or redundant information.\n",
    "\n",
    "Many classes and functions in the [sklearn.feature_selection](https://scikit-learn.org/stable/modules/feature_selection.html) module can be used for feature selection.\\\n",
    "Here we show a simple feature selection using `VarianceThreshold` in `scikt-learn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "326fe67b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we removed all zero-variance features, i.e. features that have the same value in all samples.\n",
    "from sklearn.feature_selection import VarianceThreshold\n",
    "selector = VarianceThreshold(threshold=0.0)\n",
    "features = selector.fit_transform(features)\n",
    "print(f\"Number of molecular descriptors after removing zero-variance features: {features.shape[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "138076f4",
   "metadata": {},
   "source": [
    "### Dataset split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa8675b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X = features\n",
    "# training data size : test data size = 0.8 : 0.2\n",
    "# fixed seed using the random_state parameter, so it always has the same split.\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, train_size=0.8, random_state=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa51e6b",
   "metadata": {},
   "source": [
    "### Data preprocessing\n",
    "\n",
    "As it turns out, different features have different scales and distributions.\\\n",
    "Think of molecular weight, which goes from **0 to around 1000 u.a.m. for small organic molecules**, and electrochemical potential, [typically between -3 and 3](https://par.nsf.gov/servlets/purl/10016877). These differences have a huge impact on ML models, which is the reason why we will normalize all the features.\n",
    "\n",
    "Many classes and functions in the [sklearn.preprocessing](https://scikit-learn.org/stable/modules/preprocessing.html) module can be used for preprocessing data. Here we show a simple Min-Max normalization of features using `MinMaxScaler` in `scikt-learn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b684fd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "scaler = MinMaxScaler()\n",
    "scaler.fit(X_train)\n",
    "\n",
    "# save original X\n",
    "X_train_ori = X_train\n",
    "X_test_ori = X_test\n",
    "# transform data\n",
    "X_train = scaler.transform(X_train)\n",
    "X_test = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e1e3784",
   "metadata": {},
   "source": [
    "### Q: Is there a difference between doing data preprocessing before and after data split? Which one is better and why?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394d369a",
   "metadata": {},
   "source": [
    "### Create models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15b94772",
   "metadata": {},
   "outputs": [],
   "source": [
    "# random forest regressor, and the default criterion is mean squared error (MSE)\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "ranf_reg = RandomForestRegressor(n_estimators=10, random_state=0)  # using 10 trees and seed=0\n",
    "\n",
    "# XGBoost regressor\n",
    "from xgboost import XGBRegressor\n",
    "xgb_reg = XGBRegressor(n_estimators=10, random_state=0)  # using 10 trees and seed=0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b66c0b1",
   "metadata": {},
   "source": [
    "### Train and evaluate the models\n",
    "- Mean Squared Error: $MSE$ = $\\frac{1}{n} \\Sigma_{i=1}^n({y}-\\hat{y})^2$\n",
    "- Root Mean Squared Error: $RMSE$ = $\\sqrt{MSE}$ = $\\sqrt{\\frac{1}{n} \\Sigma_{i=1}^n({y}-\\hat{y})^2}$\n",
    "\n",
    "We choose `RMSE` as the evaluation metric for this task."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68b21ecd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "def train_test_model(model, X_train, y_train, X_test, y_test):\n",
    "    \"\"\"\n",
    "    Function that trains a model, and tests it.\n",
    "    Inputs: sklearn model, train_data, test_data\n",
    "    \"\"\"\n",
    "    # Train model\n",
    "    model.fit(X_train, y_train)\n",
    "    \n",
    "    # Calculate RMSE on training\n",
    "    y_pred_train = model.predict(X_train)\n",
    "    y_pred_test = model.predict(X_test)\n",
    "    model_train_mse = mean_squared_error(y_train, y_pred_train)\n",
    "    model_test_mse = mean_squared_error(y_test, y_pred_test)\n",
    "    model_train_rmse = model_train_mse ** 0.5\n",
    "    model_test_rmse = model_test_mse ** 0.5\n",
    "    print(f\"RMSE on train set: {model_train_rmse:.3f}, and test set: {model_test_rmse:.3f}.\\n\")\n",
    "\n",
    "\n",
    "# Train and test the random forest model\n",
    "print(\"Evaluating Random Forest Model.\")\n",
    "train_test_model(ranf_reg, X_train, y_train, X_test, y_test)\n",
    "\n",
    "# Train and test XGBoost model\n",
    "print(\"Evaluating XGBoost model.\")\n",
    "train_test_model(xgb_reg, X_train, y_train, X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "712c715f",
   "metadata": {},
   "source": [
    "### Q: Compare the results. Which model is better in this case?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "209edf19",
   "metadata": {},
   "source": [
    "**Exercise**: Try to train a [SVM model](https://scikit-learn.org/stable/modules/svm.html#regression) and evaluate it. You can create a [SVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR) model using default parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7746533",
   "metadata": {},
   "outputs": [],
   "source": [
    "### YOUR CODE\n",
    "\n",
    "# import\n",
    "' your code '\n",
    "\n",
    "# create a model\n",
    "svm_reg = ' your code '\n",
    "\n",
    "# train and evaluate the model\n",
    "' your code '\n",
    "\n",
    "### END"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7de3621",
   "metadata": {},
   "source": [
    "Perfect! Now you have mastered training and evaluating the model you want."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "304d42f3",
   "metadata": {},
   "source": [
    "## Cross-validation and hyperparamter optimization\n",
    "\n",
    "Our last topic in this section is hyperparameter optimization.\n",
    "\n",
    "As you've seen, models need some parameters as input, and we need to decide which are the best parameters (e.g. `n_estimators` in tree-based models). Many classes and functions in the [sklearn.model_selection](https://scikit-learn.org/stable/model_selection.html) module can be used for cross validation and hyperparamter optimization.\\\n",
    "Here, we use cross validation and grid search methods to optimize the paramters of random forest model. In cross-validation, the original training data is further split into training and validation sets.\n",
    "\n",
    "The class [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV) does all the work for us, combining cross-validation with grid search, so we can very easily optimize the hyperparameters.\n",
    "\n",
    "<div align=\"center\">\n",
    "<img src=\"https://scikit-learn.org/stable/_images/grid_search_workflow.png\" width=\"400\"/>\n",
    "</div>\n",
    "\n",
    "\n",
    "<font size=5 color=\"red\">\n",
    "Note that when you call grid_search.fit(), you only pass the training data as argument.\n",
    "</font>\n",
    "Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73db8e66-7724-40c4-8876-39668793ff1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "param_grid = {\n",
    "    'n_estimators': [50, 100],\n",
    "    'max_depth': [5, 10, 20, 30]\n",
    "}\n",
    "\n",
    "# use 5-folds cross validation during grid searching\n",
    "grid_search = GridSearchCV(\n",
    "    RandomForestRegressor(random_state=0),\n",
    "    param_grid,\n",
    "    cv=5\n",
    ")\n",
    "grid_search.fit(X_train, y_train)\n",
    "\n",
    "# re-train a model using best hyperparameters\n",
    "rf_gs = RandomForestRegressor(**grid_search.best_params_, random_state=0)\n",
    "\n",
    "print('Best paramters: ', grid_search.best_params_)\n",
    "print('Random forests performance after hyperparamter optimization:')\n",
    "train_test_model(rf_gs, X_train, y_train, X_test, y_test)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "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.9.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "e1b1e114f4dae097b9e32029c5d22d73dc21a5dd723446d46774bd2adced9390"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
