2_Hands_on_ML_SL.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "wTkOlztYCtDJ"
   },
   "source": [
    "# Hands on Machine Learning (ML) and Sequential Learning (SL)\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/paolodeangelis/AEM/blob/main/2_Hands_on_ML_SL.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E0QS0d7u8PJ9"
   },
   "source": [
    "## Preamble\n",
    "\n",
    "In this section we install and import the necessary modules for running the first part of the excersise (importing and handling data, training and validating the ML model, interpreting the ML model); we also define useful functions for materials featurization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8gpYaHBUCtDN"
   },
   "source": [
    "Installing with `pip` the necessary libraries (force the version with ==)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CaWGH7bJ4hxU"
   },
   "outputs": [],
   "source": [
    "%pip install matminer==0.9.2\n",
    "%pip install pymatgen==2024.3.1\n",
    "%pip install scikit_learn==1.3.0\n",
    "%pip install shap==0.45.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AoyeFhdQCtDR"
   },
   "source": [
    "Importing the necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "EtxeRmoT4lVG"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import shap\n",
    "import sklearn\n",
    "from matminer.featurizers import composition as cf\n",
    "from matminer.featurizers.base import MultipleFeaturizer\n",
    "from pymatgen.core import Composition\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_absolute_error, mean_squared_error\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "X0FRvEHJCtDT"
   },
   "source": [
    "Defining functions for efficiently extracting materials features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "jf6A37ww44Dq"
   },
   "outputs": [],
   "source": [
    "def get_compostion(c):\n",
    "    \"\"\"Function to get compositions from chemical formula using pymatgen\"\"\"\n",
    "    try:\n",
    "        return Composition(c)\n",
    "    except:  # noqa: E722\n",
    "        return None\n",
    "\n",
    "\n",
    "def featurizing(data, property_interest=None):\n",
    "    # Featurizer\n",
    "    f = MultipleFeaturizer(\n",
    "        [\n",
    "            cf.Stoichiometry(),\n",
    "            cf.ElementProperty.from_preset(\"magpie\"),\n",
    "            cf.ValenceOrbital(props=[\"avg\"]),\n",
    "            cf.IonProperty(fast=True),\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    # Inputs\n",
    "    data[\"composition\"] = [get_compostion(mat) for mat in data.Components]\n",
    "\n",
    "    featurized_data = pd.DataFrame(\n",
    "        f.featurize_many(data[\"composition\"], ignore_errors=True),\n",
    "        columns=f.feature_labels(),\n",
    "        index=data[\"Components\"],\n",
    "    )\n",
    "    if property_interest:\n",
    "        featurized_data[property_interest] = data[property_interest].values\n",
    "    return featurized_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WJOd_b3Ehbzw"
   },
   "source": [
    "## Data handling\n",
    "\n",
    "In this section we import the Superconductors database (two columns: brute formula and critical Temperature), extracting composition based features; we thus drop rows where at least a NaN (i.e., Not a Number) value occurs. The features represent the *translation* of chemical formulae into numbers, corresponding to the characteristics of such materials; this is necessary, since the ML works with the numbers. It is the equivalent of length and width of sepals and petals in the iris database (already seen in the first introductory lecture to Python), with the only difference that in this case we want to perform a regression of a scalar quantity (the critical Temperature) and not a classification over the type of flower.\n",
    "\n",
    "If you want to know the meaning of those 145 features, please check here https://www.nature.com/articles/npjcompumats201628?report=reader (and in the corresponding supplementary material, pdf). \n",
    "\n",
    "In brief, they are:\n",
    "* **Stoichiometric attributes** that depend only on the fractions of elements present and not what those elements actually are\n",
    "* **Elemental property statistics**, which are defined as the mean, mean absolute deviation, range, minimum, maximum and mode of 22 different elemental properties.\n",
    "* **Electronic structure attributes**, which are the average fraction of electrons from the *s*, *p*, *d* and *f* valence shells between all present elements.\n",
    "* **Ionic compound attributes** that include whether it is possible to form an ionic compound assuming all elements are present in a single oxidation state.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nOOUds-qCtDX"
   },
   "source": [
    "Downloading the database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-xf8zxuKCtDY"
   },
   "outputs": [],
   "source": [
    "# source National Institute of Materials Science, Materials Information Station,\n",
    "# SuperCon, http://supercon.nims.go.jp/index_en.html (2011)\n",
    "\n",
    "!wget https://raw.githubusercontent.com/paolodeangelis/AEM/main/data/Supercon_data_clean.xlsx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "kTutyBek5CeZ"
   },
   "outputs": [],
   "source": [
    "data = pd.read_excel(r\"Supercon_data_clean.xlsx\")  # Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "sDyqH2PoobuY"
   },
   "outputs": [],
   "source": [
    "Featurized_data = featurizing(data, \"Tc\")  # Extract composition based features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "fyi9Z2DUoilV"
   },
   "outputs": [],
   "source": [
    "Featurized_data = Featurized_data.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Iok3sybKAnHW"
   },
   "outputs": [],
   "source": [
    "Featurized_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rv8Yr6gDAF5F"
   },
   "source": [
    "### **Data for the report**\n",
    "Each group will perform the excercise with a different dataset. Your dataset will be obtained by picking 2000 rows from the original one (containing 12914 rows). The selection procedure of those rows is randomic; to ensure the *randomness reproducibility* we set a random state by means of a SEED. The value of the SEED has to be equal to the number of your group.\n",
    "\n",
    "**In the cell below, set the SEED equal to your group number.**\n",
    "\n",
    "The 2000 random rows are picked by setting the variable N_data equal to 2000.\n",
    "\n",
    "**In the cell below, set the N_data equal to 2000.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zx2HmkWdXeSJ"
   },
   "outputs": [],
   "source": [
    "SEED = \n",
    "N_data = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "iyOQdwbD7JsZ"
   },
   "outputs": [],
   "source": [
    "Featurized_data = Featurized_data.sample(n=N_data, random_state=SEED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oVcpdeQm7Vv_"
   },
   "outputs": [],
   "source": [
    "Featurized_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lsstDKLoCNaB"
   },
   "source": [
    "## Machine Learning\n",
    "In this section we train and validate a Random Forest Regressor for constructing a model which predicts the critical temperature of materials, on the basis of their chemical composition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Qcdya8AIC2cW"
   },
   "source": [
    "#### **Training and testing sets**\n",
    "We split the 2000 rows long dataset (obtained in the previous section) in a training set (80% of 2000 rows = 1600 rows) and a testing set (20% of 2000 rows = 400 rows). The split is randomic. To ensure *randomness reproducibility* we set the usual random state. The value of the random state has to be equal to the number of your group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "fwavCrnI5dd0"
   },
   "outputs": [],
   "source": [
    "train_df, test_df = train_test_split(\n",
    "    Featurized_data, test_size=0.2, random_state=SEED\n",
    ")  # split data in training set (80% of the dataset) and testing set (20% of the dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GWKIEdwRD-Gj"
   },
   "source": [
    "#### **Machine Learning model definition**\n",
    "As Machine Learning model for doing regression we choose a Random Forest Regressor.\n",
    "\n",
    "* Since it is affected by randomness, we set the usual random state (it has to be equal to the number of your group).\n",
    "\n",
    "* The only hyperparameter we consider is the number of estimators (number of trees is the forest). It has to be equal to 100\n",
    "\n",
    "**To achieve this, in the cell below set the variable N_trees equal to 100**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "soqdwBC4ZMFZ"
   },
   "outputs": [],
   "source": [
    "N_trees = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "K81kHs0F76Of"
   },
   "outputs": [],
   "source": [
    "rf = RandomForestRegressor(random_state=SEED, n_estimators=N_trees)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JmltcRgsIvBV"
   },
   "source": [
    "#### **Model training**\n",
    "We effectively train the model rf (defined in the previous cell) with the command ```rf.fit(X_training_set, y_training_set)```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2PhT4LXP8IA6"
   },
   "outputs": [],
   "source": [
    "rf.fit(train_df.iloc[:, :-1], train_df.iloc[:, -1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HYukv4bMig_c"
   },
   "source": [
    "#### **Model performances**\n",
    "We compare the real critical temperature of the materials in the testing set, with the corresponding values predicted by the trained model.\n",
    "We check those performances by means of three measures:\n",
    "* coefficient of determination $R^2$, defined as\n",
    "\\begin{equation}\n",
    "R^2(\\mathbf{y}, {\\mathbf{\\hat{y}}}) = 1-\\frac{\\sum_{i=1}^k(y_i-\\hat{y_i})^2}{\\sum_{i=1}^k(y_i-\\overline{y})^2}\n",
    "\\end{equation}\n",
    "where $\\mathbf{y}$ is the $k-$dimensional vector of real values, $\\mathbf{\\hat{y}}$ is the $k-$dimensional vector of predicted values, $\\overline{y}=k^{-1}∑_{i=1}^ky_i$ is the average over real values. $R^2= 1$ if the prediction is perfect (i.e., all the red dots are on the blue line, with $y_i=\\hat{y_i}, \\forall i\\in[1,\\dots, k]$).\n",
    "\n",
    "* Mean Absolute Error $\\mathrm{MAE}$, defined as\n",
    "\\begin{equation}\n",
    "\\mathrm{MAE}(\\mathbf{y}, {\\mathbf{\\hat{y}}}) = \\frac{1}{k}\\sum_{i=1}^k|y_i - \\hat{y_i}|\n",
    "\\end{equation}\n",
    "with the same meaning of the notation.\n",
    "\n",
    "* Root Mean Squared Error $\\mathrm{RMSE}$\n",
    "\\begin{equation}\n",
    "\\mathrm{RMSE}(\\mathbf{y}, {\\mathbf{\\hat{y}}}) = \\left(\\frac{1}{k}\\sum_{i=1}^k(y_i - \\hat{y_i})^2 \\right)^{1/2}\n",
    "\\end{equation}\n",
    "with the same meaning of the notation.\n",
    "\n",
    "In our case, $k$ is the length (i.e., the number of materials) of the testing set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CUvMrcvv8Voc"
   },
   "outputs": [],
   "source": [
    "test_predictions = rf.predict(\n",
    "    test_df.iloc[:, :-1]\n",
    ")  # Predicted y over samples of the testing set\n",
    "test_labels = test_df.iloc[:, -1].values  # True y over samples of the testing set\n",
    "\n",
    "r2 = sklearn.metrics.r2_score(\n",
    "    test_labels, test_predictions\n",
    ")  # coefficient of determination\n",
    "mae = mean_absolute_error(test_labels, test_predictions)  # mean absolute error\n",
    "rmse = np.sqrt(\n",
    "    mean_squared_error(test_labels, test_predictions)\n",
    ")  # root mean squared error\n",
    "delta = max(test_labels) - min(test_labels)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(3, 3), dpi=190)\n",
    "plt.scatter(test_labels, test_predictions, c=\"crimson\", alpha=0.2)\n",
    "p1 = max(max(test_predictions), max(test_labels))\n",
    "p2 = min(min(test_predictions), min(test_labels))\n",
    "plt.plot([p1, p2], [p1, p2], \"b-\")\n",
    "plt.annotate(\n",
    "    \"$R^2$ = %0.3f\" % r2,\n",
    "    xy=(0.02 * delta, 0.95 * delta),\n",
    "    xytext=(0.02 * delta, 0.95 * delta),\n",
    ")\n",
    "plt.annotate(\n",
    "    \"MAE = %0.2f K\" % mae,\n",
    "    xy=(0.02 * delta, 0.85 * delta),\n",
    "    xytext=(0.02 * delta, 0.85 * delta),\n",
    ")\n",
    "plt.annotate(\n",
    "    \"RMSE = %0.2f K\" % rmse,\n",
    "    xy=(0.02 * delta, 0.75 * delta),\n",
    "    xytext=(0.02 * delta, 0.75 * delta),\n",
    ")\n",
    "plt.xlabel(r\"True $T_\\mathrm{c}$ (K)\")\n",
    "plt.ylabel(r\"Predicted $T_\\mathrm{c}$ (K)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J3rpfaXpZeKd"
   },
   "source": [
    "**WARNING**: The figure above could in principle show either very bad performances and very good performances. DON'T WORRY: it's normal, it depends on how much you are lucky in picking the random rows from the original database. You will not get worse grade for bad performances of the model.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aS5lBPy_9oIy"
   },
   "source": [
    "## Interpretability\n",
    "Thanks to the TreeSHAP algorithm, we can find the most relevant features, ranking them in terms of importance with respect to the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "id": "xv7w7oJk9qDN"
   },
   "outputs": [],
   "source": [
    "explainer = shap.TreeExplainer(rf)\n",
    "shap_values = explainer.shap_values(test_df.iloc[:, :-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cnD2_MJdwj7g"
   },
   "source": [
    "#### **Cumulative of normalized feature importances**\n",
    "We produce the cumulative curve of normalized importances. The number of features explaining the 75% of the model is shown. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lD3-Yk7-9xHO"
   },
   "outputs": [],
   "source": [
    "N = np.shape(test_df.iloc[:, :-1])[1]\n",
    "k = 0.75\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "cumsum = np.cumsum(np.sort(np.mean(abs(shap_values), axis=0))[::-1])\n",
    "normalized_cumulative = np.cumsum(np.sort(np.mean(abs(shap_values), axis=0))[::-1]) / (\n",
    "    np.max(np.cumsum(np.sort(np.mean(abs(shap_values), axis=0))[::-1]))\n",
    ")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(3, 3), dpi=190)\n",
    "ax.plot(np.arange(N), normalized_cumulative)\n",
    "ax.plot(np.arange(N), k * np.ones(N))\n",
    "ind_cross1 = np.argmin(\n",
    "    np.fabs(normalized_cumulative - k * max(normalized_cumulative) * np.ones(N))\n",
    ")\n",
    "# plt.yticks(np.array([0, 0.5, 1]))\n",
    "\n",
    "ax.annotate(\n",
    "    \"%i features\" % (ind_cross1 + 1),\n",
    "    xy=(ind_cross1 + 1, 0.01),\n",
    "    xytext=(ind_cross1 + 10, 0.2),\n",
    "    arrowprops=dict(facecolor=\"black\", shrink=0.000005, width=0.1, headwidth=4),\n",
    ")\n",
    "ax.annotate(\n",
    "    \"75% of\\nthe maximum\",\n",
    "    xy=(70, 0.73),\n",
    "    xytext=(80, 0.45),\n",
    "    arrowprops=dict(facecolor=\"black\", shrink=0.0005, width=0.1, headwidth=4),\n",
    ")\n",
    "plt.scatter(ind_cross1, normalized_cumulative[ind_cross1], color=\"orange\")\n",
    "plt.plot(\n",
    "    (ind_cross1, ind_cross1),\n",
    "    (normalized_cumulative[ind_cross1], 0),\n",
    "    color=\"orange\",\n",
    "    ls=\":\",\n",
    ")\n",
    "plt.ylim(0, 1.04)\n",
    "plt.xlabel(\"Features\")\n",
    "plt.ylabel(\"Normalized\\ncumulative importance\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ltl9wJOOaf91"
   },
   "source": [
    "The blue curve here above represents the cumulative curve of coefficients of importances. It takes all those coefficients and sorts them from the biggest to the smallest. For instance, the 15th point of this curve represents the sum of the first 15 coefficients. The overall sum is 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "k0TVTVPAxga2"
   },
   "source": [
    "#### **Ranking of normalized feature importances**\n",
    "We produce the dataset \"Output_mean_shap\" containing the list of features sorted in terms of the corresponding importance. Sum of importances is 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "id": "H6X_BjiG-SHE"
   },
   "outputs": [],
   "source": [
    "Output_shap = pd.DataFrame(\n",
    "    shap_values, index=test_df.iloc[:, :-1].index, columns=test_df.iloc[:, :-1].columns\n",
    ")\n",
    "Output_mean_shap = pd.DataFrame(\n",
    "    abs(Output_shap).describe().loc[\"mean\"]\n",
    "    / sum(abs(Output_shap).describe().loc[\"mean\"])\n",
    ").sort_values(\"mean\", ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "id": "TdH63GtL_sQF"
   },
   "outputs": [],
   "source": [
    "Output_mean_shap.to_excel(\n",
    "    \"Output_mean_shap.xlsx\"\n",
    ")  # list of the features ranked in terms of importance (importances sum up to 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cWMR8AtlCQqI"
   },
   "source": [
    "## Sequential Learning\n",
    "We compare three Sequential Learning acquisition functions to choose the next material to be tested, starting from an initial pool of known materials. \n",
    "\n",
    "As regression methodology over the training set, we consider the Random Forest Regressor by lolopy (which has nothing to do with the random forest regressor used for the predictive model in the previous sections)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5OU7UAQMCSjv"
   },
   "outputs": [],
   "source": [
    "%pip install lolopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "id": "9stSxxOVCWLz"
   },
   "outputs": [],
   "source": [
    "import tqdm\n",
    "from lolopy.learners import RandomForestRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "id": "Tb8wRb_qCavy"
   },
   "outputs": [],
   "source": [
    "def MEI(X: np.ndarray, y: np.ndarray, n_steps: int) -> int:\n",
    "    \"\"\"Acquisition functions MEI.\n",
    "\n",
    "    Args:\n",
    "        X (numpy.ndarray): matrix with n rows (number of total materials for\n",
    "            which doing the SL, in our case 100) and d columns (number of features\n",
    "            taken into account for the optimization)\n",
    "        y (numpy.ndarray): vector with n rows (target property)\n",
    "        n_steps (int): number of steps allowed for doing SL (in our case,\n",
    "            100 total materials - 50 materials in the initial training set\n",
    "            = maximum 50 steps to find the optimum)\n",
    "\n",
    "    Returns:\n",
    "        int: the index of the chosen material\n",
    "    \"\"\"\n",
    "\n",
    "    arr = y\n",
    "    minima = arr.argsort()[0:50]\n",
    "    in_train = np.zeros(len(X), dtype=bool)\n",
    "    in_train[minima] = True\n",
    "\n",
    "    all_inds = set(range(len(y)))\n",
    "    F = np.zeros(n_steps)\n",
    "    G = np.zeros(n_steps)\n",
    "    mei_train = [list(set(np.where(in_train)[0].tolist()))]\n",
    "    mei_train_inds = []\n",
    "\n",
    "    T = 10\n",
    "\n",
    "    for i in tqdm.tqdm(range(n_steps)):\n",
    "        mei_train_inds = mei_train[-1].copy()\n",
    "        mei_search_inds = list(all_inds.difference(mei_train_inds))\n",
    "\n",
    "        mei_selection_index = []\n",
    "        for j in range(T):\n",
    "            model.fit(X[mei_train_inds], y[mei_train_inds])\n",
    "            mei_y_pred_prov = model.predict(X[mei_search_inds])\n",
    "            mei_selection_index.append(np.argmax(mei_y_pred_prov))\n",
    "\n",
    "        mei_index_G = max(set(mei_selection_index), key=mei_selection_index.count)\n",
    "        mei_index = mei_search_inds[mei_index_G]  # Pick the most preferred entry\n",
    "        mei_train_inds.append(mei_search_inds[mei_index_G])\n",
    "        mei_train.append(mei_train_inds)\n",
    "        G[i] = mei_index\n",
    "        F[i] = mei_train_inds[-1]\n",
    "        if mei_train_inds[-1] == np.argmax(y):\n",
    "            break\n",
    "\n",
    "    return F\n",
    "\n",
    "\n",
    "def MLI(X: np.ndarray, y: np.ndarray, n_steps: int) -> int:\n",
    "    \"\"\"Acquisition functions MLI.\n",
    "\n",
    "    Args:\n",
    "        X (numpy.ndarray): matrix with n rows (number of total materials for\n",
    "            which doing the SL, in our case 100) and d columns (number of features\n",
    "            taken into account for the optimization)\n",
    "        y (numpy.ndarray): vector with n rows (target property)\n",
    "        n_steps (int): number of steps allowed for doing SL (in our case,\n",
    "            100 total materials - 50 materials in the initial training set\n",
    "            = maximum 50 steps to find the optimum)\n",
    "\n",
    "    Returns:\n",
    "        int: the index of the chosen material\n",
    "    \"\"\"\n",
    "    arr = y\n",
    "    minima = arr.argsort()[0:50]\n",
    "    in_train = np.zeros(len(X), dtype=bool)\n",
    "    in_train[minima] = True\n",
    "\n",
    "    all_inds = set(range(len(y)))\n",
    "    K = np.zeros(n_steps)\n",
    "    L = np.zeros(n_steps)\n",
    "    mli_train = [list(set(np.where(in_train)[0].tolist()))]\n",
    "    mli_train_inds = []\n",
    "\n",
    "    T = 10\n",
    "\n",
    "    for i in tqdm.tqdm(range(n_steps)):\n",
    "        mli_train_inds = mli_train[-1].copy()\n",
    "        mli_search_inds = list(all_inds.difference(mli_train_inds))\n",
    "\n",
    "        mli_selection_index = []\n",
    "        for j in range(T):\n",
    "            model.fit(X[mli_train_inds], y[mli_train_inds])\n",
    "            mli_y_pred_prov, mli_y_std_prov = model.predict(\n",
    "                X[mli_search_inds], return_std=True\n",
    "            )\n",
    "            mli_selection_index.append(\n",
    "                np.argmax(\n",
    "                    np.divide(\n",
    "                        mli_y_pred_prov - np.max(y[mli_train_inds]), mli_y_std_prov\n",
    "                    )\n",
    "                )\n",
    "            )\n",
    "\n",
    "        mli_index_L = max(set(mli_selection_index), key=mli_selection_index.count)\n",
    "        mli_index = mli_search_inds[mli_index_L]  # Pick the most preferred entry\n",
    "\n",
    "        mli_train_inds.append(mli_search_inds[mli_index_L])\n",
    "        mli_train.append(mli_train_inds)\n",
    "        L[i] = mli_index\n",
    "        K[i] = mli_train_inds[-1]\n",
    "        if mli_train_inds[-1] == np.argmax(y):\n",
    "            break\n",
    "\n",
    "    return K\n",
    "\n",
    "\n",
    "def MU(X: np.ndarray, y: np.ndarray, n_steps: int) -> int:\n",
    "    \"\"\"Acquisition functions MU.\n",
    "\n",
    "    Args:\n",
    "        X (numpy.ndarray): matrix with n rows (number of total materials for\n",
    "            which doing the SL, in our case 100) and d columns (number of features\n",
    "            taken into account for the optimization)\n",
    "        y (numpy.ndarray): vector with n rows (target property)\n",
    "        n_steps (int): number of steps allowed for doing SL (in our case,\n",
    "            100 total materials - 50 materials in the initial training set\n",
    "            = maximum 50 steps to find the optimum)\n",
    "\n",
    "    Returns:\n",
    "        int: the index of the chosen material\n",
    "    \"\"\"\n",
    "\n",
    "    arr = y\n",
    "    minima = arr.argsort()[0:50]\n",
    "    in_train = np.zeros(len(X), dtype=bool)\n",
    "    in_train[minima] = True\n",
    "\n",
    "    all_inds = set(range(len(y)))\n",
    "    R = np.zeros(n_steps)\n",
    "    S = np.zeros(n_steps)\n",
    "    mu_train = [list(set(np.where(in_train)[0].tolist()))]\n",
    "    mu_train_inds = []\n",
    "\n",
    "    T = 10\n",
    "\n",
    "    for i in tqdm.tqdm(range(n_steps)):\n",
    "        mu_train_inds = mu_train[-1].copy()\n",
    "        mu_search_inds = list(all_inds.difference(mu_train_inds))\n",
    "\n",
    "        mu_selection_index = []\n",
    "        for j in range(T):\n",
    "            model.fit(X[mu_train_inds], y[mu_train_inds])\n",
    "            mu_y_pred_prov, mu_y_std_prov = model.predict(\n",
    "                X[mu_search_inds], return_std=True\n",
    "            )\n",
    "            mu_selection_index.append(np.argmax(mu_y_std_prov))\n",
    "\n",
    "        mu_index_R = max(set(mu_selection_index), key=mu_selection_index.count)\n",
    "        mu_index = mu_search_inds[mu_index_R]\n",
    "\n",
    "        mu_train_inds.append(mu_search_inds[mu_index_R])\n",
    "        mu_train.append(mu_train_inds)  # Pick the most preferred entry\n",
    "        R[i] = mu_index\n",
    "        S[i] = mu_train_inds[-1]\n",
    "        if mu_train_inds[-1] == np.argmax(y):\n",
    "            break\n",
    "\n",
    "    return S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "id": "NpTfDvBdCleP"
   },
   "outputs": [],
   "source": [
    "def produce_Data_SL(\n",
    "    Data: pd.DataFrame,\n",
    "    Output_mean_shap: pd.DataFrame,\n",
    "    n_relevant: int,\n",
    "    target_property: str,\n",
    ") -> tuple:\n",
    "    \"\"\"Function to produce datasets for SL starting from the complete database (Featurized_data).\n",
    "\n",
    "    Args:\n",
    "        Data (pandas.DataFrame): complete database (Featureized_data)\n",
    "        Output_mean_shap (pandas.DataFrame): ranking of features used in terms of importance\n",
    "        n_relevant (int): number of relevant features\n",
    "        target_property (str): name of the target property in the complete database\n",
    "\n",
    "    Returns:\n",
    "        tuple: containing:\n",
    "            -  pandas.DataFrame: dataset with only the relevant features + the target property\n",
    "            -  pandas.DataFrame: dataset with relevant features + set of unrelevant features\n",
    "                (summing up to 30 columns) + the target property\n",
    "    \"\"\"\n",
    "\n",
    "    relevant_features = list(Output_mean_shap.iloc[:n_relevant].index)\n",
    "    unrelevant_features = list(\n",
    "        Output_mean_shap.sort_values(\"mean\").iloc[: int(30 - n_relevant)].index\n",
    "    )\n",
    "    all_features = relevant_features + unrelevant_features\n",
    "\n",
    "    relevant_features.append(target_property)\n",
    "    all_features.append(target_property)\n",
    "\n",
    "    Data_sampled = Data.sample(\n",
    "        n=100, random_state=SEED\n",
    "    )  # replace the random_state with your group number\n",
    "\n",
    "    Data_relevant_features = pd.DataFrame(Data_sampled, columns=relevant_features)\n",
    "    Data_all_features = pd.DataFrame(Data_sampled, columns=all_features)\n",
    "\n",
    "    return (Data_relevant_features, Data_all_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5_mmmDpG6du0"
   },
   "source": [
    "### **Definition of the SL predictor**\n",
    "The SL predictor is constructed on the basis of the Random Forest Regressor by lolopy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "id": "Ng-hmO_ODcvC"
   },
   "outputs": [],
   "source": [
    "model = RandomForestRegressor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WX6MDlyD61kc"
   },
   "source": [
    "### **Datasets production for SL**\n",
    "We produce two datasets:\n",
    "* Data_relevant_features, with 100 random rows and the relevant features explaining the 75% of the predictive model (see section Interpretability) + target property column\n",
    "* Data_all_features, with 100 random rows and 30 features, given by the relevant features (above point) + set of unrelevant features + target property column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "id": "jibeNvi3DiDq"
   },
   "outputs": [],
   "source": [
    "Data_relevant_features, Data_all_features = produce_Data_SL(\n",
    "    Featurized_data, Output_mean_shap, ind_cross1 + 1, \"Tc\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jZkCC5ouKAR1"
   },
   "outputs": [],
   "source": [
    "Data_relevant_features.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VTRCRt0xKCb5"
   },
   "outputs": [],
   "source": [
    "Data_all_features.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KGseQBYpFNOJ",
    "outputId": "667c13a1-3f02-417c-dd0a-0273dfdebdc6"
   },
   "source": [
    "![https://raw.githubusercontent.com/paolodeangelis/AEM/main/img/Image_Data.png](https://raw.githubusercontent.com/paolodeangelis/AEM/main/img/Image_Data.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SDGlsl1RFOMz"
   },
   "source": [
    "#### **Sequential Learning relevant features**\n",
    "Start with the pool of the 50 worst samples in terms of the target $y$ (critical Temperature); SL strategies suggest the next material to be evaluated. The objective is to find the optimum among the other 50 materials with as few evaluations as possible. Search stops when the material with the maximum $y$ (among those remaining 50) is chosen. \n",
    "\n",
    "Three strategies are compared: Maximum Expected Improvement (MEI), Maximum Likelihood Improvement (MLI), Maximum Uncertainty (MU)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OQYCshEHCtDq"
   },
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>WARNING: </b> \n",
    "\n",
    "Since also the Random Forest Regressor by lolopy for SL is not deterministic, if you run the code more times, you will end up with different *trajectories* of evaluations.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qWuoV0hjIKyR",
    "outputId": "a3d10649-501b-4f6c-ceff-846fb2a82410"
   },
   "source": [
    "![https://raw.githubusercontent.com/paolodeangelis/AEM/main/img/Image_SL.png](https://raw.githubusercontent.com/paolodeangelis/AEM/main/img/Image_SL.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-fhmSGpOD2x5"
   },
   "outputs": [],
   "source": [
    "MEI_index_relevant = MEI(\n",
    "    Data_relevant_features.iloc[:, :-1].values,\n",
    "    Data_relevant_features.iloc[:, -1].values,\n",
    "    50,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OLg1tkkrECfi"
   },
   "outputs": [],
   "source": [
    "MLI_index_relevant = MLI(\n",
    "    Data_relevant_features.iloc[:, :-1].values,\n",
    "    Data_relevant_features.iloc[:, -1].values,\n",
    "    50,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5MooOKyfEw8U"
   },
   "outputs": [],
   "source": [
    "MU_index_relevant = MU(\n",
    "    Data_relevant_features.iloc[:, :-1].values,\n",
    "    Data_relevant_features.iloc[:, -1].values,\n",
    "    50,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RnCvE4WBbpPq"
   },
   "source": [
    "In the cell below, write the numbers of evaluations performed in the three runs above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JUSAei7Ob2e9"
   },
   "outputs": [],
   "source": [
    "N_MEI_index_relevant = \n",
    "N_MLI_index_relevant = \n",
    "N_MU_index_relevant = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5ZVm3m4mFSjT"
   },
   "source": [
    "#### **Sequential Learning relevant + unrelevant features**\n",
    "Start with the pool of the 50 worst samples in terms of the target $y$ (critical Temperature); SL strategies suggest the next material to be evaluated. The objective is to find the optimum among the other 50 materials with as few evaluations as possible. Search stops when the material with the maximum $y$ (among those remaining 50) is chosen. \n",
    "\n",
    "Three strategies are compared: Maximum Expected Improvement (MEI), Maximum Likelihood Improvement (MLI), Maximum Uncertainty (MU)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fRSu2AL_CtDu"
   },
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>WARNING: </b> \n",
    "\n",
    "Since also the Random Forest Regressor by lolopy for SL is not deterministic, if you run the code more times, you will end up with different *trajectories* of evaluations. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6ruq1WZ7FVdd"
   },
   "outputs": [],
   "source": [
    "MEI_index_all = MEI(\n",
    "    Data_all_features.iloc[:, :-1].values, Data_all_features.iloc[:, -1].values, 50\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "psfm-u6UFkAt"
   },
   "outputs": [],
   "source": [
    "MLI_index_all = MLI(\n",
    "    Data_all_features.iloc[:, :-1].values, Data_all_features.iloc[:, -1].values, 50\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5LhovgrNFpH4"
   },
   "outputs": [],
   "source": [
    "MU_index_all = MU(\n",
    "    Data_all_features.iloc[:, :-1].values, Data_all_features.iloc[:, -1].values, 50\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QJpxThhTcPkF"
   },
   "source": [
    "In the cell below, write the numbers of evaluations performed in the three runs above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HHWxbJ5_cCfB"
   },
   "outputs": [],
   "source": [
    "N_MEI_index_all = \n",
    "N_MLI_index_all = \n",
    "N_MU_index_all = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EbvwqG20Od3B"
   },
   "source": [
    "#### **Comprehensive comparison between more strategies and different sets of features**\n",
    "We plot the performances of the SL in terms of the number of evaluations needed to find the optimum normalized with respect to the average number of evaluations needed with a *naive* random choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EG9-xWdiFt0b"
   },
   "outputs": [],
   "source": [
    "N = 50 / 2  # number of evaluations in random choice\n",
    "labels = [\"MEI\", \"MLI\", \"MU\"]\n",
    "relevant = [\n",
    "    N_MEI_index_relevant / N,\n",
    "    N_MLI_index_relevant / N,\n",
    "    N_MU_index_relevant / N,\n",
    "]  # replace numbers with the numbers of evaluations performed by SL for MEI, MLI, MU with only relevant features\n",
    "all = [\n",
    "    N_MEI_index_all / N,\n",
    "    N_MLI_index_all / N,\n",
    "    N_MU_index_all / N,\n",
    "]  # replace numbers with the numbers of evaluations performed by SL for MEI, MLI, MU with also unrelevant features\n",
    "\n",
    "x = np.arange(len(labels))  # the label locations\n",
    "width = 0.2  # the width of the bars\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 3), dpi=190)\n",
    "\n",
    "rects1 = ax.bar(x - width / 2, relevant, width, label=\"%i features\" % (ind_cross1 + 1))\n",
    "rects2 = ax.bar(x + width / 2, all, width, label=\"30 features\")\n",
    "plt.axhline(y=1, color=\"k\", linewidth=1, linestyle=\"--\")\n",
    "\n",
    "\n",
    "# Add some text for labels, title and custom x-axis tick labels, etc.\n",
    "ax.set_ylabel(\"Number of experiments/\\n number of random experiments\")\n",
    "ax.set_xticks(x)\n",
    "ax.set_xticklabels(labels)\n",
    "ax.legend()\n",
    "\n",
    "ax.annotate(\n",
    "    \"random choice\",\n",
    "    xy=(1, 1),\n",
    "    xytext=(1.3, 1.2),\n",
    "    arrowprops=dict(facecolor=\"black\", width=0.1, headwidth=4),\n",
    ")\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aGFkrohfcfmA"
   },
   "source": [
    "**WARNING**: In the figure here above you could have all the bars over the random choice. DON'T WORRY, you will not get worse grade for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XKFm_Nm5k7tG"
   },
   "source": [
    "## Assignment \n",
    "### Notebook\n",
    "* In this notebook, set the  ```SEED``` equal to the number of your group. Set the ```N_data``` equal 2000 (for the first part of the report). Set ```N_trees``` equal to 100 (for the first part of the report).\n",
    "* In this notebook, after the SL runs, in the subsection \"Comprehensive comparison between more strategies and different sets of features\", please replace the number of evaluations \n",
    "```N_MEI_index_relevant ```,\n",
    "```N_MLI_index_relevant ```,\n",
    "```N_MU_index_relevant ```,\n",
    "```N_MEI_index_all ```,\n",
    "```N_MLI_index_all ```,\n",
    "```N_MU_index_all ```\n",
    "\n",
    "with the ones that you obtain in your runs. Please notice that, since the code is not deterministic, if you run the same code more times, in general you will end up with different results. Don't worry: just report the number of evaluations of one run for each of the 6 cases (3 methodologies with 2 different sets of features).\n",
    "\n",
    "### Report\n",
    "* **First part**: with this notebook modified as prescribed above, obtain main results (model performances, relevant features, Sequential Learning performances).\n",
    "* **Second part**: redo the *Machine Learning* section and the *Interpretability* section modifying the number of materials ```N_data``` in the dataset and the number of estimators ```N_trees``` in the Random Forest Regressor (for the predictive model). Combinations of those parameters are provided for each group in the document \"Assignments for the machine learning lab\". For each combination, please, re-run the code from the SuperCon database import. In this second part DO NOT redo the Sequential Learning.\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "2-Hands_on_ML_SL.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}