8S50x_L16.ipynb

{"cells": [{"cell_type": "markdown", "id": "1ee373c5", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<hr style=\"height: 1px;\">\n", "<i>This notebook was authored by the 8.S50x Course Team, Copyright 2022 MIT All Rights Reserved.</i>\n", "<hr style=\"height: 1px;\">\n", "<br>\n", "\n", "<h1>Lesson 16: Predicting the Momentum of the Tau Using Machine Learning</h1>\n"]}, {"cell_type": "markdown", "id": "41345f3a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_0'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.0 Overview</h2>\n"]}, {"cell_type": "markdown", "id": "fce9b388", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<h3>Navigation</h3>\n", "\n", "<table style=\"width:100%\">\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_16_1\">L16.1 Higgs to Taus</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#exercises_16_1\">L16.1 Exercises</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_16_2\">L16.2 Regression Analysis of Tau Momentum</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#exercises_16_2\">L16.2 Exercises</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_16_3\">L16.3 Reconstructing the Higgs Mass</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#exercises_16_3\">L16.3 Exercises</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_16_4\">L16.4 The Full Mass Regression</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#exercises_16_4\">L16.4 Exercises</a></td>\n", "    </tr>\n", "    <tr>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#section_16_5\">L16.5 Tuning the NN Architecture</a></td>\n", "        <td style=\"text-align: left; vertical-align: top; font-size: 10pt;\"><a href=\"#exercises_16_5\">L16.5 Exercises</a></td>\n", "    </tr>\n", "</table>"]}, {"cell_type": "markdown", "id": "467638f0", "metadata": {"tags": ["learner", "catsoop_00", "md"]}, "source": ["<h3>Learning Objectives</h3>\n", "\n", "In the following sections, we are going to use deep learning regression to solve an important problem in Higgs identification and reconstruction. In particular, what we are going to study is the Higgs decay into two tau ($\\tau$) leptons. Along with the electron and muon, the tau is the third lepton, and was discovered in 1975. Similar to the other two, it comes in positive and negative charges, has spin 1/2, and only interacts via the electroweak forces. However, it is much heavier than the other two, 170 times as heavy as the muon and 3500 times as heavy as the electron, meaning that it is almost twice as heavy as a proton. As a result, its lifetime is very short, about $3\\times 10^{\u221213}$ s, or about 8 million times shorter than that for the muon.\n", "\n", "By correctly measuring and accounting for the total transverse momentum of the two tau particles (the component of their momentum perpendicular to the beam direction), we can reconstruct the Higgs mass. However, the short tau lifetime means that it cannot be detected directly, but can only be reconstructed from its decay products. One difficulty in this reconstruction is that some decay products of tau particles are neutrinos. These are not detectable due to their negligibly small interaction probability and, therefore, their momentum cannot be measured. Machine learning will help us accurately reconstruct the momentum of the tau particles and thereby help us constrain the Higgs mass.\n", "\n", "In the following sections, we will do the following:\n", "\n", "- First, we will describe the Higgs decay products and identify the decay to taus as an important channel.\n", "\n", "- Then, we will use regression analysis to improve our tau momentum reconstruction.\n", "\n", "- Next, we will work towards reconstructing the Higgs mass from our improved tau information.\n", "\n", "- When we have a setup that seems to be able to reconstruct the Higgs mass, the final challenge will be to use all the information we have in our event to do the best possible reconstruction. Ultimately, we will use as many different event and particle variables as we can from our data. Machine learning will help us construct the appropriate weights.\n", "\n", "- Finally, we will consider other neural network architectures in our solution.\n"]}, {"cell_type": "markdown", "id": "116c757a", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Installing Tools</h3>\n", "\n", "Before we do anything, let's make sure we install the tools we need."]}, {"cell_type": "code", "execution_count": null, "id": "d59f3ad0", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.0-runcell00\n", "\n", "!pip install torch\n", "!pip install imageio\n", "!pip install george\n", "!pip install uproot\n", "!pip install awkward\n", "!pip install pylorentz"]}, {"cell_type": "markdown", "id": "dff0575f", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Importing Libraries</h3>\n", "\n", "Before beginning, run the cell below to import the relevant libraries for this notebook. "]}, {"cell_type": "code", "execution_count": null, "id": "0f3a0154", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.0-runcell01\n", "\n", "import torch                        #https://pytorch.org/docs/stable/torch.html\n", "import torch.nn as nn               #https://pytorch.org/docs/stable/nn.html\n", "from torch.autograd import Variable #https://pytorch.org/docs/stable/autograd.html\n", "import torch.nn.functional as F     #https://pytorch.org/docs/stable/nn.functional.html\n", "import torch.utils.data as Data     #https://pytorch.org/docs/stable/data.html\n", "\n", "import matplotlib.pyplot as plt     #https://matplotlib.org/3.5.3/api/_as_gen/matplotlib.pyplot.html\n", "%matplotlib inline\n", "\n", "from pylorentz import Momentum4     #https://pypi.org/project/pylorentz/\n", "\n", "import numpy as np                  #https://numpy.org/doc/stable/\n", "import imageio                      #https://imageio.readthedocs.io/en/stable/"]}, {"cell_type": "markdown", "id": "3d515076", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Setting Default Figure Parameters</h3>\n", "\n", "The following code cell sets default values for figure parameters."]}, {"cell_type": "code", "execution_count": null, "id": "97d0edf8", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.0-runcell02\n", "\n", "#set plot resolution\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "#set default figure parameters\n", "plt.rcParams['figure.figsize'] = (9,6)\n", "\n", "medium_size = 12\n", "large_size = 15\n", "\n", "plt.rc('font', size=medium_size)          # default text sizes\n", "plt.rc('xtick', labelsize=medium_size)    # xtick labels\n", "plt.rc('ytick', labelsize=medium_size)    # ytick labels\n", "plt.rc('legend', fontsize=medium_size)    # legend\n", "plt.rc('axes', titlesize=large_size)      # axes title\n", "plt.rc('axes', labelsize=large_size)      # x and y labels\n", "plt.rc('figure', titlesize=large_size)    # figure title"]}, {"cell_type": "markdown", "id": "e198f21b", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_1'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.1 Higgs to Taus</h2>  \n", "\n", "| [Top](#section_16_0) | [Previous Section](#section_16_5) | [Exercises](#exercises_16_1) | [Next Section](#section_16_2) |\n"]}, {"cell_type": "markdown", "id": "00dcb7cd", "metadata": {"tags": ["learner", "md"]}, "source": ["*The material in this section is discussed in the video **<a href=\"https://courses.mitxonline.mit.edu/learn/course/course-v1:MITxT+8.S50.2x+1T2025/block-v1:MITxT+8.S50.2x+1T2025+type@sequential+block@seq_LS16/block-v1:MITxT+8.S50.2x+1T2025+type@vertical+block@vert_LS16_vid1\" target=\"_blank\">HERE</a>.** You are encouraged to watch that video and use this notebook concurrently.*"]}, {"cell_type": "markdown", "id": "4b051e7e", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Slides</h3>\n", "\n", "Run the code below to view the slides for this section, which are discussed in the related video. You can also open the slides in a separate window <a href=\"https://mitx-8s50.github.io/slides/L19/slides_L19_06.html\" target=\"_blank\">HERE</a>."]}, {"cell_type": "code", "execution_count": null, "id": "0bcbe36c", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.1-slides\n", "\n", "from IPython.display import IFrame\n", "IFrame(src='https://mitx-8s50.github.io/slides/L19/slides_L19_06.html', width=970, height=550)"]}, {"cell_type": "markdown", "id": "f44ff55b", "metadata": {"tags": ["learner", "md", "lect_01"]}, "source": ["<h3>Overview</h3>\n", "\n", "The Higgs boson decays to many different final states which contain many different types of particles. The most likely Higgs decay is into b-quarks. The 4th most likely final state is where the Higgs decays into tau leptons. Taus are effectively heavy electrons or muons, and behave in similar ways to electrons and muons. Analyzing the production of taus can provide valuable information that helps constrain the mass of the Higgs because their probability of production from a Higgs decay is directly proportional to the Higgs mass."]}, {"cell_type": "markdown", "id": "b477036e", "metadata": {"tags": ["learner", "md", "lect_01"]}, "source": ["<h3>Reconstructing the tau transverse energy</h3>\n", "\n", "The important aspect of taus that we will focus on in this lecture is that a tau can decay to neutrinos, and either an electron and a neutrino, a muon and a neutrino, or quarks. Electrons, muons, and quarks can be reconstructed in a collider detector. Note that the quarks are actually detected via jets (sprays of closely spaced particles) that they produce. However, the neutrinos cannot because they are very weakly interacting and, therefore, they will just pass through the detector without leaving any trace. \n", "\n", "The one constraint we have for neutrinos is what we call missing transverse energy (not discussed in the video, but see slides 7-9). Missing transverse energy is a constraint that originates from the fact that when we collide protons along a specific axis, the momentum transverse to that axis is conserved. In other words, for all of the particles created in an event: \n", "\n", "$$\n", "\\begin{equation}\n", "\\sum_{i=1}^{N}\\vec{p_{T,i}} = 0 \\\\\n", "\\sum_{i\\neq {\\rm neutrino}}^{N}\\vec{p_{T,i}} + \\sum_{i= {\\rm neutrino}}^{N}\\vec{p_{T,i}} = 0 \\\\\n", "\\sum_{i= {\\rm neutrino}}^{N}\\vec{p_{T,i}} = -\\sum_{i\\neq {\\rm neutrino}}^{N}\\vec{p_{T,i}} \\\\\n", "\\rm{\\vec{MET}}  = -\\sum_{i\\neq {\\rm neutrino}}^{N}\\vec{p_{T,i}} \\\\\n", "\\end{equation}\n", "$$\n", "\n", "This is what we call the MET or missing transverse energy The seemingly confusing use of \"energy\" comes about because the energy and momentum of neutrinos only differ by a simply multiplicative factor of the speed of light. It is very common in high energy physics to carelessly use the terms \"energy\" and \"momentum\" interchangeably. The MET is our only constraint on the neutrinos and it gives us a vector in the plane perpendicular to the collision axis along the beam.  \n", "\n", "As is the case with other decay reconstructions, for example W and Z bosons, the best way to find the Higgs decaying to taus is to reconstruct the mass and look for a bump. The challenge then is to first identify the taus, and, second, reconstruct the mass of the Higgs from this decay. "]}, {"cell_type": "markdown", "id": "c9aa1812", "metadata": {"tags": ["learner", "md", "lect_01"]}, "source": ["<h3>Role of machine learning</h3>\n", "\n", "As discussed previously, the Higgs actually decays to 2 taus, which means that we will get anywhere from 2$-$4 neutrinos. What we want to do in our machine learning algorithm is to identify the topology, and return back a 4-vector that is the sum of all the neutrinos together.  \n", "\n", "In the next section, we will load a dataset that we can use to study this problem."]}, {"cell_type": "markdown", "id": "d6ad3f5d", "metadata": {"tags": ["learner", "md"]}, "source": ["<a name='exercises_16_1'></a>     \n", "\n", "| [Top](#section_16_0) | [Restart Section](#section_16_1) | [Next Section](#section_16_2) |\n"]}, {"cell_type": "markdown", "id": "24b1aca7", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.1.1</span>\n", "\n", "Our objective is to reconstruct the tau momentum, and thereby probe the Higgs mass. First, we must consider the tau decay products, some of which are not visible (i.e. neutrinos).\n", "\n", "If we write down all possible tau decays, including decays to electrons, muons, charged and neutral pions, there are as many as 10 different types of decays, all with similar probabilities of happening. Each decay produces neutrinos going in a slightly different direction.  Why would machine learning (ML) be a good way to determine the directions of the neutrinos? (Also, note, we can simulate this whole process well.)\n", "\n", "A) ML can figure out the exact decay.\\\n", "B) ML is much faster than a rule based algorithm and speed is critical here.\\\n", "C) ML can come up with a weighted decision for the probability of each decay and, based on this weight and knowledge of the decays, determine the most likely MET (missing transverse energy).\\\n", "D) ML regression can simulate all decays and choose the best.\n"]}, {"cell_type": "markdown", "id": "ba8bca06", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_2'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.2 Regression Analysis of Tau Momentum</h2>  \n", "\n", "| [Top](#section_16_0) | [Previous Section](#section_16_1) | [Exercises](#exercises_16_2) | [Next Section](#section_16_3) |\n"]}, {"cell_type": "markdown", "id": "44ee885b", "metadata": {"tags": ["learner", "md"]}, "source": ["*The material in this section is discussed in the video **<a href=\"https://courses.mitxonline.mit.edu/learn/course/course-v1:MITxT+8.S50.2x+1T2025/block-v1:MITxT+8.S50.2x+1T2025+type@sequential+block@seq_LS16/block-v1:MITxT+8.S50.2x+1T2025+type@vertical+block@vert_LS16_vid2\" target=\"_blank\">HERE</a>.** You are encouraged to watch that video and use this notebook concurrently.*"]}, {"cell_type": "markdown", "id": "aefc9099", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Slides</h3>\n", "\n", "Run the code below to view the slides for this section, which are discussed in the related video. You can also open the slides in a separate window <a href=\"https://mitx-8s50.github.io/slides/L19/slides_L19_07.html\" target=\"_blank\">HERE</a>."]}, {"cell_type": "code", "execution_count": null, "id": "119ca253", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-slides\n", "\n", "from IPython.display import IFrame\n", "IFrame(src='https://mitx-8s50.github.io/slides/L19/slides_L19_07.html', width=970, height=550)"]}, {"cell_type": "markdown", "id": "65f7e19e", "metadata": {"tags": ["md", "lect_02", "learner"]}, "source": ["<h3>Overview</h3>\n", "\n", "In this section, we will work towards reconstructing the tau transverse momentum."]}, {"cell_type": "markdown", "id": "00aadbe6", "metadata": {"tags": ["md", "lect_02", "learner"]}, "source": ["<h3>Using Uproot to Process Data</h3>\n", "\n", "We are going to use the `uproot` package to process this simulated tau dataset.\n", "\n", "Uproot is a software program aimed at processing so-called \"root\" files. The Root package includes a very flexible data format, as well as extensive analysis and plotting capabilities. It is extensively used by high energy physicists and others to analyze data. While you really don't need to do know much about this data format, we should at least look at the variables in the dataset to understand what is going on. \n", "\n", "Let's go ahead and load data, and then print out the variables to see what the data includes. "]}, {"cell_type": "markdown", "id": "01afb3f9", "metadata": {"tags": ["md", "lect_02", "learner", "learner_chopped"]}, "source": ["<h3>Data</h3>\n", "\n", ">description: Higgs to tau dataset decaying to to tau leptons<br>\n", ">source: https://zenodo.org/record/8035277 <br>\n", ">attribution: Philip Harris (CMS Collaboration), DOI:10.5281/zenodo.8035277 "]}, {"cell_type": "code", "execution_count": null, "id": "58bc03db", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell01\n", "\n", "# NOTE: these files are too large to include in the original repository,\n", "# so you must download them from the sources below\n", "#\n", "# Ways to download:\n", "#     1. Copy/paste the link (replace =0 with =1 to download automatically)\n", "#     2. Use the wget commands below (works in Colab, but you may need to install wget if using locally)\n", "#\n", "# Location of files:\n", "#     Move the files to the directory data/L16\n", "#\n", "# Using wget: (works in Colab)\n", "#     Upon downloading, the code below will move them to the appropriate directory\n", "\n", "#get the data\n", "!wget -P data/L16/ https://www.dropbox.com/s/csgx8t35i3un9kr/Regression2.root?dl=0\n", "!mv data/L16/Regression2.root?dl=0 data/L16/Regression2.root"]}, {"cell_type": "code", "execution_count": null, "id": "76b590ca", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell02\n", "\n", "import uproot\n", "from collections import OrderedDict \n", "reg    = uproot.open(\"data/L16/Regression2.root\")[\"Tree\"]"]}, {"cell_type": "code", "execution_count": null, "id": "daf2fef8", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell03\n", "\n", "#what are the inputs\n", "print(reg.keys())\n", "cut=reg['genpt1'].array() >  100\n", "vals=reg['genpt1'].array(library=\"np\")[cut]\n", "np.histogram(vals)"]}, {"cell_type": "markdown", "id": "2e6de9dd", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["<h3>About the Data</h3>\n", "\n", "This is a bit of a complicated dataset.  What it includes are Higgs events with their decays into two taus. Thus, `genpt1` or `genpt2` refer to tau 1 or tau 2. The identity `gen` refers to \"true\" information that is generated from simulation.\n", "\n", "The taus then decay into neutrinos and visible particles. In light of that, what we show above are the top 5 visible particles from tau decays. So, `part1id5` correspons to particle 5 that decays from tau 1.\n", "\n", "In addition to the \"true\" generated data, we have reconstructed (or observed) data, identified by `reco`. These are the reconstructed properties of the tau particles, based on the visible decay products. It is what would be actually observed in the detector."]}, {"cell_type": "markdown", "id": "d541aab9", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["<h3>A First Plot</h3>\n", "\n", "Let's begin by plotting the histograms of two datasets.\n", "\n", "In red, we plot the reconstructed transverse momentum (pT) of one of the taus, `recopt1`, which we synonymously call the observed pT.\n", "\n", "In black we plot the true, generated transverse momentum of that same tau, `genpt1`.\n", "\n", "Note, the true transverse momentum of tau1 (black) will peak at half the Higgs mass, because the Higgs decays to two taus. This is the peak at about 62.5."]}, {"cell_type": "code", "execution_count": null, "id": "de19a875", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell04\n", "\n", "import numpy as np\n", "\n", "def plot(iVar,iMin,iMax,iColor,iLabel): \n", "    mask=(reg[iVar].array() > 0)\n", "    data=reg[iVar].array(library=\"np\")[mask]\n", "    counts, binEdges = np.histogram(data,bins=50,range=(iMin,iMax),density=False)\n", "    binCenters = (binEdges[1:]+binEdges[:-1])*.5\n", "    err = np.sqrt(counts)\n", "    plt.errorbar(binCenters, counts, yerr=err,fmt=\"o\",alpha=0.5,c=iColor,label=iLabel, ms=3)\n", "    plt.xlabel(\"mass\")\n", "    plt.ylabel(\"N events\")\n", "    return binCenters,counts,err\n", "\n", "plot(\"genpt1\",0,200,\"black\",\"true $p_{T}$\")\n", "plot(\"recopt1\",0,200,\"red\",\"observed $p_{T}$\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "3a07d5b9", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["<h3>The Goal</h3>\n", "\n", "Our goal is to go from the points shown in red to the black distribution. Additionally, our goal will be to predict the other directions as well (i.e., the full 4-vector).\n", "\n", "What we saw above is the transverse momentum, this is the momentum projected onto the plane transverse to the collision. In reality, we want to predict the full 4 vector of the tau decay. This means that we will need to predict the total momentum 3-vector of the tau, and then knowing the mass of the tau, we can calculate the energy, and hence we have the full 4-vector. \n", "\n", "To understand how we go from the reconstructed (observed) vector to the true (generated) one, let's first go ahead and make a 2D plot of the tau true pT correlated with its reconstructed value. "]}, {"cell_type": "code", "execution_count": null, "id": "0da6dd45", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell05\n", "\n", "#To visualize the whole problem let's make a 2D plot\n", "mask=np.logical_and(reg[\"genpt1\"].array() > 0, reg[\"recopt1\"].array()>0)\n", "x=reg[\"genpt1\"].array(library=\"np\")[mask]\n", "y=reg[\"recopt1\"].array(library=\"np\")[mask]\n", "plt.xlabel(\"gen $p_{T}$\")\n", "plt.ylabel(\"reco $p_{T}$\")\n", "plt.hist2d(x,y,bins=200)\n", "plt.xlim(0,100)\n", "plt.ylim(0,100)\n", "plt.show()\n", "\n", "\n", "print(\"Pre Correlation:\",np.corrcoef(y.flatten(),x.flatten())[0][1])"]}, {"cell_type": "markdown", "id": "1417d8f7", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["<h3>Defining the Network</h3>\n", "\n", "So, what we see is that the reconstructed momentum is not directly correlated with the generator momentum, there are some fluctuations around the true value. What this means is that we will eventually need more variables to actually get the true momentum of the tau system. We can't just predict true$=f(reco)$ for one variable, we need to predict the full vector with many variables as input. \n", "\n", "With all of this information, let's go ahead and build a regression. We can start from our previous setup, and try to just take in the reconstructed momentum, and predict the true momentum. \n", "\n", "Even though we know it's not optimal, let's go ahead and train a simple network to predict the true transverse momentum, `genpt1`, given the input reconstructed transverse momentum, `recopt1`. To do that, we will need to convert things into pytorch format. "]}, {"cell_type": "code", "execution_count": null, "id": "0898d46f", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell06\n", "\n", "#Let's prepare the data to be pytorch friendly\n", "x=torch.from_numpy(reg[\"recopt1\"].array(library=\"np\")[mask].reshape(len(reg[\"recopt1\"].array(library=\"np\")[mask]),1))\n", "y=torch.from_numpy(reg[\"genpt1\"].array(library=\"np\")[mask].reshape(len(reg[\"genpt1\"].array(library=\"np\")[mask]),1))\n", "x, y = Variable(x), Variable(y)\n", "#torch_dataset = Data.TensorDataset(x, y)\n", "#loader = Data.DataLoader(dataset=torch_dataset, batch_size=BATCH_SIZE, shuffle=True, num_workers=2,)\n", "\n", "#Now let's make a simple model\n", "torch.manual_seed(1)    # reproducible\n", "net = torch.nn.Sequential(\n", "        torch.nn.Linear(1, 200),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(200, 100),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(100, 1),\n", "    )\n", "optimizer = torch.optim.Adam(net.parameters(), lr=0.01)\n", "loss_func = torch.nn.MSELoss()"]}, {"cell_type": "markdown", "id": "8dd01059", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["Now, in the interest of looking at how our neural network evolves over time as we train, what we can do is make a plot of the neural output momentum compared to the true value. Below we define the `makePlot` and `train` functions."]}, {"cell_type": "code", "execution_count": null, "id": "d7ab80c6", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell07\n", "\n", "def makePlot(x,y,prediction,ax,fig,images,t,loss,ymin,ymax):\n", "    # plot and show learning process\n", "    plt.cla()\n", "    ax.set_title('Regression Analysis', fontsize=35)    \n", "    ax.set_xlabel('Mass', fontsize=24)\n", "    ax.set_ylabel('N', fontsize=24)\n", "    ax.hist(prediction.data.numpy(),color=\"red\",bins=20,range=(0,200),alpha=0.5,label='pred')\n", "    ax.hist(y.data.numpy(),color=\"black\" ,bins=20,range=(0,200),alpha=0.5,label='gen')\n", "    ax.hist(x.data.numpy(),color=\"green\",bins=20,range=(0,200),alpha=0.5,label='reco')\n", "    #ax.scatter(x.data.numpy(), y.data.numpy(), color = \"orange\")\n", "    #ax.plot(x.data.numpy(), prediction.data.numpy(), 'g-', lw=3)\n", "    ax.text(100, 2000, 'Epoch = %d' % t,fontdict={'size': 24, 'color':  'red'})\n", "    ax.text(100, 5000, 'Loss = %.4f' % loss.data.numpy(),fontdict={'size': 24, 'color':  'red'})\n", "    fig.canvas.draw()       # draw the canvas, cache the renderer\n", "    ax.legend()\n", "    image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')\n", "    image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))\n", "    images.append(image)\n", "    \n", "def train(x,y,net,loss_func,opt,nepochs,ymin,ymax):\n", "    images = []\n", "    fig, ax = plt.subplots(figsize=(12,7))\n", "    for epoch in range(nepochs):\n", "        if epoch % 50 == 0: \n", "            print(\"epoch:\",epoch)\n", "        prediction = net(x)\n", "        loss = loss_func(prediction, y) \n", "        opt.zero_grad()\n", "        loss.backward() \n", "        optimizer.step()\n", "        if epoch % 4 == 0:\n", "            makePlot(x,y,prediction,ax,fig,images,epoch,loss,ymin,ymax)\n", "    return images"]}, {"cell_type": "markdown", "id": "83b60131", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["<h3>Training</h3>\n", "\n", "Run the code below to train the neural network over 150 epochs.\n", "\n", "Like we in previous Lessons, we will plot our performance over time as a .gif. We show histograms of the generated (true) data (grey), the reconstructed (observed) data (green), and the NN prediction (red)."]}, {"cell_type": "code", "execution_count": null, "id": "c112f725", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell08\n", "\n", "from IPython.display import Image\n", "images=train(x,y,net,loss_func,optimizer,150,0,1)\n", "torch.save(net.state_dict(), 'data/L16/tau_pt_simple.pt')\n", "imageio.mimsave('data/L16/reg_1.gif', images, fps=12)\n", "Image(open('data/L16/reg_1.gif','rb').read())\n", "\n", "# Note: you may wish to visualize the results as open histograms, rather than filled-in.\n", "# Here, we will keep the same styling as the related video."]}, {"cell_type": "markdown", "id": "b35a06aa", "metadata": {"tags": ["learner", "md", "lect_02"]}, "source": ["So, we can see that the NN tries to predict the average, but it's not that great. Let's look in more detail at the properties of this. Specifically, in the plot below, we will show the following:\n", "\n", "- shown in red (labeled `true`) is the ratio of the full, generated tau momentum (including the neutrino) divided by the reconstructed tau momentum, with the latter based on a sum of the visible tau components.\n", "\n", "- shown in blue (labeled `corr`) is the ratio of the full, generated tau momentum (again) divided by the predicted momentum as found by the neural network. \n", "\n", "\n", "Note, a perfect prediction would exhibit a strong peak at 1. The red histogram shows the inability of the reconstruction to predict the true tau momentum. The blue histogram shows more peakiness than the red, demonstrating that the NN is an improvement."]}, {"cell_type": "code", "execution_count": null, "id": "2814e479", "metadata": {"tags": ["learner", "py", "lect_02", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.2-runcell09\n", "\n", "true=reg[\"genpt1\"].array(library=\"np\")[mask]\n", "reco=reg[\"recopt1\"].array(library=\"np\")[mask]\n", "pred=net(x)\n", "ratio=np.array(true/reco)\n", "ratiopred=y/pred\n", "plt.hist(ratio,color=\"red\",bins=20,range=(0,3),alpha=0.5,label=\"true\")\n", "plt.hist(ratiopred.data.numpy(),color=\"blue\",bins=20,range=(0,3),alpha=0.5,label=\"corr\")\n", "plt.legend()\n", "plt.show()\n", "print(\"True Mean: \",ratio.mean(),\"True StdDev:\",ratio.std())\n", "print(\"NN Mean: \",ratiopred.data.numpy().mean(),\"NN StdDev:\",ratiopred.data.numpy().std())"]}, {"cell_type": "markdown", "id": "00135026", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='exercises_16_2'></a>     \n", "\n", "| [Top](#section_16_0) | [Restart Section](#section_16_2) | [Next Section](#section_16_3) |\n"]}, {"cell_type": "markdown", "id": "49d0dd82", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.2.1</span>\n", "\n", "Consider the preceding plot, generated by code cell `L16.2-runcell09`. Before we discuss the neural network prediction, why does the ratio of the generated over reconstructed tau momentum (labeled `true`) exhibit a large tail extending to values bigger than 1?\n", "\n", "A) Systematic uncertainties in the measurement contribute to variations in the momentum estimation, leading to a large tail in the reconstructed tau momentum.\\\n", "B) The reconstruction process underestimates the total momentum (on average) because it only accounts for the momentum of the visible components, ignoring the momentum of any neutrinos.\\\n", "C) The reconstruction process overestimates the total momentum (on average) because it only accounts for the momentum of the visible components, ignoring the momentum of any neutrinos.\\\n", "D) The reconstruction process accounts for too many neutrinos, thus overestimating the total momentum (on average).\\\n", "E) The reconstruction process accounts for too many neutrinos, thus underestimating the total momentum (on average)."]}, {"cell_type": "markdown", "id": "9cab2351", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.2.2</span>\n", "\n", "The blue histogram generated by code cell `L16.2-runcell09` shows the ratio of the true momentum divided by that found by the neural network. Why does the neural network output yield an average ratio that is much closer to 1, and why is the width (i.e. the standard deviation of the ratio distribution) much smaller? Select all that apply:\n", "\n", "A) The neural network brings the average correction to 1 because (on average) it accurately predicts the neutrino momentum.\\\n", "B) The neural network brings the average correction to 1 because it learns to always predict the neutrino momentum perfectly.\\\n", "C) The width of the ratio distribution is partly an indication of how well the neural network can predict the neutrino momentum.\\\n", "D) As the neural network becomes more complex and incorporates additional inputs, the width of the distribution should become narrower, indicating improved prediction of the neutrino momentum.\n", "\n"]}, {"cell_type": "markdown", "id": "97cfd698", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.2.3</span>\n", "\n", "Code cell `L16.2-runcell05` generated a 2D plot of the correlation between the reconstructed and generated tau momentum and also calculated their correlation coefficient. Complete the code cell below to calculate the correlation coefficient for the momentum predicted by the neural network versus the generated value.\n", "\n", "What is the correlation of the neural network output with the true value? How does this compare to the correlation before the NN correction? Enter your answer as a list of numbers with precision 1e-2: `[corr-original, corr-NN]`\n"]}, {"cell_type": "code", "execution_count": null, "id": "1486329d", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>EXERCISE: L16.2.3\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "ytmp=y.detach().numpy()\n", "ptmp=pred.detach().numpy()\n", "print(\"Pre Correlation:\", #YOUR CODE HERE)\n", "print(\"NN Correlation:\", #YOUR CODE HERE)"]}, {"cell_type": "markdown", "id": "f846027c", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_3'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.3 Reconstructing the Higgs Mass</h2>  \n", "\n", "| [Top](#section_16_0) | [Previous Section](#section_16_2) | [Exercises](#exercises_16_3) | [Next Section](#section_16_4) |\n"]}, {"cell_type": "markdown", "id": "0f4e2c6b", "metadata": {"tags": ["learner", "md"]}, "source": ["*The material in this section is discussed in the video **<a href=\"https://courses.mitxonline.mit.edu/learn/course/course-v1:MITxT+8.S50.2x+1T2025/block-v1:MITxT+8.S50.2x+1T2025+type@sequential+block@seq_LS16/block-v1:MITxT+8.S50.2x+1T2025+type@vertical+block@vert_LS16_vid3\" target=\"_blank\">HERE</a>.** You are encouraged to watch that video and use this notebook concurrently.*"]}, {"cell_type": "markdown", "id": "b91d2bf4", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["<h3>Observing the improvements</h3>\n", "\n", "Ok, so we now have a very basic regression that can predict the change in momentum due to a single neutrino from a single tau decay. This is a very simple regression, and we see that the prediction does improve things, but perhaps doesn't work as well as we might have hoped. Also, it's a bit hard to understand exactly what the training is doing. \n", "\n", "In light of this, let's write go back to the samples and actually reconstruct the Higgs mass. Then, we can try the training again and see if we can really optimize it. \n", "\n", "If we look at the dataset, we find that we have two variables that give the mass. One is the \"True Higgs mass\" (`hmass`), which is the Higgs mass that is created by our simulation by sampling a Higgs mass distribution for each event. \n", "\n", "For the second mass, what the code that created this dataset did is to take the two reconstructed tau objects (without any neutrinos) and compute the reconstructed mass of the Higgs (`recohmass`) by adding the two tau 4-vectors. In other words, the mass can be written as a sum of the reconstructed 4-vectors of the two identified and reconstructed taus: \n", "\n", "$$\n", "\\begin{equation}\n", "m_{\\rm vis} = \\left(E_{1} + E_{2} \\right)^2 - \\left(\\vec{p}_{1} + \\vec{p}_{2}\\right)^2\n", "\\end{equation}\n", "$$\n", "\n", "We call this the visible mass, for the simple fact that it excludes all the invisible particles (aka neutrinos). We will forgo the details of how the 4-vectors are reconstructed. However, just assume it's a complicated process, for which we don't need to know the details right now. Ok, let's plot the mass. \n"]}, {"cell_type": "code", "execution_count": null, "id": "a6d6fa0d", "metadata": {"tags": ["learner", "py", "lect_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.3-runcell01\n", "\n", "#now let's construct the Higgs mass\n", "plot(\"hmass\",0,200,\"black\",\"True Higgs Mass\")\n", "plot(\"recohmass\",0,200,\"red\",\" Reconstructed\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "b7d0a081", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["This plot is a bit hard to read because the generated Higgs mass distribution is quite narrow, resulting in a very high value in one bin, while the reconstructed values are much more spread out. Now, to incorporate the effect of a correction to the tau momentum, we first need to reproduce the calculation that generated a Higgs mass from the reconstructed tau 4-vectors. To do this, we are going to use a pylorentz vector, which just does the calculation shown in the equation above. \n", "\n", "All the comparison plotted by code cell `L16.3-runcell02` shows is that our recalculation of the visible Higgs mass using the reconstructed tau information is an exact match to the values found in the dataset. So, we can be confident in using this code to find a Higgs mass from the corrected tau momenta given by the neural network."]}, {"cell_type": "code", "execution_count": null, "id": "66a8faea", "metadata": {"tags": ["learner", "py", "lect_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.3-runcell02\n", "\n", "#!pip install pylorentz\n", "from pylorentz import Momentum4\n", "\n", "#Let's compute the mass on the fly\n", "#Momentum4 calculates the 4-vector, taking as input the mass, eta and phi angles, and the transverse momentum (pT) \n", "def masscompute(iVec1,iVec2):\n", "    tau_1 = Momentum4.m_eta_phi_pt(iVec1[3], iVec1[1], iVec1[2], iVec1[0])\n", "    tau_2 = Momentum4.m_eta_phi_pt(iVec2[3], iVec2[1], iVec2[2], iVec2[0])\n", "    return (tau_1+tau_2).m\n", "    \n", "def hmass(massfunc):\n", "    mask=(reg[\"recohmass\"].array() > 0)\n", "    varlist=[\"recopt1\",\"recoeta1\",\"recophi1\",\"recomass1\",\"recopt2\",\"recoeta2\",\"recophi2\",\"recomass2\"]\n", "    arr=0\n", "    idx=0\n", "    for x in varlist:\n", "        pArr=reg[x].array(library=\"np\")[mask]\n", "        if idx == 0: \n", "            arr = pArr\n", "            idx = idx + 1\n", "        else:\n", "            arr=np.vstack((arr,pArr))\n", "    arr = arr.T\n", "    massc = lambda iarr: massfunc(iarr[0:4],iarr[4:8]) \n", "    hmasses = np.array([massc(p) for p in arr])\n", "    return hmasses\n", "\n", "rawmvis=hmass(masscompute)\n", "plt.title('Reconstructed Mass')\n", "plt.hist(rawmvis,bins=50,range=(0,200),color='blue',label=\"New mass calculation\")\n", "#plot(\"hmass\",0,200,\"black\",\"True Higgs Mass\")\n", "plot(\"recohmass\",0,200,\"red\",\"Mass from dataset\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "cca05edf", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["<h3>Including the Neural Net Computation</h3>\n", "\n", "Ok, great, looks like we can reconstruct the mass. Let's take a look at our how our regression mass looks like. To do this, we need to modify our mass computation to include the neural net momenta."]}, {"cell_type": "code", "execution_count": null, "id": "130d2efd", "metadata": {"tags": ["learner", "py", "lect_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.3-runcell03\n", "\n", "#Let's compute the mass on the fly\n", "def masscompute(iVec1,iVec2):\n", "    pt1 = torch.tensor([iVec1[0]])\n", "    pt2 = torch.tensor([iVec2[0]])\n", "    corr1 = net(pt1).data.numpy()\n", "    corr2 = net(pt2).data.numpy()\n", "    # Here, we replace the reconstructed momentum with the neural-net correcged value\n", "    tau_1 = Momentum4.m_eta_phi_pt(iVec1[3], iVec1[1], iVec1[2], corr1)\n", "    tau_2 = Momentum4.m_eta_phi_pt(iVec2[3], iVec2[1], iVec2[2], corr2)\n", "    return (tau_1+tau_2).m\n", "\n", "rawmass=hmass(masscompute)\n", "plt.hist(rawmass,bins=50,range=(0,200),color='blue',label=\"Neural net corrected mass\")\n", "#plot(\"hmass\",0,200,\"black\",\"True Higgs Mass\")\n", "plot(\"recohmass\",0,200,\"red\",\"Reconstructed mass\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "b928996d", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["Not surprisingly, using the corrected tau momenta found by the neural network makes the calculated Higgs mass distribution narrower!"]}, {"cell_type": "markdown", "id": "cf3ddd19", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["<h3>Optimized Target</h3>\n", "\n", "Now, before we go ahead and try to make a more sophisticated deep learning regression, let's explore an optimizied target. \n", "\n", "With our first neural network in this Lesson, we trained for $p_{T}^{Gen}$ of the neutrino. However, we could imagine instead of training based on a dimensionless variable: \n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{p_{T}^{Gen}}{p_{T}^{Reco}}\n", "\\end{equation}\n", "$$\n", "\n", "The advantage of this alternative target is that it is unitless, which means that if there is an invariance for the momentum (for example, Lorentz invariance), we can capture it in this quantity. Defining losses which respect invariant quantities in physics has become an industry over the past few years. Let's go ahead and repeat our training from before, but now with our training performance code updated to reconstruct the mass. \n"]}, {"cell_type": "code", "execution_count": null, "id": "dfe71127", "metadata": {"tags": ["learner", "py", "lect_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.3-runcell06\n", "\n", "#Now let's make a simple model\n", "torch.manual_seed(1)    # reproducible\n", "net = torch.nn.Sequential(\n", "        torch.nn.Linear(1, 200),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(200, 100),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(100, 1),\n", "    )\n", "optimizer = torch.optim.Adam(net.parameters(), lr=0.01)\n", "loss_func = torch.nn.MSELoss()\n", "x=torch.from_numpy(reg[\"recopt1\"].array(library=\"np\")[mask].reshape(len(reg[\"recopt1\"].array(library=\"np\")[mask]),1))\n", "y=torch.from_numpy(reg[\"genpt1\"].array(library=\"np\")[mask].reshape(len(reg[\"genpt1\"].array(library=\"np\")[mask]),1))\n", "x, y = Variable(x), Variable(y)\n", "ratio=torch.div(y,x)\n", "y=ratio\n", "\n", "#Let's compute the mass on the fly\n", "def masscomputeNN(iVec1,iVec2):\n", "    pt1 = torch.tensor([iVec1[0]])\n", "    pt2 = torch.tensor([iVec2[0]])\n", "    corr1 = net(pt1).data.numpy()*iVec1[0]\n", "    corr2 = net(pt2).data.numpy()*iVec2[0]\n", "    tau_1 = Momentum4.m_eta_phi_pt(iVec1[3], iVec1[1], iVec1[2], corr1)\n", "    tau_2 = Momentum4.m_eta_phi_pt(iVec2[3], iVec2[1], iVec2[2], corr2)\n", "    return (tau_1+tau_2).m\n", "\n", "#And let's plot the mass instead of the pT\n", "def makePlot(x,y,prediction,ax,fig,images,t,loss,ymin,ymax):\n", "    #compute the mass\n", "    rawmass=hmass(masscomputeNN)\n", "    # plot and show learning process\n", "    plt.cla()\n", "    ax.hist(rawmvis,bins=40,range=(0,250),color='blue',alpha=0.5)#,label='raw')\n", "    ax.hist(rawmass,bins=40,range=(0,250),color='red',alpha=0.5)#,label='regressed')\n", "    ax.text(150, 300, 'Epoch = %d' % t,fontdict={'size': 24, 'color':  'red'})\n", "    ax.text(150, 600, 'Loss = %.4f' % loss.data.numpy(),fontdict={'size': 24, 'color':  'red'})\n", "    ax.set_title('Regression Analysis', fontsize=35)\n", "    ax.set_xlabel('Mass', fontsize=24)\n", "    ax.set_ylabel('N', fontsize=24)\n", "    ax.set_ylim(0,2000)\n", "    fig.canvas.draw()       # draw the canvas, cache the renderer\n", "    #ax.legend()\n", "    image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')\n", "    image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))\n", "    images.append(image)"]}, {"cell_type": "markdown", "id": "f4c07ea8", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["Now, let's run the training. The code cell below runs 500 epochs, so it will take some time to finish (15min-1h, depending on available computing resources). Feel free to reduce the number of epochs or go and do something else for a while. Another option is to reduce the number of epochs and then run the code repeatedly (but, in this case, you would not be able to see the full .gif beginning from the initial state)."]}, {"cell_type": "code", "execution_count": null, "id": "43c5847c", "metadata": {"tags": ["learner", "py", "lect_03", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.3-runcell07\n", "\n", "#NOTE: if training does not complete due to timeout in Colab,\n", "#reduce the number of epochs to 250 and run this cell twice,\n", "#or reduce to 125 and run this cell four times,\n", "#for a total of 500 training epochs\n", "images=train(x,y,net,loss_func,optimizer,500,0,1)\n", "torch.save(net.state_dict(), 'data/L16/tau_pt_ratio.pt')\n", "imageio.mimsave('data/L16/reg2_long.gif', images, fps=12)\n", "Image(open('data/L16/reg2_long.gif','rb').read())"]}, {"cell_type": "markdown", "id": "17b05e5f", "metadata": {"tags": ["learner", "md", "lect_03"]}, "source": ["<h3>Challenge Questions</h3>\n", "\n", "After how many epochs does the network stop changing? Does it make sense to train for more epochs (say 1000 total)?\n", "\n", "Repeat the above regression, but now make the true mass `genhmass` as the target variable, what happens to the regression? Why will this not work? "]}, {"cell_type": "markdown", "id": "2adf5fce", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='exercises_16_3'></a>     \n", "\n", "| [Top](#section_16_0) | [Restart Section](#section_16_3) | [Next Section](#section_16_4) |\n"]}, {"cell_type": "markdown", "id": "103a9b46", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.3.1</span>\n", "\n", "The Higgs boson has a mass of 125 GeV, which means that for a Higgs at rest, the momentum of each tau will be 125/2 = 62.5 GeV. This also means that the visible components of the tau will be less than 62.5. For instance, if the visible components take up half the tau momentum, then the visible momentum is 31.25 GeV. The correction to account for the neutrino will be large or small depending on whether the tau has low  or high momentum, respectively. In this problem, we will verify that the neural network has learned this momentum dependence by looking at the correction that it makes for a visible tau momentum of 20 GeV (low), compared to a visible tau momentum of 200 GeV(high). \n", "\n", "Now, you might ask, how do you get 200 GeV taus as decay products from the Higgs? This can happen because the Higgs is typically **not** produced at rest. Instead, it usually has some non-zero momentum that is a result of the production process inside the proton collision. Basically, other quarks in the proton recoil off the Higgs, giving it momentum. As a result, we can get higher momentum Higgs bosons decaying into higher momentum taus. However, the rate of these types of events rapidly drops off with increasing momentum, with the highest typical Higgs momentum being 50 GeV.\n", "\n", "Compute the relative NN correction for a 20 GeV input, compared to the correction for a 200 GeV input (use the state of the network after training for 500 epochs, i.e., after having run `L16.3-runcell07`). Express these results in terms of a fractional correction (the ratio of final over initial momenta), in order to compare the relative scale properly. Report your answer as a list of two numbers with precision 1e-3: `[fractional correction for 20 GeV input , fractional correction for 200 GeV input]`\n"]}, {"cell_type": "code", "execution_count": null, "id": "a14c4715", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>EXERCISE: L16.3.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "pVal=torch.tensor([20]).float()\n", "corr_20GeV = net(pVal).data.numpy()\n", "print(\"Fraction Correction [20 GeV Input]: \", #YOUR CODE HERE)\n", "\n", "pVal=torch.tensor([200]).float()\n", "corr_200GeV = net(pVal).data.numpy()\n", "print(\"Fraction Correction [200 GeV Input]: \", #YOUR CODE HERE)"]}, {"cell_type": "markdown", "id": "0ce2a772", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.3.2</span>\n", "\n", "In this section, we showed that we could reconstruct the Higgs mass based on tau momenta. We then applied the NN to the reconstructed (observed) tau momenta and used the new NN-corrected momenta to calculate the Higgs mass, which exhibited a nice peak at the expected value.\n", "\n", "Before moving on, let's consider sources of bias in this analysis. Select all options below that characterize a source of bias:\n", "\n", "A) Our NN is potentially biased because the features are all corrleated.\\\n", "B) Our NN is potentially biased because we gave it one Higgs sample at mass of 125 and so it will assume all taus no matter what energy they have came from a Higgas with a mass of 125.\\\n", "C) Our NN is potentially biased because some variables have clear correlations with our regression target, when the tau has a specific energy, but not all energies.\\\n", "D) Our NN is potentially biased because our simulation is not as good at simulating taus with hadronic neutral pion decays.\n", "\n"]}, {"cell_type": "markdown", "id": "8d118b99", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_4'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.4 The Full Mass Regression</h2>     \n", "\n", "| [Top](#section_16_0) | [Previous Section](#section_16_3) | [Exercises](#exercises_16_4) | [Next Section](#section_16_5) |\n"]}, {"cell_type": "markdown", "id": "4720831c", "metadata": {"tags": ["learner", "md"]}, "source": ["*The material in this section is discussed in the video **<a href=\"https://courses.mitxonline.mit.edu/learn/course/course-v1:MITxT+8.S50.2x+1T2025/block-v1:MITxT+8.S50.2x+1T2025+type@sequential+block@seq_LS16/block-v1:MITxT+8.S50.2x+1T2025+type@vertical+block@vert_LS16_vid4\" target=\"_blank\">HERE</a>.** You are encouraged to watch that video and use this notebook concurrently.*"]}, {"cell_type": "markdown", "id": "903ff83e", "metadata": {"tags": ["learner", "md"]}, "source": ["<h3>Slides</h3>\n", "\n", "Run the code below to view the slides for this section, which are discussed in the related video. You can also open the slides in a separate window <a href=\"https://mitx-8s50.github.io/slides/L19/slides_L19_08.html\" target=\"_blank\">HERE</a>."]}, {"cell_type": "code", "execution_count": null, "id": "88015a5d", "metadata": {"tags": ["learner", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-slides\n", "\n", "from IPython.display import IFrame\n", "IFrame(src='https://mitx-8s50.github.io/slides/L19/slides_L19_08.html', width=970, height=550)"]}, {"cell_type": "markdown", "id": "26ec5257", "metadata": {"tags": ["learner", "md", "lect_04"]}, "source": ["<h3>Overview</h3>\n", "\n", "Ok, now that we have a setup that seems to be working, the final challenge is to use more of the information in our events to do the best possible reconstruction. In this case, what we will do is use as many different variables as we possibly can to predict the neutrino properties. \n", "\n", "The way we are going to do this is by digging into the dataset a little bit. A tau is reconstructed in the same way a jet is reconstructed at the LHC. We first look for a clump of particles that have a large amount of energy in a cone. In this case, we also look for what we called an isolated clump. Namely, we see a few particles in the cone that have very large momentum and, when we remove these particles, the rest of the momentum in the cone is small (typically less that 10 percent of the total momentum). This final collection of particles gives us a tau. \n", "\n", "Interestingly, the decay components of the tau in the cone can tell us information about the energy of the neutrino. The reason for this comes from the properties of the tau decay. A tau has a fixed mass at 1.75 GeV. As a consequence, the invariant mass of the all of the particles within a cone should sum up to this value. In reality, we can't get this perfectly because our measurements are not perfect. However, we know that \n", "\n", "$$\n", "\\begin{eqnarray}\n", "\\left(\\sum_{i} \\vec{p_{i}}\\right)^2 & \\approx & m^{2}_{\\tau} \\\\\n", "\\left(\\sum_{i\\neq{\\rm neutrino}} \\vec{p_{i}} + \\sum_{i={\\rm neutrino}} \\vec{p_{i}} \\right)^2 & \\approx & m^{2}_{\\tau} \\\\\n", "\\end{eqnarray}\n", "$$\n", "\n", "Or, in other words, the neutrino position and direction are constrained by the other particles in the tau decay. We don't need to expand this equation and solve it, we just need to know there is a correlation in order to use the tau decay components to solve for the neutrino properties.  \n", "\n", "This all means that we can build a neural network with the particles as input to predict the neutrinos. Let's go ahead and look at a few taus, and then we build our network and train it.  "]}, {"cell_type": "code", "execution_count": null, "id": "daee5201", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell01\n", "\n", "#Let's compute the mass on the fly\n", "def makedataset(iMask,iPart=\"part1\"):\n", "    varlist=[iPart+\"pt1\",iPart+\"eta1\",iPart+\"phi1\",iPart+\"id1\",iPart+\"pt2\",iPart+\"eta2\",iPart+\"phi2\",iPart+\"id2\",iPart+\"pt3\",iPart+\"eta3\",iPart+\"phi3\",iPart+\"id3\",iPart+\"pt4\",iPart+\"eta4\",iPart+\"phi4\",iPart+\"id4\",iPart+\"pt5\",iPart+\"eta5\",iPart+\"phi5\",iPart+\"id5\"]\n", "    arr=0\n", "    idx=0\n", "    for x in varlist:\n", "        pArr=reg[x].array(library=\"np\")[iMask]\n", "        if idx == 0: \n", "            arr = pArr\n", "            idx = idx + 1\n", "        else:\n", "            arr=np.vstack((arr,pArr))\n", "    arr = arr.T\n", "    return arr\n", "\n", "mask1=(reg[\"genpt1\"].array(library=\"np\") > 0)\n", "mask2=(reg[\"recopt1\"].array(library=\"np\") > 0)\n", "mask3=(reg[\"genpt2\"].array(library=\"np\") > 0)\n", "mask4=(reg[\"recopt2\"].array(library=\"np\") > 0)\n", "mask = np.logical_and.reduce([mask1,mask2,mask3,mask4])\n", "x=torch.from_numpy(makedataset(mask))\n", "yb=torch.from_numpy(reg[\"recopt1\"].array(library=\"np\")[mask].reshape(len(reg[\"recopt1\"].array(library=\"np\")[mask]),1))\n", "y=torch.from_numpy(reg[\"genpt1\"].array(library=\"np\")[mask].reshape(len(reg[\"genpt1\"].array(library=\"np\")[mask]),1))\n", "ratio=torch.div(y,yb)\n", "y=ratio\n", "x,y = Variable(x),Variable(y)\n", "torch_dataset = Data.TensorDataset(x, y)\n", "#print(x)"]}, {"cell_type": "markdown", "id": "beaebe2a", "metadata": {"tags": ["learner", "md", "lect_04"]}, "source": ["<h3>Plotting Events in Position Space</h3>\n", "\n", "To get a feel for how these corrections work, we can visualize a few events by making some displays using the angular coordinates eta ($\\eta$) and phi ($\\phi$) to look at the directions of the top 5 particles. Recall that $\\phi$ is the azimuthal angle around the beam direction and $\\eta$ is the relativistic equivalent of $\\theta$, the polar angle from the beam direction. We will also make the size of the circles proportional to their energies. Each display will also show its correction factor (the ratio of the reconstructed over true tau momenta)."]}, {"cell_type": "code", "execution_count": null, "id": "af81fadc", "metadata": {"scrolled": false, "tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell02\n", "\n", "ds=makedataset(mask)\n", "colors = ['g','r','b','y','orange']\n", "for i0 in range(5):\n", "    for ipart in range(5):\n", "        plt.scatter(ds[i0,4*ipart+1], ds[i0,4*ipart+2], s=ds[i0,4*ipart]*5000/yb[i0], c=colors[ipart], alpha=0.5)\n", "    plt.xlim(-0.5,0.5)\n", "    plt.ylim(-0.5,0.5)\n", "    plt.xlabel(\"$\\eta$\")\n", "    plt.ylabel(\"$\\phi$\")\n", "    plt.text(-0.3,0.4,\"Correction Factor \"+str(ratio[i0].numpy()[0]))\n", "    plt.show()"]}, {"cell_type": "markdown", "id": "1d0ede7e", "metadata": {"tags": ["learner", "md", "lect_04"]}, "source": ["Broadly speaking, what we observe is that events that are spread out tend to have larger corrections than events that are not. This means we can use this shape to actually predict the correction. There is some more subtlety coming from the fact that taus decay into different particle types. That will allow us to understand the type of decay that occurred. For example, tau decays that have a lot of photons, will have a different energy distribution than taus that don't have any photons. \n", "\n", "Alright, let's go ahead and train our network. "]}, {"cell_type": "code", "execution_count": null, "id": "2fe551cc", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell03\n", "\n", "#now let's see if we can improve this with something more complicated\n", "torch.manual_seed(1)    # reproducible\n", "net = torch.nn.Sequential(\n", "        torch.nn.Linear(20, 200),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(200, 200),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(200, 50),\n", "        torch.nn.LeakyReLU(),\n", "        torch.nn.Linear(50, 1),\n", "    )\n", "optimizer = torch.optim.Adam(net.parameters(), lr=0.01)\n", "loss_func = torch.nn.MSELoss()"]}, {"cell_type": "code", "execution_count": null, "id": "92bd4e77", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell04\n", "\n", "p1=torch.from_numpy(makedataset(mask,\"part1\"))\n", "p2=torch.from_numpy(makedataset(mask,\"part2\"))\n", "\n", "def masscomputeNN(iC1,iC2,iVec1,iVec2):\n", "    corr1 = iC1*iVec1[0]\n", "    corr2 = iC2*iVec2[0]\n", "    tau_1 = Momentum4.m_eta_phi_pt(iVec1[3], iVec1[1], iVec1[2], corr1)\n", "    tau_2 = Momentum4.m_eta_phi_pt(iVec2[3], iVec2[1], iVec2[2], corr2)\n", "    return (tau_1+tau_2).m\n", "\n", "#now let's compute the corrected mass on the data set\n", "def hmass(masscomputeNN):\n", "    corr1=net(p1)\n", "    corr2=net(p2)\n", "    varlist=[\"recopt1\",\"recoeta1\",\"recophi1\",\"recomass1\",\"recopt2\",\"recoeta2\",\"recophi2\",\"recomass2\"]\n", "    arr=np.vstack((corr1.data.numpy().T,corr2.data.numpy().T))\n", "    for x in varlist:\n", "        pArr=reg[x].array(library=\"np\")[mask]\n", "        arr=np.vstack((arr,pArr))\n", "    arr = arr.T\n", "    massc = lambda iarr: masscomputeNN(iarr[0],iarr[1],iarr[2:6],iarr[6:10]) \n", "    hmasses = np.array([massc(p) for p in arr])\n", "    return hmasses\n", "\n", "#now update to add history\n", "history_lr = {'loss':[], 'val_loss':[]}\n", "def train(x,y,net,loss_func,opt,nepochs,ymin,ymax):\n", "    images = []\n", "    fig, ax = plt.subplots(figsize=(12,7))\n", "    for epoch in range(nepochs):\n", "        if epoch % 50 == 0: \n", "            print(\"epoch:\",epoch)\n", "        prediction = net(x)\n", "        loss = loss_func(prediction, y) \n", "        opt.zero_grad()\n", "        loss.backward() \n", "        optimizer.step()\n", "        with torch.no_grad():#disable updating gradient\n", "            if epoch % 50 == 0:\n", "                print('[%d] loss: %.4f ' % (epoch + 1, loss ))\n", "            history_lr['loss'].append(loss)\n", "        if epoch % 5 == 0:\n", "            makePlot(x,y,prediction,ax,fig,images,epoch,loss,ymin,ymax)\n", "    return images\n", "\n", "rawmass=hmass(masscomputeNN)\n", "plt.hist(rawmass,bins=50,range=(0,500),color='blue',alpha=0.5,label=\"Adding the NN\")\n", "plt.xlabel(\"mass(GeV)\")\n", "plt.ylabel(\"N$_{events}$\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "id": "f6f36bae", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell05\n", "\n", "#NOTE: if training does not complete due to timeout in Colab,\n", "#reduce the number of epochs to 250 and run this cell twice,\n", "#or reduce to 125 and run this cell four times,\n", "#for a total of 500 training epochs\n", "images=train(x,y,net,loss_func,optimizer,500,0,1)\n", "torch.save(net.state_dict(), 'data/L16/tau_reg_fullpart.pt')\n", "imageio.mimsave('data/L16/full_reg2.gif', images, fps=12)\n", "Image(open('data/L16/full_reg2.gif','rb').read())"]}, {"cell_type": "markdown", "id": "21251f2b", "metadata": {"tags": ["learner", "md", "lect_04"]}, "source": ["Let's look at how the loss evolves as a function of epoch."]}, {"cell_type": "code", "execution_count": null, "id": "9e8c5def", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell06\n", "\n", "tmploss=[]\n", "for i0 in range(len(history_lr['loss'])):\n", "    tmploss.append(history_lr['loss'][i0].detach().numpy())\n", "plt.semilogy(tmploss, label='loss')\n", "plt.legend(loc=\"upper right\")\n", "plt.xlabel('epoch')\n", "plt.ylabel('loss (binary crossentropy)')\n", "plt.show()"]}, {"cell_type": "markdown", "id": "ef431e4f", "metadata": {"tags": ["learner", "md", "lect_04"]}, "source": ["And, finally, we can compare the generated, reconstructed, and neural-network corrected mass distributions."]}, {"cell_type": "code", "execution_count": null, "id": "2a4df160", "metadata": {"tags": ["learner", "py", "lect_04", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.4-runcell07\n", "\n", "def plotcorr(iVar,iNN,iMin,iMax,iColor,iLabel,iCorr=True): \n", "    corr1=iNN(p1)\n", "    data=reg[iVar].array(library=\"np\")[mask]\n", "    if iCorr:\n", "        data=data*corr1.data.numpy().T\n", "    counts, binEdges = np.histogram(data,bins=50,range=(iMin,iMax),density=False)\n", "    binCenters = (binEdges[1:]+binEdges[:-1])*.5\n", "    err = np.sqrt(counts)\n", "    plt.errorbar(binCenters, counts, yerr=err,fmt=\"o\",c=iColor, ms=3,label=iLabel)\n", "    \n", "plotcorr(\"genpt1\" ,net,0,200,\"black\",\"gen\",False)\n", "plotcorr(\"recopt1\",net,0,200,\"red\",\"reco\",False)\n", "plotcorr(\"recopt1\",net,0,200,\"blue\",\"corrected\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "fdbe6b95", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='exercises_16_4'></a>     \n", "\n", "| [Top](#section_16_0) | [Restart Section](#section_16_4) | [Next Section](#section_16_5) |\n"]}, {"cell_type": "markdown", "id": "2f0dda41", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.4.1</span>\n", "\n", "Let's see how well this works by looking at the correction for the second particle, which was not included in the training. Run the code cell below in order to compare the NN correction with the true correction (use the state of the network after training for 500 epochs, i.e., after having run `L16.4-runcell05`). Report your answer as a list of two numbers with precision 1e-2: `[NN correction, true correction]'"]}, {"cell_type": "code", "execution_count": null, "id": "a5473421", "metadata": {"tags": ["py", "draft", "learner_chopped"]}, "outputs": [], "source": ["#>>>EXERCISE: L16.4.1\n", "\n", "def plotcorrp2(iTrueVar,iRecoVar,iNN,iMin,iMax,iColor,iLabel,iCorr=True): \n", "    corr1=iNN(p2).detach().numpy()    \n", "    true=reg[iTrueVar].array(library=\"np\")[mask]\n", "    reco=reg[iRecoVar].array(library=\"np\")[mask]\n", "    corr2=true/reco\n", "    counts, binEdges = np.histogram(corr1,bins=50,range=(iMin,iMax),density=False)\n", "    binCenters = (binEdges[1:]+binEdges[:-1])*.5\n", "    err = np.sqrt(counts)\n", "    plt.errorbar(binCenters, counts, yerr=err,fmt=\"o\",c=iColor, ms=3,label=\"NN\")\n", "    counts, binEdges = np.histogram(corr2,bins=50,range=(iMin,iMax),density=False)\n", "    plt.errorbar(binCenters, counts, yerr=err,fmt=\"o\",c=\"Black\", ms=3,label=\"True\")\n", "    print(\"NN   Mean : \",np.mean(corr1),\"\\t RMS: \",corr1.std())\n", "    print(\"True Mean : \",np.mean(corr2),\"\\t RMS: \",corr2.std())\n", "    \n", "plotcorrp2(\"genpt2\",\"recopt2\" ,net,0,4,\"green\",\"gen\",False)\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "cfbf9b64", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='section_16_5'></a>\n", "<hr style=\"height: 1px;\">\n", "\n", "## <h2 style=\"border:1px; border-style:solid; padding: 0.25em; color: #FFFFFF; background-color: #90409C\">L16.5 Tuning the NN Architecture</h2>     \n", "\n", "| [Top](#section_16_0) | [Previous Section](#section_16_4) | [Exercises](#exercises_16_5) |\n"]}, {"cell_type": "markdown", "id": "a113055f", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["<h3>Overview</h3>\n", "\n", "*Note: There is no corresponding video for this section.*\n", "\n", "Now, one last thing that we can consider is to use a different type of architecture for the neutral network. This is where deep learning has become very valuable over the last 10 years. By making it such that networks take only a subset of inputs at a time, we can effectively organize the training of the data, and guide the deep learning architecture to a better overall minimum. \n", "\n", "For the tau regression, the ideal tool for this is a recurrent neural network or an <a href=\"https://en.wikipedia.org/wiki/Recurrent_neural_network\" target=\"_blank\">RNN</a>. These networks are designed to take in data piece by piece to build up a prediction. They are often used on time series datasets to process data and infer the likelihood of an event or to apply a correction. Perhaps their most common use has been in language translation to translate word by word, or phrase by phrase. \n", "\n", "In recent times, high energy physicists have taken to using RNNs to process data by feeding the RNN individual particles one at a time. The nice thing is that since this takes in a series of particles, we could, in principle, give it a variable number of particles. \n", "\n", "In this example, we will use a long-short term memory (LSTM) network. Here is the pytorch for how it's defined. "]}, {"cell_type": "code", "execution_count": null, "id": "d38db36b", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell01\n", "\n", "torch.manual_seed(1)    # reproducible\n", "\n", "class LSTM(nn.Module):\n", "\n", "    def __init__(self, num_classes, input_size, hidden_size, num_layers):\n", "        super(LSTM, self).__init__()\n", "        \n", "        self.num_classes = num_classes\n", "        self.num_layers = num_layers\n", "        self.input_size = input_size\n", "        self.hidden_size = hidden_size\n", "        self.lstm = nn.LSTM(input_size=input_size, hidden_size=hidden_size,\n", "                            num_layers=num_layers, batch_first=True)\n", "        \n", "        self.fc1 = nn.Linear(hidden_size, 20)\n", "        self.fc2 = nn.Linear(20, num_classes)\n", "\n", "    def forward(self, x):\n", "        h_0 = Variable(torch.zeros(\n", "            self.num_layers, x.size(0), self.hidden_size))\n", "        \n", "        c_0 = Variable(torch.zeros(\n", "            self.num_layers, x.size(0), self.hidden_size))\n", "        \n", "        # Propagate input through LSTM\n", "        ula, (h_out, _) = self.lstm(x, (h_0, c_0))\n", "        \n", "        h_out = h_out.view(-1, self.hidden_size)\n", "        \n", "        out = self.fc1(h_out)\n", "        out = F.relu(out)\n", "        out = self.fc2(out)\n", "        return out\n", "\n", "input_size = 4 # take in 4 vectors\n", "hidden_size = 128 # hidden layers\n", "num_layers = 1 # output layers\n", "num_classes = 1 # output values (just 1 the correction)\n", "lstm = LSTM(num_classes, input_size, hidden_size, num_layers)\n", "criterion = torch.nn.MSELoss()    # mean-squared error for regression\n", "optimizer = torch.optim.Adam(lstm.parameters(), lr=0.01)\n", "lstm.train()"]}, {"cell_type": "markdown", "id": "d0c15590", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Now, to process the Higgs-to-tau data, we need to format it slightly differently than was done before, where we just took all of our particles and made them into one giant vector. What we are going to do here is to feed one particle at a time into the network. \n", "\n", "The network will then run the its inference and send an output (the state) to itself along with an output out of the network. We can keep iterating by giving more particles to the network and also reading in the input state from the previous network inference. \n", "\n", "To build this input dataset, we will take our previous collection of 5 particles, and reshape it to correspond to \n", "\n", " * 5x4 vectors per event\n", " \n", "This is easily done using numpy commands"]}, {"cell_type": "code", "execution_count": null, "id": "aacf38df", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell02\n", "\n", "def makedatasetrnn(iMask,iPart=\"part1\"):\n", "    arr=makedataset(iMask,iPart)\n", "    return arr.reshape(len(arr),5,4)\n", "\n", "x=torch.from_numpy(makedatasetrnn(mask))"]}, {"cell_type": "markdown", "id": "1f29946c", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Now, just like we did before, we can make an updated dataset for each of the tau legs in the Higgs decays. We can then run our LSTM on each of the legs and compute the mass on the fly. Moreover, the way pytorch is setup, it will automatically run the neural network on all of the particles. Thanks Pytorch! \n", "\n", "Let's code up the mass calculation and try it on our untrained network. "]}, {"cell_type": "code", "execution_count": null, "id": "c9a53ad1", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell03\n", "\n", "p1=torch.from_numpy(makedatasetrnn(mask,\"part1\"))\n", "p2=torch.from_numpy(makedatasetrnn(mask,\"part2\"))\n", "\n", "#now let's compute the corrected mass on the data set\n", "def hmass(masscomputeNN):\n", "    mask=(reg[\"recohmass\"].array(library=\"np\") > 0)\n", "    corr1=lstm(p1)\n", "    corr2=lstm(p2)\n", "    varlist=[\"recopt1\",\"recoeta1\",\"recophi1\",\"recomass1\",\"recopt2\",\"recoeta2\",\"recophi2\",\"recomass2\"]\n", "    arr=np.vstack((corr1.data.numpy().T,corr2.data.numpy().T))\n", "    for x in varlist:\n", "        pArr=reg[x].array(library=\"np\")[mask]\n", "        arr=np.vstack((arr,pArr))\n", "    arr = arr.T\n", "    massc = lambda iarr: masscomputeNN(iarr[0],iarr[1],iarr[2:6],iarr[6:10]) \n", "    hmasses = np.array([massc(p) for p in arr])\n", "    return hmasses\n", "\n", "outmass=hmass(masscomputeNN)\n", "plt.hist(outmass,bins=40,range=(0,250),color='blue')\n", "plt.xlabel(\"mass\")\n", "plt.ylabel(\"N$_{events}$\")\n", "plt.show()"]}, {"cell_type": "markdown", "id": "2da99d93", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Alright, now we can go ahead and train our recursive neural network. Note that this network is a bit smaller than our previous network.\n", "\n", "One of the reasons that deep learning has become so mainstream is because we now have fancy architectures that only take in an organized set of inputs. This allows us to be much more expressive in the style of deep learning that we do! "]}, {"cell_type": "code", "execution_count": null, "id": "b48eeff9", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell04\n", "\n", "#NOTE: if training does not complete due to timeout in Colab,\n", "#reduce the number of epochs to 250 and run this cell twice,\n", "#or reduce to 125 and run this cell four times,\n", "#for a total of 500 training epochs\n", "history_lr = {'loss':[], 'val_loss':[]}\n", "images=train(x,y,lstm,criterion,optimizer,500,0,1)\n", "torch.save(lstm.state_dict(), 'data/L16/tau_reg_lstm.pt')\n", "imageio.mimsave('data/L16/reg_lstm.gif', images, fps=12)\n", "Image(open('data/L16/reg_lstm.gif','rb').read())\n"]}, {"cell_type": "markdown", "id": "8e6d5710", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Now, let's look at the performance of our loss over each of the steps so we can see how the convergence looks and how the fit is doing. Sometimes, we get jumps in our losses, which occur because our minimizer will sometimes go on large excursions from the minimum. "]}, {"cell_type": "code", "execution_count": null, "id": "99fdcd0b", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell05\n", "\n", "tmploss=[]\n", "\n", "for i0 in range(len(history_lr['loss'])):\n", "    tmploss.append(history_lr['loss'][i0].detach().numpy())\n", "plt.semilogy(tmploss, label='loss')\n", "plt.legend(loc=\"upper right\")\n", "plt.xlabel('epoch')\n", "plt.ylabel('loss (binary crossentropy)')\n", "plt.show()"]}, {"cell_type": "markdown", "id": "07d19689", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Now, we can check our dataset by looking at the reconstructed momentum like we did before. "]}, {"cell_type": "code", "execution_count": null, "id": "6fee7b6a", "metadata": {"tags": ["learner", "py", "lect_05", "learner_chopped"]}, "outputs": [], "source": ["#>>>RUN: L16.5-runcell06\n", "\n", "def plotcorr(iVar,iNN,iMin,iMax,iColor,iLabel,iCorr=True):\n", "    corr1=iNN(p1)\n", "    data=reg[iVar].array(library=\"np\")[mask]\n", "    if iCorr:\n", "        data=data*corr1.data.numpy().T\n", "    counts, binEdges = np.histogram(data,bins=50,range=(iMin,iMax),density=False)\n", "    binCenters = (binEdges[1:]+binEdges[:-1])*.5\n", "    err = np.sqrt(counts)\n", "    plt.errorbar(binCenters, counts, yerr=err,fmt=\"o\",c=iColor, ms=3,label=iLabel)\n", "    \n", "plotcorr(\"genpt1\" ,lstm,0,200,\"black\",\"gen\",False)\n", "plotcorr(\"recopt1\",lstm,0,200,\"red\",\"reco\",False)\n", "plotcorr(\"recopt1\",lstm,0,200,\"blue\",\"corrected\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "id": "778a9ef3", "metadata": {"tags": ["learner", "md", "lect_05"]}, "source": ["Our mass definitely looks a bit peakier than it did with the previous algorithm."]}, {"cell_type": "markdown", "id": "fe8da95a", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["<a name='exercises_16_5'></a>   \n", "\n", "| [Top](#section_16_0) | [Restart Section](#section_16_5) |\n"]}, {"cell_type": "markdown", "id": "df9f7e93", "metadata": {"tags": ["learner", "md", "learner_chopped"]}, "source": ["### <span style=\"border:3px; border-style:solid; padding: 0.15em; border-color: #90409C; color: #90409C;\">Exercise 16.5.1</span>\n", "\n", "Complete the code cell below to compute the correlation coefficient with the LSTM that we have made with the gen momentum. How does this compare to our initial correlation before using the neural network? Use the state of the network after training for 1000 epochs, i.e., after having run `L16.5-runcell04`.\n", "\n", "Enter your answer as a list of numbers with precision 1e-2: `[corr-original, corr-NN]`"]}, {"cell_type": "code", "execution_count": null, "id": "5c040ac1", "metadata": {"tags": ["draft", "py", "learner_chopped"]}, "outputs": [], "source": ["#>>>EXERCISE: L16.5.1\n", "# Use this cell for drafting your solution (if desired),\n", "# then enter your solution in the interactive problem online to be graded.\n", "\n", "pred=lstm(p1)\n", "ytmp=y.detach().numpy().flatten()*(yb.detach().numpy().flatten())\n", "ptmp=pred.detach().numpy().flatten()*(yb.detach().numpy().flatten())\n", "\n", "print(\"Pre Correlation:\", #YOUR CODE HERE)\n", "print(\"NN Correlation:\", #YOUR CODE HERE)"]}], "metadata": {"celltoolbar": "Tags", "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.9.9"}}, "nbformat": 4, "nbformat_minor": 5}