{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tOfOEuH5gqx8"
   },
   "source": [
    "# 2025 Exercise 9: Kernels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "fPOidkAYmdXX"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn\n",
    "from sklearn import metrics\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.datasets import make_moons\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WccG7Hgyhbsd"
   },
   "source": [
    "## Exercise 1: Kernel Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ysAkIn0_hgBK"
   },
   "source": [
    "1. Consider an artificial dataset with random input points $X$ uniformly-distributed in the range $[−4, +4]$\n",
    "and ground-truth function $f^*(X) = \\sin(X) + 2 \\cos(4X) + 10 \\sin(X + 1)$. We are going to try to learn\n",
    "this dataset from the observations $y(X) = f^∗(X) + \\eta$, with $\\eta \\sim N(0, 1)$. Plot $n = 50$ training points\n",
    "with their labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "jXY3Ji8ogzwX"
   },
   "outputs": [],
   "source": [
    "n = 50\n",
    "X_train = np.random.rand(n, 1) * 8 - 4\n",
    "\n",
    "# Data for sklearn methods should be of the shape (n_samples, n_features).\n",
    "# Thus, in 1d, (n_samples, 1)\n",
    "\n",
    "def ground_truth(x):\n",
    "    y = np.sin(x) + 2 * np.cos(4 * x) + 10 * np.sin(x + 1)\n",
    "    return y.squeeze()\n",
    "\n",
    "# Labels, instead, should be of shape (n_samples,).\n",
    "# array.squeeze() removes all array's dimensions of length 1.\n",
    "\n",
    "y_train = ground_truth(X_train) + np.random.rand(n)\n",
    "X_test = np.linspace(-4, 4, 100).reshape(-1, 1)\n",
    "y_test = ground_truth(X_test)\n",
    "\n",
    "# You can check the shapes of X's and y's with array.shape!\n",
    "\n",
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hhDA2RM0h61V"
   },
   "source": [
    "2. Try to learn this function using `sklearn.kernel_ridge.KernelRidge` (https://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html) with a linear \"kernel\", and regularization strength $\\alpha = 0.01$. Plot the predictor on the test set and compare with the ground truth.\n",
    "Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lIYdn1woOS1n"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ucxUj27siDKE"
   },
   "source": [
    "3. Repeat substituting the linear kernel with the RBF kernel (*do not use metrics.pairwise.rbf_kernel as it was indicated before, just pass 'rbf' as an argument as shown in the KernelRidge documentation*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nuLsHlYdiLAw"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NAttEdeEiMoF"
   },
   "source": [
    "## Exercise 2: Kernelized Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fxfgx0NRiQLz"
   },
   "source": [
    "What about classification? Previously, in the class, we have seen that a linear model for classification is logistic regression, which solves the following optimization problem $$\\min_w \\sum_{i=1}^n \\log (1 + \\exp (-y_i(X_i \\cdot w))).$$ Can we use it with kernels to classify non-linearly-separable data? The answer is yes, by mapping the data in the feature space, i.e., $X_i \\to \\phi(X_i)$, and applying the representer theorem, i.e., $w = \\sum_{j=1}^n \\alpha_j \\phi(X_j)$, we obtain\n",
    "$$\\min_\\alpha \\sum_{i=1}^n \\log \\left(1 + \\exp \\left(-y_i \\left(\\phi(X_i)\\sum_{j=1}^n \\alpha_j \\phi(X_j)\\right)\\right)\\right),$$\n",
    "$$\\min_\\alpha \\sum_{i=1}^n \\log \\left(1 + \\exp \\left(-y_i \\left(\\sum_{j=1}^n \\alpha_j K(X_i, X_j)\\right)\\right)\\right).$$\n",
    "Therefore, it is sufficient to substitute the data matrix $X$ with the kernel matrix $K_{ij} = K(X_i, X_j)$ in the sklearn algorithm!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RSS7WCTWiwE1"
   },
   "source": [
    "1. Using `sklearn.linear_model.LogisticRegression` (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html), fit the two-moons dataset with the RBF kernel and predict the value of the model on the test points. You can use `metrics.pairwise.rbf_kernel` to compute the kernel and pass the kernel bandwidth $\\gamma=10$ as an argument. Use the code in the notebook to plot the training points and the decision boundary.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dfF7T3eTiPkf"
   },
   "outputs": [],
   "source": [
    "X_train, y_train = make_moons(n_samples=50, noise=0.1)\n",
    "plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5Kb7p5xYohUN"
   },
   "outputs": [],
   "source": [
    "def test_meshgrid(X_train):\n",
    "    '''\n",
    "    This function returns an array of test points X_test and their coordinates\n",
    "    on a meshgrid (x_mesh, y_mesh). numpy.meshgrid creates coordinate values to\n",
    "    construct a 2d grid. We will pass the output X_test to the kernel and use\n",
    "    x_mesh, y_mesh to plot the decision boundary with matplotlib's plt.contourf.\n",
    "    '''\n",
    "\n",
    "    x_min, x_max = X_train[:, 0].min() - 0.5, X_train[:, 0].max() + 0.5\n",
    "    y_min, y_max = X_train[:, 1].min() - 0.5, X_train[:, 1].max() + 0.5\n",
    "    x_mesh, y_mesh = np.meshgrid(\n",
    "        np.arange(x_min, x_max, 0.02), np.arange(y_min, y_max, 0.02))\n",
    "\n",
    "    X_test = np.c_[x_mesh.ravel(), y_mesh.ravel()] # Concatenation along second axis.\n",
    "\n",
    "    return X_test, x_mesh, y_mesh\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "L8EN86jsnobS"
   },
   "outputs": [],
   "source": [
    "X_test, x_mesh, y_mesh = test_meshgrid(X_train)\n",
    "\n",
    "# The best way to understand these 3 arrays is printing their shapes and values!\n",
    "\n",
    "y_pred = ### YOUR CODE ###\n",
    "y_pred = y_pred.reshape(x_mesh.shape)\n",
    "plt.contourf(x_mesh, y_mesh, y_pred, alpha=0.8)\n",
    "plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DZLjF6NMi883"
   },
   "source": [
    "2. Repeat with the Laplacian kernel (`metrics.pairwise.laplacian_kernel`). How do the decision boundaries learned with the two kernels differ? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Z2FGzbini8X5"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IvIMuNPRjCaR"
   },
   "source": [
    "## Exercise 3: Predicting Atomization Energies with Kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1_xNqPVwjCrM"
   },
   "source": [
    "Predicting ground-state molecular properties is a critical task in chemistry. However, solving the Schrodinger equation is computationally intractable even for small systems. Instead, in this homework, you will use a machine learning approach. In particular, you will predict the atomization energies of small organic molecules using kernels. You will use a subset of the QM7 dataset. Each molecule in this dataset is represented by an $N \\times N$ *Coulomb matrix* $M$, where $N$ is the number of atoms in the largest molecule.  The matrix $M$ is specified by the set of nuclear coordinates $\\{R_i\\}_{i\\leq N}$ and the corresponding charges $\\{Z_i\\}_{i \\leq N}$,\n",
    "\\begin{equation}\n",
    "    M_{ij} = \\frac{1}{2} \\, Z_i^{2.4 }\\, \\delta_{ij} + \\frac{Z_i Z_j}{\\|R_i-R_j\\|} \\, (1 - \\delta_{ij}).\n",
    "\\end{equation}\n",
    "If a molecule has less than $N$ atoms, the remaining elements of $M$ are set to 0.\n",
    "Notice that the off-diagonal elements correspond to the Coulomb repulsion between different atoms in the molecule, which gives the name to this molecular representation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "v_wtT_OFqYoO"
   },
   "source": [
    " 1. The dataset is stored as a dictionary-like object in the file `molecules_dataset.npz` and is already divided into train, validation, and test sets, with keys `X_train`, `y_train`, `X_val`, `y_val`, and `X_test`, `y_test` respectively. The numpy arrays $X$'s contain the upper triangular parts of the normalized Coulomb matrices, and the numpy arrays $y$'s the atomization energies of the molecules in units of kcal/mol. Open the dataset using `numpy.load`.\n",
    " ```\n",
    " dataset = numpy.load('molecules_dataset.pt')\n",
    " X_train = dataset['X_train']\n",
    " ...\n",
    " ```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WmoU6zor2z_n"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zmrICfUM3Dx4"
   },
   "source": [
    " 2. Do kernel ridge regression with the Laplacian kernel (https://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html). Set the regularization $\\alpha=10^{-8}$ and select the scale parameter $\\gamma\\in\\{10^{-6},10^{-5},10^{-4},10^{-3},10^{-2},10^{-1}\\}$ resulting in the best Mean Squared Error on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d2BxkAv43pCU"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dXGaN4Vh3FEe"
   },
   "source": [
    " 3. For the selected value of $\\gamma$, plot your predictions on the test set vs the true test labels and compute the root mean square error in kcal/mol.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "T9hcB0Ua3oYX"
   },
   "outputs": [],
   "source": [
    "### YOUR CODE ###"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}