{ "cells": [ { "cell_type": "markdown", "id": "6370a49e", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ddmms/camml-tutorials/blob/main/notebooks/01-intro/tutorial.ipynb)" ] }, { "cell_type": "markdown", "id": "ca147b2d", "metadata": {}, "source": [ "# Tutorial on Basics of Machine Learning in Computational Molecular Science" ] }, { "cell_type": "markdown", "id": "2de94b30", "metadata": {}, "source": [ "## Prelude\n", "\n", "\n", "This tutorial is a basic introduction into machine learning (ML) workflows in the context of molecular simulation and computational materials research. The notebook introduces workflows and concepts related to\n", "\n", "- Data preparation and analysis\n", "- Generating chemical representations and features\n", "- Dimensionality Reduction\n", "- Clustering\n", "- Kernel-based model fitting\n", "- Hyperparameter Optimization\n", "- Uncertainty Quantification\n", "\n", "The required dependencies include\n", "- [Atomic Simulation Environment (ase)](https://wiki.fysik.dtu.dk/ase/): We will use this to store molecular structures and properties\n", "- [scikit-learn](https://scikit-learn.org/stable/): machine learning library\n", "- [dscribe](https://singroup.github.io/dscribe/stable/): library to generate molecular representations (descriptors)\n", "- [openTSNE](https://opentsne.readthedocs.io/en/stable/)\n", "\n", "As you go through the notebook and work through the tasks, look through the documentation pages of those packages if you get stuck.\n", "\n", "We will be working on two datasets of molecules, one that explores composition space (dataset of different molecules), and one that explores configurational space (data from molecular dynamics):\n", "1. Dataset of five short molecular dynamics trajectories of cyclohexane, Courtesy of \"Quantum Chemistry in the Age of Machine Learning\", edited by Pavlo Dral (2022)\n", "2. The QM7 dataset of organic molecules with up to 7 non-hydrogen atoms, Courtesy of Rupp et al., Phys. Rev. Lett. 108, 058301 (2012) and qmlcode.org\n", "\n", "Both datasets feature structures and energies.\n", "\n", "\n", "Reinhard Maurer, University of Warwick (2025)" ] }, { "cell_type": "code", "execution_count": null, "id": "78e3c6c3", "metadata": {}, "outputs": [], "source": [ "# only if you run in google colab\n", "\n", "# import locale\n", "# locale.getpreferredencoding = lambda: \"UTF-8\"\n", "\n", "# ! pip install ase scikit-learn dcribe opentsne data-tutorials\n", "# get_ipython().kernel.do_shutdown(restart=True)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7c2e781b", "metadata": {}, "outputs": [], "source": [ "# get the data \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "276949ab", "metadata": {}, "outputs": [], "source": [ "#basic stuff\n", "import os\n", "import sys\n", "from functools import partial\n", "import numpy as np\n", "import ase\n", "import random\n", "from ase.io import read, write\n", "from ase.visualize import view\n", "from matplotlib import pyplot as plt\n", "import matplotlib as mpl\n", "import pandas\n", "import seaborn as sns\n", "from ase.build import molecule\n", "from weas_widget import WeasWidget\n", "\n", "#ML stuff\n", "import dscribe\n", "import sklearn\n", "from sklearn.metrics.pairwise import pairwise_kernels\n", "from sklearn.model_selection import train_test_split\n", "from tqdm.auto import tqdm # progress bars for loops\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "c0d9bd4d", "metadata": {}, "source": [ "## Part 1: Data Preparation and Analysis" ] }, { "cell_type": "markdown", "id": "13f16b8b", "metadata": {}, "source": [ "Let's start with a dataset of five short MD simulations of different conformers of cyclohexane.\n", "\n", "We start five independent simulations\n", "initialized within each of the known cyclohexane conformers shown in the below figure. The MD simulations were run for 10,000 time\n", "steps each. The data contains the atom positions, velocities, energies and forces. These types of MD simulations will explore the energy landscape and settle within\n", "local and/or global minima.\n", "\n", "![cyclehexane conformers](./conformers.png)" ] }, { "cell_type": "markdown", "id": "14247a58", "metadata": {}, "source": [ "### Read and analyse Molecular Dynamics Data" ] }, { "cell_type": "code", "execution_count": null, "id": "6e10ccd5", "metadata": {}, "outputs": [], "source": [ "# read in the frames from each MD simulation\n", "traj = []\n", "names = ['chair', 'twist-boat', 'boat', 'half-chair', 'planar']\n", "rgb_colors = [(0.13333333333333333, 0.47058823529411764, 0.7098039215686275),\n", " (0.4588235294117647, 0.7568627450980392, 0.34901960784313724),\n", " (0.803921568627451, 0.6078431372549019, 0.16862745098039217),\n", " (0.803921568627451, 0.13725490196078433, 0.15294117647058825),\n", " (0.4392156862745098, 0.2784313725490196, 0.611764705882353),]\n", "\n", "ranges = np.zeros((len(names), 2), dtype=int)\n", "conf_idx = np.zeros(len(names), dtype=int)\n", "\n", "for i, n in enumerate(names):\n", " frames = read(f'./cyclohexane_data/MD/{n}.xyz', '::')\n", "\n", " for frame in frames:\n", " # wrap each frame in its box\n", " frame.wrap(eps=1E-10)\n", "\n", " # mask each frame so that descriptors are only centered on carbon (#6) atoms\n", " mask = np.zeros(len(frame))\n", " mask[np.where(frame.numbers == 6)[0]] = 1\n", " frame.arrays['center_atoms_mask'] = mask\n", "\n", " ranges[i] = (len(traj), len(traj) + len(frames)) #list of data ranges\n", " conf_idx[i] = len(traj) # handy list to indicate the index of the first frame for each trajectory\n", " traj = [*traj, *frames] # full list of frames, 50000 entries" ] }, { "cell_type": "markdown", "id": "3450f70f", "metadata": {}, "source": [ "After this call, all the MD frames of the 5 different runs sit in the object `traj`. \n", "`traj` is a list of ASE `Atoms` objects, which contain the structure definition, positions, velocities, energies etc.\n", "\n", "If you are unfamiliar with ASE, please explore the online documentation and play around a bit with it in the following cells.\n", "\n", "You can see, we have 50,000 configurations, each labeled with an energy in electronvolt (eV)." ] }, { "cell_type": "markdown", "id": "fa829141", "metadata": {}, "source": [ "**Task for you**\n", "
\n", "\n", "- Visualize the trajectories of the 5 different molecules and study their different configurations and dynamics.\n", " \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "81d9e5cb", "metadata": {}, "outputs": [], "source": [ "#pick a molecule to visualize its trajectory\n", "# ['chair', 'twist-boat', 'boat', 'half-chair', 'planar']\n", "molecule = 'chair'\n", "########\n", "i = names.index(molecule)\n", "r = ranges[i]\n", "mol_traj = traj[r[0]:r[1]]\n", "\n", "print(\"Visualizing '{0}' trajectory\".format(molecule))\n", "# visualize the trajectory of the molecule\n", "viewer = WeasWidget()\n", "viewer.from_ase(mol_traj)\n", "viewer.avr.model_style = 1\n", "viewer.avr.show_hydrogen_bonds = True\n", "viewer\n" ] }, { "cell_type": "markdown", "id": "2c347238", "metadata": {}, "source": [ "Next, we plot the energy as a function of time along the trajectory." ] }, { "cell_type": "code", "execution_count": null, "id": "def3caeb", "metadata": {}, "outputs": [], "source": [ "# energies of the simulation frames\n", "energy = np.array([a.info['energy_eV'] for a in traj])\n", "\n", "# energies of the known conformers\n", "c_energy = np.array([traj[c].info['energy_eV'] for c in conf_idx])\n", "\n", "# extrema for the energies\n", "max_e = max(energy)\n", "min_e = min(energy)\n", "\n", "print('energy range goes from {0:10.6f} to {1:10.6f} eV'.format(min_e, max_e))" ] }, { "cell_type": "code", "execution_count": null, "id": "ed2606ac", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(3*4.8528, 3*1.2219))\n", "\n", "for n, c, r, rgb in zip(names, c_energy, ranges, rgb_colors):\n", " ax.plot(range(0, r[1] - r[0]),\n", " energy[r[0]:r[1]] - min_e,\n", " label=n,\n", " c=rgb,\n", " zorder=-1)\n", " \n", "ax.legend()\n", "ax.set_xlabel(\"Simulation Timestep\")\n", "ax.set_ylabel(\"Energy in eV\")\n", "\n", "ax.set_xlim([0, len(energy)//5])\n", "ax.set_ylim([-0.1, 1.25 * (max_e - min_e)])\n", "ax.set_yticklabels([])\n", "\n", "plt.tight_layout()\n", "plt.savefig('energy.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5bd805f7", "metadata": {}, "source": [ "We can see that most trajectories appear relatively stable in their configuration, but the 'planar' configuration seems to change its configuration relatively quickly to one of the others." ] }, { "cell_type": "markdown", "id": "bf4c47e3", "metadata": {}, "source": [ "\n", "**Task for you**\n", "
\n", " \n", "- The five simulations convergence into two different conformers, the chair conformer and the twist-boat conformer. Can you tell which simulation converges into which configuration? (This is doable the old-fashioned way by visually inspecting the MD simulations or by using machine learning methods.)\n", "\n", "- Make a note of this as it will help you to rationalise some of the later findings.\n", "
\n" ] }, { "cell_type": "markdown", "id": "8450b8d0", "metadata": {}, "source": [ "### Read and analyse QM7 molecule dataset" ] }, { "cell_type": "markdown", "id": "6befb02b", "metadata": {}, "source": [ "The QM7 dataset contains 7101 molecules with up to 7 non-hydrogen atoms. It is one of the first datasets that was studied when the idea of \"quantum machine learning\" came up, so machine learning regression of quantum chemistry properties based on feature representations that correlate well with molecular structure and properties.\n", "\n", "For each of the 7101 molecules, we have 2 labels:\n", "- PBE0/def2-TZVP atomization energy in eV (`energy_eV`)\n", "- and the difference between PBE0 and DFTB3 atomization energy in eV (`delta_energy_eV`)\n", "\n", "Reading the database can take up to two minutes, so please be patient." ] }, { "cell_type": "code", "execution_count": null, "id": "61665806", "metadata": {}, "outputs": [], "source": [ "def get_qm7_energies(filename, key=\"dft\"):\n", " \"\"\" Returns a dictionary with heats of formation for each xyz-file.\n", " \"\"\"\n", "\n", " f = open(filename, \"r\")\n", " lines = f.readlines()\n", " f.close()\n", "\n", " energies = []\n", "\n", " for line in lines:\n", " tokens = line.split()\n", "\n", " xyz_name = tokens[0]\n", " hof = float(tokens[1])\n", " dftb = float(tokens[2])\n", "\n", " if key==\"delta\":\n", " energies.append(hof - dftb)\n", " else:\n", " energies.append(hof)\n", "\n", " return energies\n" ] }, { "cell_type": "code", "execution_count": null, "id": "53d9b802", "metadata": {}, "outputs": [], "source": [ "# Import QM7, already parsed to QML\n", "from ase.io import read\n", "\n", "qm7_dft_energy = get_qm7_energies(\"hof_qm7.txt\", key=\"dft\")\n", "qm7_delta_energy = get_qm7_energies(\"hof_qm7.txt\", key = \"delta\")\n", "\n", "qm7 = [read(\"qm7/\"+f) for f in sorted(os.listdir(\"qm7/\"))]\n", "\n" ] }, { "cell_type": "markdown", "id": "f4e216bf", "metadata": {}, "source": [ "Let's visualize some of those more than 7000 molecules. You will quickly see that the database also contains some rather strange ones ..." ] }, { "cell_type": "code", "execution_count": null, "id": "40234a44", "metadata": {}, "outputs": [], "source": [ "mol = qm7[21]\n", "print(mol)\n", "viewer = WeasWidget()\n", "viewer.from_ase(mol)\n", "viewer.avr.model_style = 1\n", "viewer.avr.show_hydrogen_bonds = True\n", "viewer" ] }, { "cell_type": "code", "execution_count": null, "id": "1bae28bf", "metadata": {}, "outputs": [], "source": [ "#lets attach the labels to the ASE objects as we did in the previous dataset\n", "for i, mol in enumerate(qm7):\n", " mol.info[\"energy_eV\"] = qm7_dft_energy[i] # PBE0/def2-TZVP atomization energy in eV\n", " mol.info[\"energy_delta_eV\"] = qm7_delta_energy[i] # difference between PBE0 and DFTB3 atomization energy in eV\n", "\n", "random.seed(42)\n", "random.shuffle(qm7)\n" ] }, { "cell_type": "markdown", "id": "fc9f827b", "metadata": {}, "source": [ "Let's analyse this data.\n", "\n", "Let's first look at the molecular weight distribution and the elemental distribution" ] }, { "cell_type": "code", "execution_count": null, "id": "8c1231f6", "metadata": {}, "outputs": [], "source": [ "\n", "mol_weights = np.array([a.get_masses().sum() for a in qm7], dtype=np.float64)\n", "number_of_atoms = np.array([len(a) for a in qm7], dtype=np.float64)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "daad5db5", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import scipy.stats as st\n", "\n", "### plotting histograms\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 4)) # 1 row, 2 columns\n", "data1, data2 = mol_weights, number_of_atoms\n", "\n", "#Freedman-Diaconis rule to find optimal number of bins for histogram\n", "q25, q75 = np.percentile(data1, [25, 75])\n", "bin_width = 2 * (q75 - q25) * len(data1) ** (-1/3)\n", "bins = round((data1.max() - data1.min()) / bin_width)\n", "print(\"Freedman–Diaconis number of bins for molecular weight:\", bins)\n", "\n", "axes[0].hist(data1, bins=bins, density=True, alpha=0.7)#, histtype='step')\n", "mn, mx = axes[0].get_xlim()\n", "kde_xs = np.linspace(mn, mx, 300)\n", "kde = st.gaussian_kde(data1)\n", "axes[0].plot(kde_xs, kde.pdf(kde_xs), label=\"PDF\")\n", "axes[0].set_title('Molecular weight distribution')\n", "axes[0].set_xlabel('Molecular weight')\n", "axes[0].set_ylabel('Frequency')\n", "\n", "q25, q75 = np.percentile(data2, [25, 75])\n", "bin_width = 2 * (q75 - q25) * len(data2) ** (-1/3)\n", "bins = round((data2.max() - data2.min()) / bin_width)\n", "print(\"Freedman–Diaconis number of bins for molecule size:\", bins)\n", "\n", "# Plot the second histogram\n", "axes[1].hist(data2, bins=bins, density=True, alpha=0.7)#, histtype='step')\n", "mn, mx = axes[1].get_xlim()\n", "kde_xs = np.linspace(mn, mx, 300)\n", "kde = st.gaussian_kde(data2)\n", "axes[1].plot(kde_xs, kde.pdf(kde_xs), label=\"PDF\")\n", "axes[1].set_title('Molecule size distribution')\n", "axes[1].set_xlabel('Molecule size')\n", "axes[1].set_ylabel('Frequency')\n", "\n", "# Adjust layout and show the plot\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6b04576b", "metadata": {}, "source": [ "\n", "**Task for you**\n", "
\n", " \n", "- Plot the PBE0 and PBE0-DFTB3 \"delta\" energy distributions as histograms. \n", "- Do molecular weight or molecule size correlate with the PBE0 atomization energy?\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f7f7bf27", "metadata": {}, "outputs": [], "source": [ "edata = np.array(qm7_dft_energy)\n" ] }, { "cell_type": "markdown", "id": "2cc7a5f1", "metadata": {}, "source": [ "### Data Prep Summary" ] }, { "cell_type": "markdown", "id": "0c661234", "metadata": {}, "source": [ "Now we basically understand our datasets and have them set up. \n", "\n", "- The cyclohexane MD data arrays are\n", "`names` (list of 5 conformation names), `ranges` (list of data ranges), `conf_idx` (handy list to indicate the index of the first frame for each trajectory), `traj` (full list of 50,000 structures)\n", "\n", "- The QM7 data sits in `qm7`" ] }, { "cell_type": "markdown", "id": "9a37bc79", "metadata": {}, "source": [ "## Part 2: Generate Descriptors\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "93389931", "metadata": {}, "source": [ "We will now play around with two different examples of descriptors to capture the key features of the dynamics. First, we will look at global descriptors and then at an atom-centered descriptor.\n", "\n", "This will include\n", "- Coulomb Matrix (global)\n", "- MBTR (global)\n", "- SOAP (atom-centered and global)" ] }, { "cell_type": "markdown", "id": "b0837dcf", "metadata": {}, "source": [ "### Coulomb Matrix" ] }, { "cell_type": "markdown", "id": "29972108", "metadata": {}, "source": [ "\n", "Coulomb Matrix (CM) is a simple global descriptor which mimics the\n", "electrostatic interaction between nuclei.\n", "\n", "Coulomb matrix is calculated with the equation below.\n", "\n", "$$\n", " \\begin{equation}\n", " M_{ij}^\\mathrm{Coulomb}=\\left\\{\n", " \\begin{matrix}\n", " 0.5 Z_i^{2.4} & \\text{for } i = j \\\\\n", " \\frac{Z_i Z_j}{R_{ij}} & \\text{for } i \\neq j\n", " \\end{matrix}\n", " \\right.\n", " \\end{equation}\n", "$$\n", "\n", "The diagonal elements can be seen as the interaction of an atom with itself. The strange exponent arises from an attempt to correlate the trend of the atomic energy w.r.t $Z_i$. If the diagonal elements can closely capture that trend, it goes along way in correctly capturing molecular stability. The off-diagonal elements represent the Coulomb repulsion between nuclei $i$ and $j$. Therefore, all elements of the representation are physically motivated." ] }, { "cell_type": "code", "execution_count": null, "id": "d48faf80", "metadata": {}, "outputs": [], "source": [ "import collections.abc\n", "collections.Iterable = collections.abc.Iterable # this is unfortunately necessary for compatbility of describe with Py>=3.10\n", "from dscribe.descriptors import CoulombMatrix\n", "\n", "atomic_numbers = [1, 6] # H and C atoms\n", "\n", "# Setting up the CM descriptor\n", "cm = CoulombMatrix(\n", " n_atoms_max=18, # maximum no. atoms in the studied molecules\n", " permutation = 'eigenspectrum',\n", " #permutation=\"sorted_l2\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "46a50ca0", "metadata": {}, "outputs": [], "source": [ "# Create CM output for the system\n", "coulomb_matrices =cm.create(traj)\n", "print(\"flattened\", coulomb_matrices.shape)\n", "print('50000 MD frames, each has a descriptor of length 324 (18x18 matrix flattened)')" ] }, { "cell_type": "markdown", "id": "910c7044", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- Take a look at one of the CM descriptors for a given frame. Compare to the list of atoms. Does it make sense? \n", "- Why does the descriptor have the length that it has?\n", " \n", "
\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "332ed32d", "metadata": {}, "outputs": [], "source": [ "def matprint(mat, fmt=\"g\", size=18):\n", " if len(mat.shape)>1:\n", " #mat = mat.reshape([size,size])\n", " col_maxes = [max([len((\"{:\"+fmt+\"}\").format(x)) for x in col]) for col in mat.T]\n", " for x in mat:\n", " for i, y in enumerate(x):\n", " print((\"{:\"+str(col_maxes[i])+fmt+\"}\").format(y), end=\" \")\n", " print(\"\")\n", " else:\n", " words = [\"{0:8.4f}\".format(x) for x in mat]\n", " print(\" \".join(words))\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "d0131da4", "metadata": {}, "source": [ "\n", "\n", "A good descriptor should be invariant with respect to atom permutations, translations, and rigid rotations. However, for the CM, it matters in which order the atoms appear. Changing of atom ordering leads to swapping of rows and columns in the matrix. This should raise alarm bells as it means that it isn't invariant to atom index permutations and therefore not suitable for general machine learning." ] }, { "cell_type": "markdown", "id": "e293a40e", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- The CM construction above allows for multiple different options for `permutation`. Read the DScribe documentation and try different ones. Which one satisfies permutational invariance? This is the one we should be using.\n", "
\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "78b6bbee", "metadata": {}, "outputs": [], "source": [ "#Pick first molecule frame\n", "first_frame = frames[0]\n", "cm_original = cm.create(first_frame)\n", "print('Original CM')\n", "#print(cm_original)\n", "matprint(cm_original)\n", "\n", "# Translation\n", "first_frame.translate((5, 7, 9))\n", "cm_translated = cm.create(first_frame)\n", "print('Translated CM')\n", "matprint(cm_translated)\n", "\n", "# Rotation\n", "first_frame.rotate(90, 'z', center=(0, 0, 0))\n", "cm_rotated = cm.create(first_frame)\n", "print(\"Rotated CM\")\n", "matprint(cm_rotated)\n", "\n", "# Permutation\n", "upside_down = first_frame[::-1]\n", "cm_upside_down = cm.create(upside_down)\n", "print(\"upside down CM\")\n", "matprint(cm_upside_down)\n", "\n", "print('Differences due to translation')\n", "trans_diff = (cm_translated-cm_original)\n", "matprint(trans_diff)\n", "print('Differences due to rotation')\n", "rot_diff = (cm_rotated-cm_original)\n", "matprint(rot_diff)\n", "print('Differences due to permutation')\n", "permut_diff = (cm_upside_down-cm_original)\n", "matprint(permut_diff)" ] }, { "cell_type": "markdown", "id": "db0a489a", "metadata": {}, "source": [ "The benefit of a global descriptor is that it always has the same length, so we can directly compare between any frame. The downside is that it scales badly with system size.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dab429b5", "metadata": {}, "outputs": [], "source": [ "counts, bins = np.histogram(coulomb_matrices.flatten(),bins=200)\n", "#plt.stairs(counts, bins)\n", "plt.hist(bins[:-1], bins, weights=counts)\n", "plt.xlabel('CM feature value')\n", "plt.ylabel('log(count)')\n", "plt.yscale('log')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "806f724a", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- How does the CM construction work for the QM7 dataset where different molecules have different numbers of atoms? Generate CM features for the QM7 dataset and plot a few of them as a histogram\n", "\n", "(Hint: You will need to figure out what elements are contained in the database and what the maximum number of atoms in the largest molecule is.)\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "id": "167ff348", "metadata": {}, "outputs": [], "source": [ "# space to play around\n", "\n", "import collections.abc\n", "collections.Iterable = collections.abc.Iterable # this is unfortunately necessary for compatbility of describe with Py>=3.10\n", "from dscribe.descriptors import CoulombMatrix\n", "\n", "\n", "max_atoms_qm7 = int(number_of_atoms.max())\n", "\n", "atomic_numbers = []\n", "for mol in qm7:\n", " atomic_numbers.extend(mol.get_atomic_numbers())\n", "\n", "atomic_numbers = np.unique(np.array(atomic_numbers))\n", "print(atomic_numbers)\n", "\n", "species = []\n", "for mol in qm7:\n", " species.extend(mol.get_chemical_symbols())\n", "\n", "species = np.unique(np.array(species))\n", "print(species)\n", "\n", "# Setting up the CM descriptor\n", "cm_qm7 = CoulombMatrix(\n", " n_atoms_max=max_atoms_qm7, # maximum no. atoms in the studied molecules\n", " #permutation = 'eigenspectrum',\n", " permutation=\"sorted_l2\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "e6f59bc4", "metadata": {}, "outputs": [], "source": [ "# Create CM output for the system\n", "coulomb_matrices_qm7 =cm_qm7.create(qm7)\n", "print(\"flattened\", coulomb_matrices_qm7.shape)\n", "\n", "\n", "counts, bins = np.histogram(coulomb_matrices_qm7.flatten(),bins=200)\n", "#plt.stairs(counts, bins)\n", "plt.hist(bins[:-1], bins, weights=counts)\n", "plt.xlabel('CM feature value')\n", "plt.ylabel('log(count)')\n", "plt.yscale('log')\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "2ca738fc", "metadata": {}, "source": [ "### MBTR" ] }, { "cell_type": "markdown", "id": "85349c1a", "metadata": {}, "source": [ "The many-body tensor representation (MBTR) encodes a structure by using a distribution of different structural motifs. It can be used directly for both finite and periodic systems. MBTR is especially suitable for applications where interpretability of the input is important because the features can be easily visualized and they correspond to specific structural properties of the system (distances, angles, dihedrals).\n", "\n", "The MBTR representation can be used as a global representation (one feature vector per system) or as an atom-centered or local one (one feature vector per atom). As a global descriptor, the individual contributions are discretised on a single spectrum for the molecule as shown in the below figure." ] }, { "cell_type": "markdown", "id": "e8eb8d8a", "metadata": {}, "source": [ "![bla](https://singroup.github.io/dscribe/latest/_images/mbtr.jpg)" ] }, { "cell_type": "code", "execution_count": null, "id": "29c79559", "metadata": {}, "outputs": [], "source": [ "import collections.abc\n", "collections.Iterable = collections.abc.Iterable # this is unfortunately necessary for compatbility of describe with Py>=3.10\n", "from dscribe.descriptors import MBTR\n", "\n", "min_dist = 0.8 #smallest distance in Angstrom\n", "max_dist = 12.0 #largest distance in Angstrom\n", "#we are constructing MBTR with inverse distances\n", "grid_min = 1.0/max_dist # 1/Angstrom\n", "grid_max = 1.0/min_dist # 1/Angstrom\n", "grid_n = 100 # number of grid points\n", "sigma = 0.05\n", "\n", "#k=1, 1-body\n", "#\"atomic_number\": The atomic number.\n", "#k=2, 2-body\n", "#\"distance\": Pairwise distance in angstroms.\n", "#\"inverse_distance\": Pairwise inverse distance in 1/angstrom.\n", "#k=3, 3-body\n", "#\"angle\": Angle in degrees.\n", "#\"cosine\": Cosine of the angle.\n", "\n", "# MBTR setup for cyclohexane dataset\n", "mbtr = MBTR(\n", " species=[\"H\", \"C\"],\n", " geometry={\"function\": \"inverse_distance\"}, # only 2-body terms, no 3-body\n", " grid={\"min\": grid_min, \"max\": grid_max, \"n\": grid_n, \"sigma\": sigma},\n", " weighting={\"function\": \"exp\", \"scale\": 0.5, \"threshold\": 1e-3},\n", " periodic=False,\n", " normalization=\"l2\",\n", " sparse=False\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "16c8c8a2", "metadata": {}, "outputs": [], "source": [ "mbtrs = mbtr.create(traj)\n", "#mbtrs = mbtr.create(traj[:10])\n", "print(mbtrs.shape)" ] }, { "cell_type": "code", "execution_count": null, "id": "9e2752e2", "metadata": {}, "outputs": [], "source": [ "plt.plot(mbtrs[0])\n", "plt.title('Full MBTR vector for one molecule')\n", "\n", "\n", "# Create the mapping between an index in the output and the corresponding\n", "# chemical symbol\n", "n_elements = len(mbtr.species)\n", "x = np.linspace(grid_min, grid_max, grid_n)\n", "\n", "# Plot k=2\n", "if mbtr.k==2:\n", " fig, ax = plt.subplots()\n", " for i in range(n_elements):\n", " for j in range(n_elements):\n", " if j >= i:\n", " i_species = mbtr.species[i]\n", " j_species = mbtr.species[j]\n", " loc = mbtr.get_location((i_species, j_species))\n", " plt.plot(x, mbtrs[0][loc], label=\"{}-{}\".format(i_species, j_species))\n", " if mbtr.geometry['function'] == 'inverse_distance':\n", " ax.set_xlabel(\"Inverse distance (1/angstrom)\")\n", " else:\n", " ax.set_xlabel(\"distance (angstrom)\")\n", " plt.title('Individual components')\n", "\n", "#plot k=3\n", "if mbtr.k==3:\n", " fig, ax = plt.subplots()\n", " for i in range(n_elements):\n", " for j in range(n_elements):\n", " for k in range(n_elements):\n", " i_species = mbtr.species[i]\n", " j_species = mbtr.species[j]\n", " k_species = mbtr.species[k]\n", " loc = mbtr.get_location((i_species, j_species, k_species))\n", " plt.plot(x, mbtrs[0][loc], label=\"{}-{}-{}\".format(i_species, j_species, k_species))\n", " if mbtr.geometry['function'] == 'angle':\n", " ax.set_xlabel(\"angle in rad\")\n", " else:\n", " ax.set_xlabel(\"cos(angle)\")\n", "\n", " plt.title('Individual components')\n", "\n", "\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "e93ab896", "metadata": {}, "source": [ "**Task for you**\n", "
\n", "\n", "\n", "- Can you explain why the MBTR descriptor has the length it does? \n", "- Is this the optimal distance range to be considered? What does the value of sigma do to the descriptor?\n", "- What impact does it have to use 1-body, 2-body, 3-body components? How large is the descriptor in each case?\n", "- Play around with MBTR parameters a bit to find the optimal settings, then revert to the original settings provided.\n", "\n", "- Does MBTR satisfy all our symmetry requirements?\n", "\n", "(Hint: pick a small set of geometries traj[0:10] to play around with settings to reduce computing time)\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "cc06aff7", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- Set up an MBTR feature representation for the QM7 dataset and play around with the parameters\n", "(Make sure to give the MBTR object a different name, e.g. mbtr_qm7)\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "id": "301176a3", "metadata": {}, "outputs": [], "source": [ "# space to play around" ] }, { "cell_type": "code", "execution_count": null, "id": "670f3d4c", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "3d14cce4", "metadata": {}, "source": [ "### SOAP" ] }, { "cell_type": "markdown", "id": "caf4ba34", "metadata": {}, "source": [ "Smooth Overlap of Atomic Positions (SOAP) is a descriptor that encodes atomic environments by using a local expansion of a gaussian smeared atomic density with orthonormal functions based on spherical harmonics and radial basis functions.\n", "\n", "$$\n", " p^{Z_1 Z_2}_{n n' l} = \\pi \\sqrt{\\frac{8}{2l+1}}\\sum_m {c^{Z_1}_{n l m}}^*c^{Z_2}_{n' l m}\n", "$$\n", "\n", "here $n$ and $n'$ are indices for the different radial basis\n", "functions up to $n_\\mathrm{max}$, $l$ is the angular degree of the\n", "spherical harmonics up to $l_\\mathrm{max}$ and $Z_1$ and $Z_2$\n", "are atomic species.\n", "\n", "The coefficients $c^Z_{nlm}$ are defined as the following inner\n", "products:\n", "\n", "$$\n", " c^Z_{nlm} =\\iiint_{\\mathcal{R}^3}\\mathrm{d}V g_{n}(r)Y_{lm}(\\theta, \\phi)\\rho^Z(\\mathbf{r}).\n", "$$\n", "\n", "where $\\rho^Z(\\mathbf{r})$ is the gaussian smoothed atomic density for\n", "atoms with atomic number $Z$ defined as\n", "\n", "$$\n", " \\rho^Z(\\mathbf{r}) = \\sum_i^{\\lvert Z_i \\rvert} e^{-1/2\\sigma^2 \\lvert \\mathbf{r} - \\mathbf{R}_i \\rvert^2}\n", "$$\n", "\n", "$Y_{lm}(\\theta, \\phi)$ are the real spherical harmonics, and\n", "$g_{n}(r)$ is the radial basis function.\n", "\n", "For the radial degree of freedom, $g_{n}(r)$, multiple approaches may be used. By\n", "default the DScribe implementation uses spherical gaussian type orbitals as\n", "radial basis functions, as they allow faster analytic computation.\n", "\n", "The SOAP similarity kernel between two atomic environments can be retrieved\n", "as a normalized polynomial kernel of the partial power spectrums:\n", "$$\n", " K^\\mathrm{SOAP}(\\mathbf{p}, \\mathbf{p'}) = \\left( \\frac{\\mathbf{p} \\cdot \\mathbf{p'}}{\\sqrt{\\mathbf{p} \\cdot \\mathbf{p}~\\mathbf{p'} \\cdot \\mathbf{p'}}}\\right)^{\\xi}\n", "$$\n", "\n", "Starting point for parameters:\n", "- `r_cut=3.5`: we are considering 3.5A radius of sphere for describing each atom environments. For cyclohexane, this includes the whole ring, which is necessary to differentiate conformers. For QM7, this will likely be quite different.\n", "- `n_max=4`: expand over 4 radial GTO bases\n", "- `l_max=4`: expand over the first 4 spherical harmonics\n", "- `sigma=0.3`: assume each atom has a gaussian of width sigma=0.3 imposed on its lattice site\n", "\n", "SOAP is a descriptor centred around atoms. For the construction of machine learning interatomic potentials, this is beneficial. " ] }, { "cell_type": "code", "execution_count": null, "id": "6f2047bf", "metadata": {}, "outputs": [], "source": [ "import collections.abc\n", "collections.Iterable = collections.abc.Iterable # this is unfortunately necessary for compatbility of describe with Py>=3.10\n", "from dscribe.descriptors import SOAP\n", "\n", "species = [\"H\", \"C\"]\n", "r_cut = 3.5 #cutoff in Angstrom\n", "n_max = 4 # \n", "sigma = 0.3 #stdev of gaussians\n", "l_max = 4\n", "\n", "# Setting up the SOAP descriptor\n", "soap = SOAP(\n", " species=species,\n", " periodic=False,\n", " r_cut=r_cut,\n", " n_max=n_max,\n", " l_max=l_max,\n", " sigma=sigma,\n", ")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2dc57b73", "metadata": {}, "outputs": [], "source": [ "# Create SOAP output for all frames\n", "soaps = soap.create(traj[:10])\n", "#soaps = soap.create(traj)\n", "print(soaps.shape)" ] }, { "cell_type": "markdown", "id": "23bcc7f5", "metadata": {}, "source": [ "We get a descriptor array with 3 dimensions: no. of frames, no. of atoms, and length of SOAP feature vector. \n", "\n", "This is the input we need to fit machine learning interatomic potentials." ] }, { "cell_type": "code", "execution_count": null, "id": "b7d619f3", "metadata": {}, "outputs": [], "source": [ "plt.plot(soaps[0,0,:]) #plot the soap vector for the first geometry, first atom\n", "plt.title('Full SOAP feature vector for 1st geometry, 1st atom')" ] }, { "cell_type": "markdown", "id": "aec2c1a1", "metadata": {}, "source": [ "**Task for you**\n", "\n", "
\n", " \n", "- Play around with r_cut, n_max, l_max and study how the length of the descriptor changes.\n", "\n", "- Show that the SOAP descriptors are also invariant w.r.t. to translations, permutations, and rotations\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "3c5bdd53", "metadata": {}, "outputs": [], "source": [ "# space to play around" ] }, { "cell_type": "markdown", "id": "b50b00de", "metadata": {}, "source": [ "Measuring the similarity of structures becomes easy when the feature vectors represent the whole structure, such as in the case of Coulomb matrix or MBTR. In these cases the similarity between two structures is directly comparable with different kernels, e.g. the linear or Gaussian kernel.\n", "\n", "Sometimes, we are focussed on comparative analysis between systems, then it is more useful to have a global descriptor. We can transform SOAP (and other atom-centered descriptors) into a global descriptor. We can achieve this by generating an average kernel. DScribe enables this with the option `average='inner'`.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a0a5d025", "metadata": {}, "outputs": [], "source": [ "# Setting up the SOAP descriptor\n", "soap = SOAP(\n", " species=species,\n", " periodic=False,\n", " r_cut=r_cut,\n", " n_max=n_max,\n", " l_max=l_max,\n", " sigma=sigma,\n", " average='inner',\n", ")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "33be2c9f", "metadata": {}, "outputs": [], "source": [ "# Create SOAP output for all frames\n", "soaps = soap.create(traj[:10])\n", "#soaps = soap.create(traj)\n", "print(soaps.shape)\n" ] }, { "cell_type": "markdown", "id": "336dc4ef", "metadata": {}, "source": [ "Now we get only one SOAP descriptor per frame as we have averaged over the atoms" ] }, { "cell_type": "code", "execution_count": null, "id": "bd3216a9", "metadata": {}, "outputs": [], "source": [ "plt.plot(soaps[0,:]) #plot the soap vector for the first geometry, first atom\n", "plt.title('Full SOAP feature vector for 1st geometry, 1st atom')" ] }, { "cell_type": "code", "execution_count": null, "id": "a5d21a82", "metadata": {}, "outputs": [], "source": [ "counts, bins = np.histogram(soaps.flatten(),bins=100)\n", "#plt.stairs(counts, bins)\n", "plt.hist(bins[:-1], bins, weights=counts)\n", "plt.xlabel('SOAPS')\n", "plt.ylabel('count')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "a57c7149", "metadata": {}, "source": [ "**Tasks for you**\n", "
\n", " \n", "- When SOAP is an atom-centered descriptor, 3.5 Angstrom cutoff radius around each atom pretty much covers all relevant information for cyclohexane. When we average and create a global descriptor, is this still the case? \n", "- Try to generate SOAP descriptors for the QM7 dataset. We will need this later.\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "5c59eb7a", "metadata": {}, "outputs": [], "source": [ "#space to play around" ] }, { "cell_type": "markdown", "id": "44d329e0", "metadata": {}, "source": [ "We will now generate several different SOAP descriptors for the cyclohexane dataset (probably a good time to stretch your legs, grab a drink or take a break)" ] }, { "cell_type": "code", "execution_count": null, "id": "6e1547bf", "metadata": {}, "outputs": [], "source": [ "#Let's generate two different soap descriptors with different settings so we hvae them handy\n", "import collections.abc\n", "collections.Iterable = collections.abc.Iterable # this is unfortunately necessary for compatbility of describe with Py>=3.10\n", "from dscribe.descriptors import SOAP\n", "\n", "####DEFAULT SOAP\n", "\n", "species = [\"H\", \"C\"]\n", "r_cut = 3.5 #cutoff in Angstrom\n", "n_max = 4 # \n", "sigma = 0.3 #stdev of gaussians\n", "l_max = 4\n", "\n", "# Setting up the SOAP descriptor\n", "soap = SOAP(\n", " species=species,\n", " periodic=False,\n", " r_cut=r_cut,\n", " n_max=n_max,\n", " l_max=l_max,\n", " sigma=sigma,\n", " average='inner',\n", ")\n", "\n", "soaps_default = soap.create(traj)\n", "\n", "####SOAP with high n and l\n", "\n", "species = [\"H\", \"C\"]\n", "r_cut = 3.5 #cutoff in Angstrom\n", "n_max = 6 # \n", "sigma = 0.3 #stdev of gaussians\n", "l_max = 6\n", "\n", "# Setting up the SOAP descriptor\n", "soap = SOAP(\n", " species=species,\n", " periodic=False,\n", " r_cut=r_cut,\n", " n_max=n_max,\n", " l_max=l_max,\n", " sigma=sigma,\n", " average='inner',\n", ")\n", "\n", "soaps_high = soap.create(traj)\n", "\n", "#### SOAP with even higher n and l and larger cutoff\n", "\n", "#soap_big = SOAP(\n", "# species=species,\n", "# periodic=False,\n", "# r_cut=4.0,\n", "# n_max=8,\n", "# l_max=6,\n", "# sigma=0.3,\n", "# average='inner',\n", "#)\n", "\n", "#soaps_big=soap.create(traj)#, n_jobs=6) # we can parallelize over cores with n_jobs to speed things up" ] }, { "cell_type": "markdown", "id": "1fa8468c", "metadata": {}, "source": [ "**Tasks for you**\n", "
\n", " \n", "- Which of the three descriptors with default parameters is likely to be the best? Rank them by their computational expense for the ML tasks that follow below?\n", "- Before moving on, please go back to CM and MBTR and make sure you reset all values to default and generate descriptors for the full cyclohexane dataset and the QM7 dataset. \n", "\n", "
" ] }, { "cell_type": "markdown", "id": "75e6f4b0", "metadata": {}, "source": [ "## Part 3: Similarity Analysis" ] }, { "cell_type": "markdown", "id": "2b28caf5", "metadata": {}, "source": [ "We can use our new descriptors to find the configurational change of the planar configuration without knowing anything about the structure. In the following box, we use a pairwise kernel from scikit-learn to compare structures for their similarity. For each trajectory, we compare the last frame with all other frames. We have the CM, the MBTR (inverse distance), and the three SOAP vectors to compare initially. Of course, you can play around with settings of all descriptors to find even better representations.\n", "\n", "The ability of a descriptor to measure \"distances\" in the space of the dataset is a crucial prerequisite for being able to build reliable machine learning models." ] }, { "cell_type": "code", "execution_count": null, "id": "3fa3f7e2", "metadata": {}, "outputs": [], "source": [ "#descriptors = soaps_default\n", "#descriptors = soaps_high\n", "#descriptors = coulomb_matrices\n", "descriptors = mbtrs\n", "\n", "#descriptors = soaps_big\n", "\n", "similarities = np.zeros((len(traj), 5))\n", "kernel = partial(pairwise_kernels, metric='rbf', gamma=2.0) #allows us to call pairwise_kernel like a function later\n", "\n", "ranges2 = ranges.copy()\n", "ranges2[4,1] = -1\n", "for i, ci in enumerate(tqdm(conf_idx)):\n", " ki = kernel(descriptors[ranges2[i,0]].reshape(1,-1), descriptors)\n", " similarities[:, i] = ki\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d85595f3", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize=(3*4.8528, 3*1.2219))\n", "\n", "for i, (n, c, r, rgb) in enumerate(zip(names, c_energy, ranges, rgb_colors)):\n", " ax.plot(range(0, r[1] - r[0]),\n", " similarities[r[0]:r[1],i],\n", " label=n,\n", " c=rgb,\n", " zorder=-1)\n", " \n", "ax.legend()\n", "ax.set_xlabel(\"Simulation Timestep\")\n", "ax.set_ylabel(\"Similarity\")\n", "\n", "ax.set_xlim([0, len(energy)//5])\n", "#ax.set_ylim([-0.1, 1.25 * (max_e - min_e)])\n", "#ax.set_yticklabels([])\n", "\n", "plt.tight_layout()\n", "plt.savefig('similarity_soap_big.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4a881d1e", "metadata": {}, "source": [ "A low similarity means that the descriptor considers the structures to be maximally unsimilar. \n", "\n", "- As we can see, the coulomb matrix is not good at identifying similarity. (See what happens with 'sorted_l2' and 'eigenspectrum'. What is the origin of the difference?)\n", "- MBTR and SOAP do a lot better. \n", "- Using the SOAP descriptors, we see that all configurations retain similarity with the initial structure during the dynamics except for the planar structure, which starts in a structure that is disimilar to the final frame." ] }, { "cell_type": "markdown", "id": "592bb009", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- Go back up and retry this with the CM and MBTR descriptors with different settings. Then try tuning the SOAP with higher n_max and l_max. Can you improve the similarity description with MBTR and SOAP? Which descriptor performs best? The similarity difference between the initial and final structure of the planar geometry is a good indicator of the \"resolving power\" of the descriptor. Before we do any machine learning, we should make sure that the resolving power and the structural covariance that our descriptor captures is as large as possible. The descriptor that achieves the biggest disimilarity between different structures should be our descriptor of choice.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "93309b22", "metadata": {}, "outputs": [], "source": [ "#space to play" ] }, { "cell_type": "markdown", "id": "f2f12c40", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- If time permits, try to do the same for the QM7 dataset. Just look at the similarity between two molecules. Which descriptor best captures (dis)similarity between structures?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "64c8cbbc", "metadata": {}, "outputs": [], "source": [ "#space to play" ] }, { "cell_type": "markdown", "id": "1208549d", "metadata": {}, "source": [ "## Part 4: Dimensionality Reduction" ] }, { "cell_type": "markdown", "id": "9dbbbd4e", "metadata": {}, "source": [ "Let's look at dimensionality reduction and clustering techniques. Can we come up with representations that clearly resolve the structural differences in the data?" ] }, { "cell_type": "markdown", "id": "015ebbeb", "metadata": {}, "source": [ "### Principal Component Analysis (PCA)" ] }, { "cell_type": "markdown", "id": "82fb9146", "metadata": {}, "source": [ "PCA is the simplest linear dimensionality reduction technique based on Singular value Decomposition. The idea is to find a small set of principal components that represent the maximum amount of variance (diversity) in the data. We then project the data into that lower dimensional space. There are many variations on this technique including non-linear ones (Kernel PCA)." ] }, { "cell_type": "code", "execution_count": null, "id": "83cdf1d4", "metadata": {}, "outputs": [], "source": [ "from sklearn.decomposition import PCA\n", "\n", "#pick either \n", "descriptors = soaps_default\n", "#descriptors = soaps_high\n", "#descriptors = soaps_big\n", "#descriptors = coulomb_matrices\n", "#descriptors = mbtrs\n", "\n", "n_components = 5\n", "\n", "pca = PCA(n_components=n_components) #looking at the 5 largest eigenvalues\n", "pca.fit(descriptors)\n", "print(pca.explained_variance_ratio_)\n", "print(pca.singular_values_)\n", "# this is the part where we transform the data into the PCA representation. Now each element of the descriptor corresponds to the respective PCA eigenvalue\n", "X_new = pca.transform(descriptors) \n", "\n", "plt.xlabel(\"PCA components\")\n", "plt.ylabel(\"Variance ratio\")\n", "\n", "plt.bar(np.arange(5),pca.explained_variance_ratio_, edgecolor='black')\n", "\n", "print('Total variance of {0} lowest eigenvalues: {1}'.format(n_components, np.sum(pca.explained_variance_ratio_)))\n", "print('Total variance of 2 lowest eigenvalues {0}'.format(np.sum(pca.explained_variance_ratio_[:2])))" ] }, { "cell_type": "markdown", "id": "8929047f", "metadata": {}, "source": [ "The above plot shows us the variance captured by the 5 largest eigenvalues. Together the 5 components dominantly represent the data. The 2 lowest components cover a significant part of the variance." ] }, { "cell_type": "markdown", "id": "9098104f", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- Which descriptor gives the highest variance for the first two PCA components? Can you tune that descriptor to maximize that number?\n", "
" ] }, { "cell_type": "markdown", "id": "c8e693c3", "metadata": {}, "source": [ "We pick the two components with the highest variance over the data and plot them against each other" ] }, { "cell_type": "code", "execution_count": null, "id": "8b0f2db2", "metadata": {}, "outputs": [], "source": [ "print(X_new.shape)\n", "descriptors = X_new[:,:2]\n", "print(descriptors.shape)" ] }, { "cell_type": "code", "execution_count": null, "id": "540e80bc", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "for (r, rgb) in zip(ranges,rgb_colors):\n", " ax.scatter(descriptors[r[0]:r[1],0], descriptors[r[0]:r[1],1], edgecolors='white', color=rgb, linewidths=0.4, label=\"datapoints\", alpha=0.8)\n", "\n", "plt.xlabel(\"PC1\")\n", "plt.ylabel(\"PC2\")\n", "fig.set_figheight(6.0)\n", "fig.set_figwidth(6.0)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "03ed01b9", "metadata": {}, "source": [ "Let's only plot the 'planar' trajectory and colour the dots according to the trajectory time. Red are early times, blue are late times." ] }, { "cell_type": "code", "execution_count": null, "id": "f8fb712f", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "cm = plt.get_cmap('RdYlBu')\n", "colors = range(10000)\n", "\n", "sc = ax.scatter(descriptors[40000:,0], descriptors[40000:,1], edgecolors='white', \n", " c=colors, linewidths=0.4, label=\"datapoints\", alpha=0.8, cmap=cm)\n", "fig.colorbar(sc)\n", "\n", "\n", "plt.xlabel(\"PC1\")\n", "plt.ylabel(\"PC2\")\n", "fig.set_figheight(6.0)\n", "fig.set_figwidth(6.0)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d20fc8a6", "metadata": {}, "source": [ "We can see that the deformation of the molecule out of the planar geometry can be captured well with the two principal components. However, the subtle differences between structures are not well resolved. We will try another method below." ] }, { "cell_type": "markdown", "id": "1a94a6a8", "metadata": {}, "source": [ "**Task for you**\n", "
\n", " \n", "- Compare the PCA plots for CM, MBTR, SOAP(n_max=4, l_max=4) (`soaps_default`), SOAP(n_max=6,l_max=6) (`soaps_high`). \n", "\n", "
" ] }, { "cell_type": "markdown", "id": "7d2f0ca8", "metadata": {}, "source": [ "### t-SNE" ] }, { "cell_type": "markdown", "id": "df4d9cc9", "metadata": {}, "source": [ " t-Distributed Stochastic Neighbor Embedding (t-SNE) is a popular dimensionality-reduction algorithm for visualizing high-dimensional data sets. It can be very informative (and sometimes also a bit misleading). Both methods are nonlinear (in contrast to PCA). t-SNE focuses on preserving the pairwise similarities between data points in a lower-dimensional space. t-SNE is concerned with preserving small pairwise distances whereas, PCA focuses on maintaining large pairwise distances to maximize variance. \n", " \n", " So PCA preserves the variance in the data, whereas t-SNE preserves the relationships between data points in a lower-dimensional space, making it quite a good algorithm for visualizing complex high-dimensional data and relationships between individual datapoints. \n", "\n", " Popular libraries for t-SNE (or an alternative algorithm called UMAP) are\n", " - [openTSNE](https://opentsne.readthedocs.io/en/stable/)\n", " - [umap-learn](https://umap-learn.readthedocs.io/en/latest/)" ] }, { "cell_type": "code", "execution_count": null, "id": "f21f05bc", "metadata": {}, "outputs": [], "source": [ "from openTSNE import TSNE" ] }, { "cell_type": "code", "execution_count": null, "id": "cfb9be97", "metadata": {}, "outputs": [], "source": [ "#list of colors where each point is colored according to its trajectory\n", "\n", "subsampling = 10 #consider only every 10th datapoint\n", "\n", "data_length = 10000//subsampling\n", "\n", "colors = []\n", "for c in rgb_colors:\n", " for i in range(data_length):\n", " colors.append(c)\n", "\n", "#list of colors for one trajectory which indicates time\n", "\n", "\n", "##pick a descriptor\n", "descriptors = soaps_default[::subsampling,:]\n", "#descriptors = soaps_high[::subsampling,:]\n", "#descriptors = soaps_big[::subsampling,:]\n", "#descriptors = coulomb_matrices[::subsampling,:] \n", "#descriptors = mbtrs[::subsampling,:]" ] }, { "cell_type": "markdown", "id": "8bd69e06", "metadata": {}, "source": [ "The **perplexity** is an important tunable parameter for t-SNE. \n", "Here we can see how increasing the perplexity (number of expected neighbors) changes the layout of the projection." ] }, { "cell_type": "markdown", "id": "485b323b", "metadata": {}, "source": [ "
\n", " \n", "This algorithm may take several minutes. Be sure to parallelize it (`n_jobs`) and be patient. Above, we have already sliced the datasets to only take every 10th datapoint (1000 data points for each trajectory).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "060319cf", "metadata": {}, "outputs": [], "source": [ "perplexities = np.logspace(0, 2.69, 6, dtype=int)\n", "fig, ax = plt.subplots(1,\n", " len(perplexities),\n", " figsize=(4 * len(perplexities), 4),\n", " )\n", "\n", "for i, perp in enumerate(tqdm(perplexities)):\n", " tsne = TSNE(\n", " n_components=2, # number of components to project across\n", " perplexity= perp, \n", " metric=\"euclidean\", # distance metric\n", " n_jobs=6, # parallelization\n", " random_state=42,\n", " verbose=False,\n", " )\n", " t_tsne = tsne.fit(descriptors)\n", " ax[i].scatter(*t_tsne.T, c=colors, s=2)\n", " ax[i].axis('off')\n", " ax[i].set_title(\"Perplexity = {}\".format(perp))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4a9079b6", "metadata": {}, "source": [ "For MBTR and SOAP and for large enough perplexity values, we find that the data clusters into 2 very large clusters and 1 very little one. The colours we are using are still the colours associated with the five cyclohexane MD runs. The t-SNE is proof that the 5 conformers start in different places but equilibrate into only two structures, namely the `boat` (green, red, yellow trajectories) and the `chair` conformation (blue and purple trajectories). \n", "\n" ] }, { "cell_type": "markdown", "id": "aab208b2", "metadata": {}, "source": [ "Now let's take the most suitable perplexity value and look at the t-SNE, but now we color it according to the trajectory time. Red points are early, purple points are towards the end of the trajectory. We represent the different trajectory types with different markers. FOcus on the circles (planar trajectory)." ] }, { "cell_type": "code", "execution_count": null, "id": "0f7fcf07", "metadata": {}, "outputs": [], "source": [ "perplexity = 50\n", "plt.figure(figsize=(4,4),dpi=300)\n", "\n", "cm = plt.cm.get_cmap('RdYlBu')\n", "colors_time = []\n", "for r in range(5):\n", " for i in np.linspace(0,1,data_length):\n", " colors_time.append(cm(i))\n", "\n", "markers = []\n", "for c in ['v','^','<','>','o']:\n", " for i in range(data_length):\n", " markers.append(c)\n", "\n", "tsne = TSNE(\n", " n_components=2, # number of components to project across\n", " perplexity=perplexity, # amount of neighbors one point is posited to have... play around with this!\n", " metric=\"euclidean\", # distance metric\n", " n_jobs=4, # parallelization\n", " random_state=42,\n", " verbose=False,\n", ")\n", "\n", "size = 4\n", "t_tsne = tsne.fit(descriptors)\n", "for i,m in enumerate(['v','^','<','>','o']):\n", " r = ranges[i]//subsampling\n", " tmp = t_tsne[r[0]:r[1]]\n", " plt.scatter(*tmp.T, c=colors_time[r[0]:r[1]], marker=m, s=size)\n", "\n", "plt.axis('off')\n", "plt.title(\"Perplexity = {}\".format(perplexity))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4a4ce34b", "metadata": {}, "source": [ "The planar trajectory (circles) starts off in a place that is structurally similar to the initial configurations of the `boat` (yellow previousy, '>' signs here) trajectory. While the `boat` trajectory converges into the same space as the half-chair and twist-boat, the `planar` trajectory moves through the `half-chair` configuration and the `boat` configuration before it ends up in the `chair` configuration." ] }, { "cell_type": "markdown", "id": "3ab8bbb0", "metadata": {}, "source": [ "**Task for you**:\n", "
\n", " \n", "- A big part of ML and data science is data visualization. Can you come up with a better way to visualize the dynamical changes of conformations along the trajectories based on what we know now?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "62efa79c", "metadata": {}, "outputs": [], "source": [ "#space to play around" ] }, { "cell_type": "markdown", "id": "c9ad45ce", "metadata": {}, "source": [ "**Task for you**:\n", "
\n", " \n", "- Surely, now you might be interested to make PCA plots of QM7. Can you plot the molecules into a 2D PCA descriptor space?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "630e518f", "metadata": {}, "outputs": [], "source": [ "#space to play around\n" ] }, { "cell_type": "markdown", "id": "6420564c", "metadata": {}, "source": [ "## Part 5: Clustering" ] }, { "cell_type": "markdown", "id": "c7965469", "metadata": {}, "source": [ "Let's assume we use PCA or t-SNE to perform dimensionality reduction and in the process, we have found a way to visualize the landscape of the data and relevant similarities. If we trust the space, we might want to perform clustering to identify the $n$ most structurally diverse geometries from this space. \n", "\n", "This would allow us for example to pick a set of structures that best represent the diversity of this space to perform subsequent electronic structure calculations.\n", "\n", "We will use [k-means clustering](https://scikit-learn.org/stable/modules/clustering.html#k-means) based on the previous PCA to achieve this. We will use the computationally more efficient batched variant of KMeans." ] }, { "cell_type": "code", "execution_count": null, "id": "83b3afba", "metadata": {}, "outputs": [], "source": [ "from sklearn.decomposition import PCA\n", "\n", "#pick either \n", "descriptors = soaps_default\n", "descriptors = soaps_high\n", "#descriptors = soaps_big\n", "#descriptors = coulomb_matrices\n", "#descriptors = mbtrs\n", "\n", "###regenerate PCA\n", "n_components = 5\n", "pca = PCA(n_components=n_components)\n", "pca.fit(descriptors)\n", "X_new = pca.transform(descriptors) \n", "print(X_new.shape)\n", "descriptors = X_new[:,:2]" ] }, { "cell_type": "code", "execution_count": null, "id": "4acff3cd", "metadata": {}, "outputs": [], "source": [ "from sklearn.cluster import MiniBatchKMeans\n", "\n", "niter = 1000000\n", "minibatch = 1800 # Batch size (how many data points are treated in the algorithm at once, reduce if you run out of memory)\n", "ncluster= 20 # number of clusters\n", "\n", "kmeans = MiniBatchKMeans(n_clusters=ncluster,\n", " init=\"k-means++\",\n", " max_iter=niter,\n", " n_init=3,\n", " batch_size=minibatch)\n", "# Fit the clustering model\n", "km = kmeans.fit(np.array(descriptors).reshape(-1,1))\n", "# Use the clustering model to predict the cluster number of each data point\n", "indices = km.fit_predict(descriptors)" ] }, { "cell_type": "markdown", "id": "4ca23e20", "metadata": {}, "source": [ "Once we have performed the clustering, we can calculate the centres of clusters, that is, the indices of the datapoints that are closest to the centroids of the clusters." ] }, { "cell_type": "code", "execution_count": null, "id": "50773396", "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics.pairwise import pairwise_distances_argmin\n", "\n", "centroid = kmeans.cluster_centers_\n", "b = np.inf\n", "ind = pairwise_distances_argmin(centroid, descriptors) # measure which structure is closest to the centroids of the clusters\n", "print(ind)\n", "\n", "X_centers = []\n", "for i in ind:\n", " X_centers.append(descriptors[i])\n", "X_centers = np.array(X_centers)" ] }, { "cell_type": "markdown", "id": "79a48fb5", "metadata": {}, "source": [ "Now you could save the indices and you can identify which conformations belong to thse clusters. Let's look at them" ] }, { "cell_type": "code", "execution_count": null, "id": "3160a14a", "metadata": {}, "outputs": [], "source": [ "mol_traj = []\n", "for i in ind:\n", " mol_traj.append(traj[i])\n", "\n", "# visualize the trajectory of the clustered molecules\n", "viewer = WeasWidget()\n", "viewer.from_ase(mol_traj)\n", "viewer.avr.model_style = 1\n", "viewer.avr.show_hydrogen_bonds = True\n", "viewer" ] }, { "cell_type": "markdown", "id": "3bd03d31", "metadata": {}, "source": [ "Now we replot the PCA plot while coloring the different clusters and highlighting their centres." ] }, { "cell_type": "code", "execution_count": null, "id": "9f123ff1", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "ax.scatter(X_new[:,0], X_new[:,1],c=indices, edgecolors='white', linewidths=0.4, label=\"datapoints\", alpha=0.8)\n", "ax.scatter(X_centers[:,0],X_centers[:,1],color=\"white\",label=\"centers\",edgecolors='black',linewidths=0.7,marker=\"X\", alpha=0.8)\n", "\n", "plt.legend(fancybox=True,framealpha=1,edgecolor='black',handletextpad=0.05,borderpad=0.3,handlelength=1.2,columnspacing=0.4,labelspacing=0.2,ncol=1,loc=1, fontsize=\"medium\") #, bbox_to_anchor=(1.75, 1.02))\n", "plt.xlabel(\"PC1\")\n", "plt.ylabel(\"PC2\")\n", "fig.set_figheight(6.0)\n", "fig.set_figwidth(6.0)\n", "\n", "plt.show()\n", "# Two distinct parts of PC1/2 space, so pick enough clusters to avoid clusters spanning both parts. \n" ] }, { "cell_type": "markdown", "id": "299911a3", "metadata": {}, "source": [ "We can also see if our clusters cover a wide range of values for PC1" ] }, { "cell_type": "code", "execution_count": null, "id": "7ff2a8db", "metadata": {}, "outputs": [], "source": [ "from matplotlib.ticker import MaxNLocator, FuncFormatter\n", "\n", "ax = plt.axes()\n", "\n", "ax.scatter(X_new[:,0], indices, edgecolors='white', linewidths=0.4, alpha=0.8)\n", "\n", "plt.xlabel(\"PC1\")\n", "plt.ylabel(\"Cluster\")\n", "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n", "ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b8d1952b", "metadata": {}, "source": [ "**Task for you**:\n", "
\n", " \n", "- Try playing around with the number of clusters. How big a number would still be a sensible set of clusters?\n", "
" ] }, { "cell_type": "markdown", "id": "460396e9", "metadata": {}, "source": [ "
\n", " \n", "This is a very simplified example of clustering and if the PCA and the descriptors are not resolving the diversity of the space really well, you are running the risk of selecting \"diverse\" structures with the wrong measure. There are many more advanced methods for clustering, which you should explore. I recommend the resource on Unsupervised Learning quoted at the top of the notebook and the links at the end of the notebook.\n", "
" ] }, { "cell_type": "markdown", "id": "2fbc0fc3", "metadata": {}, "source": [ "## Part 6: Supervised Learning / Model Fitting" ] }, { "cell_type": "markdown", "id": "9f9ab496", "metadata": {}, "source": [ "In a later tutorial, you will learn all about neural networks and deep learning models for interatomic potentials. Here, we will look at models based on linear regression and Kernel Ridge Regression as they form a good starting point. Our aim is to create a model that can predict the energy of any cyclohexane conformation. The previously calculated feature vectors (`descriptors`) are our input quantities ($\\mathbf{x}$), the energies are our labels.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8e1becc0", "metadata": {}, "outputs": [], "source": [ "#pick either \n", "#descriptors = soaps_default\n", "descriptors = soaps_high\n", "#descriptors = soaps_big\n", "#descriptors = coulomb_matrices\n", "#descriptors = mbtrs" ] }, { "cell_type": "markdown", "id": "df1b0b1f", "metadata": {}, "source": [ "### Data Preparation" ] }, { "cell_type": "markdown", "id": "a68178de", "metadata": {}, "source": [ "First we have to split our data into train and test data. The train data is what we will use to train the model. The test data will be \"unseen\" and used to evaluate the prediction capabilities of the model.\n", "\n", "We won't be training on the full training data set, but only on a small subset. We will take 1000 data points.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2fb889d2", "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import train_test_split\n", "from sklearn.utils import resample\n", "\n", "#we know the `energy`` of each frame, the energy of the starting conformers (c_energy), max_e, and min_e\n", "\n", "X = descriptors\n", "#our labels are the energies\n", "y = energy\n", "\n", "X_train_full, X_test, y_train_full, y_test = train_test_split(\n", " X, y, #inputs and outputs\n", " test_size=0.20, #20% of the data is test data\n", " random_state=42, #a random seed\n", " ) \n", "\n", "#only train on a small subset of samples rather than the full training dataset\n", "X_train, y_train = resample(X_train_full, y_train_full, replace=False, n_samples=1000, random_state=42)\n" ] }, { "cell_type": "markdown", "id": "f734006d", "metadata": {}, "source": [ "Next, we standardize the training and testing sets. This means scaling and normalizing the data. This way, the model mean is centred around zero and the model does not have to cover large value ranges. We do this with the inputs and the outputs.\n", "\n", "For the labels/outputs, this is as simple as\n", "$$\n", "\\tilde{\\mathbf{y}} = \\frac{\\mathbf{y} - u}{s}\n", "$$\n", "where u and s are the mean and the standard deviation of the labels." ] }, { "cell_type": "code", "execution_count": null, "id": "47be463a", "metadata": {}, "outputs": [], "source": [ "from sklearn import preprocessing\n", "\n", "#scikit-learn has a utility to help with the feature vectors\n", "\n", "scaler = preprocessing.StandardScaler()\n", "scaler.fit(X_train)\n", "X_train = scaler.transform(X_train)\n", "X_test = scaler.transform(X_test) #? This needs to be the same scaler as for the training data, otherwise the test data would be scaled into the wrong space\n", "\n", "#rescaling the labels\n", "\n", "u = np.mean(y_train)\n", "s = np.std(y_train)\n", "\n", "y_train_scaled = (y_train-u)/s\n", "y_test_scaled = (y_test-u)/s\n", "\n", "def unscale(y):\n", " # this is for transforming the predictions back into proper energies\n", " return (y*s)+u\n" ] }, { "cell_type": "markdown", "id": "cc47b74e", "metadata": {}, "source": [ "### Multiple Linear Regression" ] }, { "cell_type": "markdown", "id": "a8eb7e5b", "metadata": {}, "source": [ "In multiple linear regression, we try to find a set of regression coefficients $\\mathbf{\\beta}$ to approximate a $p$-dimensional function\n", "\n", "$$\n", "f(\\mathbf{x};\\mathbf{\\beta}) = \\sum_{j=1}^{p} \\beta_j x_j = \\mathbf{x}^T \\mathbf{\\beta}\n", "$$" ] }, { "cell_type": "markdown", "id": "f4b5f805", "metadata": {}, "source": [ "In our case, we would like to predict the energy for an arbitrary conformation of cyclohexane. We have already acknowledged that this would not be possible based simply on the Cartesian positions of atoms, which is why we have created various representations that provide us with high dimensional feature vectors (e.g. SOAP).\n", "\n", "The feature vectors form the inputs $\\mathbf{x}$ and the dimensionality of the space ($p$) in which we want to fit a linear relationship is the length of the descriptor. The quality of the fit will depend on the ability of the high dimensional space spanned by the descriptor to express complex relationships as linear trends. The more chemically meaningful (and higher dimensional) the descriptor, the better the fit. \n", "\n", "Let's test this for the energy prediction.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7e56ea31", "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression\n", "\n", "multiple_linear_regression = LinearRegression()\n", "multiple_linear_regression.fit(X = X_train, y = y_train_scaled)\n", "# This is the training part" ] }, { "cell_type": "markdown", "id": "30f112ef", "metadata": {}, "source": [ "Determine $R^2$, MAE, and RMSD of the fit\n", "\n", "$R^2$ metric (Coefficient of determination)\n", "$$\n", "R^2 = 1 - \\frac{\\sum_i(y_i-f_i)^2}{\\sum_i(y_i-u)^2}\n", "$$\n", "where $u$ is the mean of the observed data, $y_i$ is the data value (\"ground truth\") and $f_i$ is the prediction.\n", "\n", "Mean Absolute Error (MAE)\n", "$$\n", "\\mathrm{MAE} = \\frac{\\sum_i |y_i - f_i|}{n}\n", "$$\n", "\n", "Root Mean Square Error (RMSE)\n", "$$\n", "\\mathrm{RMSE} = \\sqrt{\\frac{\\sum_{i} (y_i - f_i)^2}{n} }\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "cba8cf15", "metadata": {}, "outputs": [], "source": [ "# This is the prediction part\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\n", "\n", "pred_train_lr = unscale(multiple_linear_regression.predict(X_train))\n", "pred_test_lr = unscale(multiple_linear_regression.predict(X_test))\n", "\n", "r2 = r2_score(y_test, pred_test_lr)\n", "print('R2 train = ', r2_score(y_train, pred_train_lr))\n", "print('R2 test = ', r2_score(y_test, pred_test_lr))\n", "\n", "mae = mean_absolute_error(y_test, pred_test_lr)\n", "print('MAE test = ', mae)\n", "rmse = np.sqrt(mean_squared_error(y_test,pred_test_lr ))\n", "print('RSME test = ', rmse)" ] }, { "cell_type": "code", "execution_count": null, "id": "730e0fd8", "metadata": {}, "outputs": [], "source": [ "plt.scatter(y_train,pred_train_lr, label='Training Set')\n", "plt.scatter(y_test,pred_test_lr, label='Test Set')\n", "\n", "plt.xlabel('Real')\n", "plt.ylabel('Predicted')\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "76562e5b", "metadata": {}, "source": [ "**Task for you**\n", "\n", "
\n", " \n", "- Which of the descriptors that we have prepared gives the best correlation and test error with a simple linear fit?\n", "- What value of MAE or RMSE on the test set should we consider to be a low error? What are the energy differences between conformers that we need to be able to resolve?\n", "
" ] }, { "cell_type": "markdown", "id": "b0f3e410", "metadata": {}, "source": [ "### Kernel Ridge Regression (KRR) and Gaussian Process Regression (GPR)" ] }, { "cell_type": "markdown", "id": "6072b170", "metadata": {}, "source": [ "**Kernel Ridge Regression** provides a rigorous mathematical framework to extend linear statistical learning techniques to model nonlinear dependencies in high-dimensional data. The method is based on measuring similarity betwen pairs of data points and then using this similarity to unfold nonlinearities in the original data space. The similarity transformation enables the mapping of the data into a higher dimensional space spanned by the training data where the data might be representable in terms of linear functions.\n", "\n", "If we start from linear regression\n", "$$\n", "f(\\mathbf{x};\\mathbf{\\beta}) = \\sum_{j=1}^{p} \\beta_j x_j = \\mathbf{x}^T \\mathbf{\\beta}\n", "$$\n", "we can expand the regression coefficients $\\beta$ via different regression coefficients $\\alpha$ in a much higher dimensional space spanned by the size of the training data ($N_{tr}$). \n", "$$\n", "\\beta_j = \\sum_{i=1}^{N_{\\mathrm{tr}}} \\alpha_i x_{ij}\n", "$$\n", "$x_{ij}$ are elements of the input vectors (feature vectors based on our descriptors). By inserting into the above we get\n", "$$\n", "f(\\mathbf{x}') = \\sum_{i=1}^{N_{tr}} \\alpha_i \\braket{\\mathbf{x}_i,\\mathbf{x}'}\n", "$$\n", "where $\\braket{\\mathbf{x}_i,\\mathbf{x}'}$ is a dot-product of two feature vectors, one associated with data point $i$ and the other one associated with the point for which we would like to make a prediction ($\\mathbf{x}'$). The same can be written in some choice of nonlinearly transformed feature space:\n", "$$\n", "f(\\phi(\\mathbf{x}')) = \\sum_{i=1}^{N_{tr}} \\alpha_i k(\\mathbf{x}_i,\\mathbf{x}')\n", "$$\n", "where $k(\\mathbf{x}_i,\\mathbf{x}')$ is the dot product in feature space (the so-called kernel). The kernel is the similarity measure with which we compare training datapoint $\\mathbf{x}_i$ and the point to be predicted $\\mathbf{x}'$.\n", "\n", "If $k(\\mathbf{x}_i,\\mathbf{x}') = \\braket{\\mathbf{x}_i,\\mathbf{x}'}$, it is called a linear kernel or `dot-product kernel`. This would be equivalent to linear regression.\n", "\n", "We can also choose a Gaussian kernel:\n", "$$ \n", "k(\\mathbf{x}_i,\\mathbf{x}') = \\exp \\left(-\\frac{1}{2\\sigma^2}\\sum_j^{p} (x_{ij} - x_j^{'})^2 \\right) \n", "$$\n", "This is equivalent to mapping the input vector to a space spanned by Gaussian functions (with width $\\sigma$) centred at the training data points.\n", "\n", "The term **Ridge Regression** comes from the fact that a simple inversion of the equation $\\mathbf{y}=\\mathbf{K}\\mathbf{\\alpha}$ to determine $\\mathbf{\\alpha}$ is numerically unstable. This is often connected with individual parameters becoming extremely large. In Ridge regression, we use a regularized form of linear regression, where we minimize the L2 loss between prediction and training data under the constraint that the coefficients remain small (Tikhonov regularization). That means, instead of solving\n", "$$ \\mathbf{\\alpha} = \\mathbf{K}^{-1}\\mathbf{y} $$\n", "we solve\n", "$$ \n", "\\mathbf{\\alpha} = (\\mathbf{K}^{-1}+\\lambda\\mathbf{I})\\mathbf{y}\n", "$$\n", "where $\\mathbf{I}$ is the identity matrix and $\\lambda$ is a small nonnegative smoothing parameter (regularization parameter). Regularization is crucial to reduce the risk of overfitting.\n", "\n", "\n", "Gaussian kernel fits are defined by two hyperparameters, $\\sigma$ and $\\lambda$. If $\\sigma$ is too small (Gaussians might not overlap), the model might not generalize well, but if $\\sigma$ is too large, the model will fail to learn from the data. \n", "\n", "
\n", " \n", "Typically, hyperparameters within the ML model need to be optimised to achieve optimal accuracy and generalisability.\n", "
\n", "\n", "**Gaussian Process Regresssion** is related to KRR as it generates fits to the same mathematical form of the regression function, but it does so with a Bayesian approach that automatically gives the model variance as a byproduct. Another byproduct is that it often can automatically find optimal hyperparameters (via the log-marginal likelihood).\n" ] }, { "cell_type": "markdown", "id": "07221c93", "metadata": {}, "source": [ "Let's start by reproducing the linear fit with a linear kernel" ] }, { "cell_type": "code", "execution_count": null, "id": "3e5343db", "metadata": {}, "outputs": [], "source": [ "from sklearn.kernel_ridge import KernelRidge\n", "\n", "krr = KernelRidge(alpha=0.00001, kernel='linear') # alpha here is lambda in the explanation above\n", "krr.fit(X_train, y_train_scaled)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9d6f97c7", "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\n", "\n", "pred_train_kr = unscale(krr.predict(X_train))\n", "pred_test_kr = unscale(krr.predict(X_test))\n", "\n", "r2 = r2_score(y_test, pred_test_kr)\n", "print('R2 train = ', r2_score(y_train, pred_train_kr))\n", "print('R2 test = ', r2_score(y_test, pred_test_kr))\n", "\n", "mae = mean_absolute_error(y_test, pred_test_kr)\n", "print('MAE test = ', mae)\n", "rmse = np.sqrt(mean_squared_error(y_test,pred_test_kr ))\n", "print('RSME test = ', rmse)\n", "\n", "plt.scatter(y_train,pred_train_kr, label='Training Set')\n", "plt.scatter(y_test,pred_test_kr, label='Test Set')\n", "\n", "plt.xlabel('Real')\n", "plt.ylabel('Predicted')\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "71448fdb", "metadata": {}, "source": [ "Next we try the Gaussian (RBF) kernel. Here we have to set the width of the RBFs" ] }, { "cell_type": "code", "execution_count": null, "id": "a51ac10f", "metadata": {}, "outputs": [], "source": [ "from sklearn.kernel_ridge import KernelRidge\n", "\n", "krr = KernelRidge(alpha=0.00001, kernel='rbf', gamma=0.005)\n", "krr.fit(X_train, y_train_scaled)" ] }, { "cell_type": "code", "execution_count": null, "id": "9ee4e00d", "metadata": {}, "outputs": [], "source": [ "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\n", "\n", "pred_train_kr = unscale(krr.predict(X_train))\n", "pred_test_kr = unscale(krr.predict(X_test))\n", "\n", "r2 = r2_score(y_test, pred_test_kr)\n", "print('R2 train = ', r2_score(y_train, pred_train_kr))\n", "print('R2 test = ', r2_score(y_test, pred_test_kr))\n", "\n", "mae = mean_absolute_error(y_test, pred_test_kr)\n", "print('MAE test = ', mae)\n", "rmse = np.sqrt(mean_squared_error(y_test,pred_test_kr ))\n", "print('RSME test = ', rmse)\n", "\n", "plt.scatter(y_train,pred_train_kr, label='Training Set')\n", "plt.scatter(y_test,pred_test_kr, label='Test Set')\n", "\n", "plt.xlabel('Real')\n", "plt.ylabel('Predicted')\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "8fc19456", "metadata": {}, "source": [ "Note how the model fits the trained data almost perfectly and much better than the linear model, but does not capture the test data quite as well. This means we are overfitting. This can have several reasons:\n", "\n", "- The data is simple and the model is too complex\n", "- The feature space is too large\n", "- Too much noise and outliers in training data (can be excluded here)\n", "- Not enough data\n", "- The model is not sufficiently constrained and the hyperparameters are not optimal.\n", "\n", "We need to find a way to properly evaluate and fine-tune the model." ] }, { "cell_type": "markdown", "id": "99fb715b", "metadata": {}, "source": [ "
\n", "\n", "In this case, the linear kernel is the right choice because the data is simple enough. There is no need to go to a Gaussian kernel, but we'll still try to optimise what we can get out of the Gaussian kernel (for educational purposes).\n", "
" ] }, { "cell_type": "markdown", "id": "c00326e5", "metadata": {}, "source": [ "### Cross Validation" ] }, { "cell_type": "markdown", "id": "002f1678", "metadata": {}, "source": [ "Cross validation (CV) is a resampling method that uses different portions of data to test and train models. It is a way to assess how well a model generalises and to robustly evaluate model performance. We can also use cross validation to check for outliers in the dataset. As shown below, the data is split into training and test data. The test data is then further split into $k$ \"folds\" (here 5). Over $k$ possible permutations, we can change which part of the training data we \"hold out\". That means we can train $k$ different models on different parts of the data. \n", "If the standard deviation over several folds is low, the model is robust. If the MAE on the hold-out set is very similar, it likely also has few outliers. The standard deviation can also provide us with a measure of the uncertainty of the model, which could be used in the context of committee-based active learning. \n", "\n", "CV is also used when the space of hyperparameters is explored to find the optimal model layout. For KRR, this is $\\sigma$ and $\\lambda$, for neural networks, it can be parameters that define the layout of the network. " ] }, { "cell_type": "markdown", "id": "2ccb8dfb", "metadata": {}, "source": [ "\n", "![crossval](https://scikit-learn.org/stable/_images/grid_search_cross_validation.png\n", ")" ] }, { "cell_type": "markdown", "id": "f588ab75", "metadata": {}, "source": [ "Let's evaluate the mean prediction error and standard deviation for the linear model and the KRR model with given settings." ] }, { "cell_type": "code", "execution_count": null, "id": "e78a05df", "metadata": {}, "outputs": [], "source": [ "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.model_selection import cross_validate\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer\n", "\n", "krr = KernelRidge(alpha=0.00001, kernel='linear')\n", "#krr.fit(X_train, y_train_scaled) #? This would be fitting twice - Check documentation for the cross validation tool in question\n", "\n", "#this is a dictionary which we can use to pass multiple scoring functions into CV (even custom-made ones)\n", "scoring = {\n", " 'r2_score': make_scorer(r2_score),\n", " 'MAE': make_scorer(mean_absolute_error),\n", " 'RMSE': make_scorer(mean_squared_error),\n", " }\n", "\n", "scores = cross_validate(\n", " krr, #model\n", " X_train, #training data\n", " y_train_scaled, #training data\n", " scoring=scoring, # types of model scores\n", " cv=5, # k-foldedness of cross validation\n", " return_estimator=False, #this function can return all five models if needed\n", " return_train_score=True\n", " )\n", "\n", "print('Linear Model')\n", "print('Train Errors:')\n", "print(' The mean and stdev of train MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_MAE']*s), np.std(np.sqrt(scores['train_MAE'])*s)))\n", "print(' The mean and stdev of train RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_RMSE']*s), np.std(np.sqrt(scores['train_RMSE'])*s)))\n", "print('Test Errors:')\n", "print(' The mean and stdev of test MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_MAE']*s), np.std(np.sqrt(scores['test_MAE'])*s)))\n", "print(' The mean and stdev of test RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_RMSE']*s), np.std(np.sqrt(scores['test_RMSE'])*s)))" ] }, { "cell_type": "code", "execution_count": null, "id": "d21f5382", "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_validate\n", "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer\n", "\n", "krr = KernelRidge(alpha=0.00001, kernel='rbf', gamma=0.005)\n", "\n", "#this is a dictionary which we can use to pass multiple scoring functions into CV (even custom-made ones)\n", "scoring = {\n", " 'r2_score': make_scorer(r2_score),\n", " 'MAE': make_scorer(mean_absolute_error),\n", " 'RMSE': make_scorer(mean_squared_error),\n", " }\n", "\n", "scores = cross_validate(\n", " krr, \n", " X_train, \n", " y_train_scaled, \n", " scoring=scoring, \n", " cv=5, \n", " return_estimator=False, #this function can return all five models if needed\n", " return_train_score=True\n", " )\n", "\n", "print('RBF Model')\n", "print('Train Errors:')\n", "print(' The mean and stdev of train MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_MAE']*s), np.std(np.sqrt(scores['train_MAE'])*s)))\n", "print(' The mean and stdev of train RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_RMSE']*s), np.std(np.sqrt(scores['train_RMSE'])*s)))\n", "print('Validation Errors:')\n", "print(' The mean and stdev of test MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_MAE']*s), np.std(np.sqrt(scores['test_MAE'])*s)))\n", "print(' The mean and stdev of test RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_RMSE']*s), np.std(np.sqrt(scores['test_RMSE'])*s)))\n" ] }, { "cell_type": "markdown", "id": "63a7a8a2", "metadata": {}, "source": [ "**Task for you**\n", "\n", "
\n", "\n", "- What have we learned and what do we have to do next? Which one of the listed causes of overfitting might be the case here?\n", "
" ] }, { "cell_type": "markdown", "id": "3cb1c6db", "metadata": {}, "source": [ "### Hyperparameter Optimisation" ] }, { "cell_type": "markdown", "id": "162c9f46", "metadata": {}, "source": [ "We can use the CV to perform a systematic hyperparameter optimisation. In the following we use a scikit-learn feature to perform a hyperparameter grid search." ] }, { "cell_type": "code", "execution_count": null, "id": "d1c99520", "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_validate, GridSearchCV\n", "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer\n", "import pandas as pd\n", "\n", "j = 0.000001\n", "alphas = [ j*(10**i) for i in range(8)]\n", "gammas = [ j*(10**(i/8)) for i in range(64)]\n", "\n", "#vary alphas\n", "#parameters = {'kernel':['rbf'], 'alpha':alphas, 'gamma':[0.005]}\n", "#vary gammas\n", "parameters = {'kernel':['rbf'], 'alpha':[0.0001], 'gamma':gammas}\n", "\n", "scoring = {\n", " 'r2_score': make_scorer(r2_score),\n", " 'MAE': make_scorer(mean_absolute_error),\n", " 'RMSE': make_scorer(mean_squared_error),\n", " }\n", "\n", "krr = KernelRidge()\n", "gridsearch = GridSearchCV(krr, \n", "param_grid=parameters, \n", "scoring=scoring, \n", "cv=5, \n", "n_jobs=6,\n", "refit=False, \n", "return_train_score=True\n", ")\n", "\n", "gridsearch.fit(X_train, y_train_scaled)\n", "\n", "cvdata = pd.DataFrame(gridsearch.cv_results_)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9bb5d1dc", "metadata": {}, "outputs": [], "source": [ "#try also printing the full dataframe by uncommenting the next line and commenting the other one out\n", "#cvdata\n", "cvdata[['param_alpha','param_gamma','param_kernel','mean_test_MAE', 'std_test_MAE', 'mean_train_MAE', 'std_train_MAE']]" ] }, { "cell_type": "code", "execution_count": null, "id": "5d415eec", "metadata": {}, "outputs": [], "source": [ "cvdata['mean_test_MAE'] = cvdata['mean_test_MAE'].multiply(s) #unscale\n", "cvdata['mean_train_MAE'] = cvdata['mean_train_MAE'].multiply(s) #unscale\n", "#plotframe = cvdata[['param_alpha','mean_test_MAE','mean_train_MAE']]\n", "#plotframe.plot('param_alpha')\n", "plotframe = cvdata[['param_gamma','mean_test_MAE','mean_train_MAE']]\n", "plotframe.plot('param_gamma',logx=True)\n", "\n", "\n", "#? \"train\" should be \"validation\" again" ] }, { "cell_type": "markdown", "id": "8724fedc", "metadata": {}, "source": [ "We can exactly see at which RBF width value the train and test errors start to deviate. This is the point above which overfitting occurs.\n", "An alpha value of 0.0001 and a gamma value of 0.001 (RBF width) seem to provide improved results. Let's use these to reevaluate the model and compare with the previous results. " ] }, { "cell_type": "code", "execution_count": null, "id": "02422542", "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import cross_validate\n", "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer\n", "\n", "krr = KernelRidge(alpha=0.0001, kernel='rbf', gamma=0.001)\n", "\n", "#this is a dictionary which we can use to pass multiple scoring functions into CV (even custom-made ones)\n", "scoring = {\n", " 'r2_score': make_scorer(r2_score),\n", " 'MAE': make_scorer(mean_absolute_error),\n", " 'RMSE': make_scorer(mean_squared_error),\n", " }\n", "\n", "scores = cross_validate(\n", " krr, \n", " X_train, \n", " y_train_scaled, \n", " scoring=scoring, \n", " cv=5, \n", " return_estimator=False, #this function can return all five models if needed\n", " return_train_score=True\n", " )\n", "\n", "print('Train Errors:')\n", "print(' The mean and stdev of train MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_MAE']*s), np.std(np.sqrt(scores['train_MAE'])*s)))\n", "print(' The mean and stdev of train RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['train_RMSE']*s), np.std(np.sqrt(scores['train_RMSE'])*s)))\n", "print('Test Errors:')\n", "print(' The mean and stdev of test MAE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_MAE']*s), np.std(np.sqrt(scores['test_MAE'])*s)))\n", "print(' The mean and stdev of test RMSE over the folds are {0:8.5f} and {1:8.5f} eV'.format(np.mean(scores['test_RMSE']*s), np.std(np.sqrt(scores['test_RMSE'])*s)))" ] }, { "cell_type": "markdown", "id": "96d0ae9a", "metadata": {}, "source": [ "**Task for you**\n", "\n", "
\n", "\n", "- Use this grid search code as a starting point to identify the optimal parameters that will provide the lowest test error and will not lead to overfitting of the training data\n", "
" ] }, { "cell_type": "markdown", "id": "8bca89b4", "metadata": {}, "source": [ "### Learning Curves\n", "\n", "Maybe we need more training data. Studying how quickly the prediction error goes down as we add training data means assessing the rate of learning of the model. Learning curves are important to undertsand the ability of models to learn from little data. We can use them to also compare across the four descriptors we have used.\n", "\n", "Generating the learning curve data will take a little while." ] }, { "cell_type": "markdown", "id": "f6dfdff5", "metadata": {}, "source": [ "**Task for you**\n", "\n", "
\n", "\n", "- Study the plot below. Which descriptor achieves the highest error for the smallest training data set?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "bbbd4d8a", "metadata": {}, "outputs": [], "source": [ "#only train on a small subset of samples rather than the full training dataset\n", "from sklearn import preprocessing\n", "from sklearn.utils import resample\n", "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\n", "from sklearn.model_selection import cross_validate\n", "from sklearn.kernel_ridge import KernelRidge\n", "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, make_scorer\n", "\n", "descriptors = [coulomb_matrices, mbtrs, soaps_default, soaps_high]\n", "ndatas = [100, 500, 1000, 2500, 5500]\n", "\n", "learning_curve = np.zeros([len(descriptors),len(ndatas),2])\n", "\n", "for i,X in tqdm(enumerate(descriptors)):\n", " X_train_full, X_test, y_train_full, y_test = train_test_split(\n", " X, y, #inputs and outputs\n", " test_size=0.20, #20% of the data is test data\n", " random_state=42, #a random seed\n", " )\n", " for j,ndata in tqdm(enumerate(ndatas)):\n", "\n", " X_train, y_train = resample(X_train_full, y_train_full, replace=False, n_samples=ndata, random_state=42)\n", " scaler = preprocessing.StandardScaler()\n", " scaler.fit(X_train)\n", " X_train = scaler.transform(X_train)\n", " X_test = scaler.transform(X_test)\n", "\n", " #rescaling the labels\n", " u = np.mean(y_train)\n", " s = np.std(y_train)\n", " y_train_scaled = (y_train-u)/s\n", " y_test_scaled = (y_test-u)/s\n", "\n", " def unscale(y):\n", " # this is for transforming the predictions back into proper energies\n", " return (y*s)+u\n", "\n", " scoring = {\n", " 'r2_score': make_scorer(r2_score),\n", " 'MAE': make_scorer(mean_absolute_error),\n", " 'RMSE': make_scorer(mean_squared_error),\n", " }\n", "\n", " krr = KernelRidge(alpha=0.0001, kernel='rbf', gamma=0.001)\n", " scores = cross_validate(\n", " krr, \n", " X_train, \n", " y_train_scaled, \n", " scoring=scoring, \n", " cv=5, \n", " n_jobs=6,\n", " return_estimator=False,\n", " return_train_score=True\n", " )\n", "\n", " learning_curve[i,j,0] = np.mean(scores['test_MAE'])*s\n", " learning_curve[i,j,1] = np.std(scores['test_MAE'])*s\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "679f3b70", "metadata": {}, "outputs": [], "source": [ "descriptor_names = ['CM', 'MBTR', 'SOAP_default', 'SOAP_high'] #, soaps_default, soaps_high]\n", "\n", "for i, d in enumerate(descriptors):\n", " plt.errorbar(ndatas, learning_curve[i,:,0], learning_curve[i,:,1],label=descriptor_names[i])\n", "plt.xlabel('number of training data points')\n", "plt.ylabel('MAE in eV')\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "f448057c", "metadata": {}, "source": [ "## Part 7: The Challenge" ] }, { "cell_type": "markdown", "id": "0e5bf76f", "metadata": {}, "source": [ "
\n", "\n", "Use the code above and your newly acquired skills to create a KRR model to predict the PBE0 atomization energy and the PBE0-DFTB3 delta energy of the QM7 dataset. \n", "\n", "Try to use a tailored descriptor and a set of hyperparameters that provide you with the lowest possible MAE on the test set.\n", "\n", "- Use a 80%/20% train/test split with `random_state=42` for reproducibility\n", "- You can use any of the descriptors that we have covered with any setting. The model just needs to feasibly train and predict in a reasonable timeframe.\n", "\n", "- Can you beat the MAE achieved in the original paper? (Phys. Rev. Lett. 108, 058301) ca. 10kcal/mol = 0.43 eV\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "8a685b6c", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "8ea72174", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "c8718f62", "metadata": {}, "source": [ "## Conclusion\n", " " ] }, { "cell_type": "markdown", "id": "0ef3991d", "metadata": {}, "source": [ "\n", "This is it. Hopefully, this gave a brief overview over some workflows and methods that are useful for ML applications in chemistry and materials research.\n", "\n", "Where next for deeper insights into methods?\n", "\n", "- Check out the AI lectures and tutorials on the https://nomad-lab.eu/aitoolkitlist webpage\n", "Link: [NOMAD AI Tutorials](https://nomad-lab.eu/prod/v1/gui/analyze/course)\n", "\n", "- \"Quantum Chemistry in the Age of Machine Learning\", edited by Pavlo Dral (2022) with Github tutorials for each chapter\n", "\n", "- Work on your own project. The best way to learn is to work on a problem and to make mistakes!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.12.0" } }, "nbformat": 4, "nbformat_minor": 5 }