8S50x_L16_runcells.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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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": "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}