lecture04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dd28da6c-e768-47a1-978f-d5631e1b47f2",
   "metadata": {},
   "source": [
    "# Lecture 4 - Jet images and transfer learning with CNNs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b97e80e-4a9d-4fef-ad66-cc79fead059a",
   "metadata": {},
   "source": [
    "> A look at converting jets into images and classifying them with convolutional neural networks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae64594c-550f-4227-ad36-305939a29778",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "\n",
    "* Understand how to convert jets into images\n",
    "* Know what transfer learning is and it's advantages compared to training a neural network from scratch\n",
    "* Understand the main steps needed to fine-tune a convolutional neural network"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b5c1de1-edd1-4e79-85d5-a33d19411a65",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* Chapter 1 of [_Deep Learning for Coders with fastai & PyTorch_](https://github.com/fastai/fastbook) by Jeremy Howard and Sylvain Gugger.\n",
    "* [_The Machine Learning Landscape of Top Taggers_](https://arxiv.org/abs/1902.09914) by G. Kasieczka et al.\n",
    "* [Jet-Images -- Deep Learning Edition](https://arxiv.org/abs/1511.05190) by L. de Oliviera et al."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db23e4a5",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "85b47ca2-e93a-4f8a-81f9-85a50deb061b",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ed1c15-faea-42dd-a0fa-05d5305c010f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment and run this cell if using Colab, Kaggle etc\n",
    "# %pip install fastai==2.6.0 datasets energyflow huggingface_hub"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a24705d7-e5a4-4adf-bb16-aca47aa3c48e",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "e5387606-9d8e-40d5-9826-e2cf20620766",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from datasets import load_dataset\n",
    "from energyflow.utils import (\n",
    "    center_ptyphims,\n",
    "    phi_fix,\n",
    "    pixelate,\n",
    "    ptyphims_from_p4s,\n",
    "    reflect_ptyphims,\n",
    "    rotate_ptyphims,\n",
    ")\n",
    "from fastai.vision.all import *\n",
    "from huggingface_hub import from_pretrained_fastai\n",
    "from sklearn.metrics import accuracy_score, auc, roc_auc_score, roc_curve\n",
    "from torch.utils.data import Dataset, TensorDataset\n",
    "from torchvision.transforms import ToTensor\n",
    "from scipy.interpolate import interp1d\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7af6929a-9e1d-493f-9c84-281c706d9467",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datasets\n",
    "\n",
    "# Suppress logs to keep things tidy\n",
    "datasets.logging.set_verbosity_error()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11426624-9951-4ae9-81cc-1a8035453120",
   "metadata": {},
   "source": [
    "## Creating jet images"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dce7dc6c-bd79-4944-b299-0fb0dc57d729",
   "metadata": {},
   "source": [
    "As we saw in lectures 1 and 2, there are various ways to represent jet data for neural networks. So far, we've focused on $N$-subjettiness features $\\tau_N^{\\beta}$, which represent jets in a _tabular_ format. Today, we'll examine another popular representation that treats _jets as 2D images_. An example from the top tagging review is shown below:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d630810d-5da3-4f35-bc4a-46c0be9ea93d",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img alt=\"jet-images\" caption=\"Jet images\" src=\"images/jet-images.png\" id=\"jet-images\" width=800/>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1964576-21e1-4e5a-a828-c2dd3a0c367f",
   "metadata": {},
   "source": [
    "But why consider images in the first place? The motivation here is partly historical: computer vision is where the deep learning revolution started in 2011, with [DanNet](https://people.idsia.ch/~juergen/DanNet-triggers-deep-CNN-revolution-2011.html) and [AlexNet](https://en.wikipedia.org/wiki/AlexNet) smashing the state-of-the-art on popular benchmarks. Both approaches were based on a type of architecture called a [_convolutional neural network_](https://en.wikipedia.org/wiki/Convolutional_neural_network) (CNN), which is especially well suited for analysing images. In the next lecture we'll examine how CNNs work in detail, so today we'll focus on creating jet images and training a CNN with the high-level API of fastai."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ab99be-3859-4de2-83b5-9dbc7fa0b18b",
   "metadata": {},
   "source": [
    "To get started, let's take a look at how we can construct images from the 4-vectors associated with the consituents in a jet. We'll use the raw events from the [_Top Quark Tagging_ dataset](https://huggingface.co/datasets/dl4phys/top_tagging) in lecture 1, so let's download it and grab a sample of 10,000 events to play with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b81e9af9-b51f-4a8a-95ee-a0d9d5aff9c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>E_0</th>\n",
       "      <th>PX_0</th>\n",
       "      <th>PY_0</th>\n",
       "      <th>PZ_0</th>\n",
       "      <th>E_1</th>\n",
       "      <th>PX_1</th>\n",
       "      <th>PY_1</th>\n",
       "      <th>PZ_1</th>\n",
       "      <th>E_2</th>\n",
       "      <th>PX_2</th>\n",
       "      <th>...</th>\n",
       "      <th>E_199</th>\n",
       "      <th>PX_199</th>\n",
       "      <th>PY_199</th>\n",
       "      <th>PZ_199</th>\n",
       "      <th>truthE</th>\n",
       "      <th>truthPX</th>\n",
       "      <th>truthPY</th>\n",
       "      <th>truthPZ</th>\n",
       "      <th>ttv</th>\n",
       "      <th>is_signal_new</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>130.172562</td>\n",
       "      <td>114.353065</td>\n",
       "      <td>-59.476486</td>\n",
       "      <td>18.188450</td>\n",
       "      <td>58.261734</td>\n",
       "      <td>56.383465</td>\n",
       "      <td>-6.442986</td>\n",
       "      <td>-13.184177</td>\n",
       "      <td>43.596382</td>\n",
       "      <td>37.914547</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>496.599976</td>\n",
       "      <td>418.331940</td>\n",
       "      <td>-202.994308</td>\n",
       "      <td>43.685955</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>144.244232</td>\n",
       "      <td>44.394894</td>\n",
       "      <td>71.800972</td>\n",
       "      <td>-116.962021</td>\n",
       "      <td>176.547729</td>\n",
       "      <td>61.791088</td>\n",
       "      <td>51.010948</td>\n",
       "      <td>-157.317657</td>\n",
       "      <td>134.554794</td>\n",
       "      <td>47.264194</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1061.465576</td>\n",
       "      <td>240.111099</td>\n",
       "      <td>465.869263</td>\n",
       "      <td>-906.836731</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>164.264938</td>\n",
       "      <td>29.799849</td>\n",
       "      <td>-155.011093</td>\n",
       "      <td>45.458752</td>\n",
       "      <td>154.540726</td>\n",
       "      <td>34.444710</td>\n",
       "      <td>-140.643768</td>\n",
       "      <td>53.997475</td>\n",
       "      <td>102.357109</td>\n",
       "      <td>19.766624</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>190.001175</td>\n",
       "      <td>-59.107372</td>\n",
       "      <td>-180.212479</td>\n",
       "      <td>11.411892</td>\n",
       "      <td>36.530270</td>\n",
       "      <td>-16.578011</td>\n",
       "      <td>-32.546757</td>\n",
       "      <td>0.582238</td>\n",
       "      <td>28.827248</td>\n",
       "      <td>-10.152600</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>197.350937</td>\n",
       "      <td>-190.924622</td>\n",
       "      <td>-28.691614</td>\n",
       "      <td>-40.889717</td>\n",
       "      <td>65.894730</td>\n",
       "      <td>-63.749012</td>\n",
       "      <td>-9.580021</td>\n",
       "      <td>-13.652922</td>\n",
       "      <td>53.116070</td>\n",
       "      <td>-52.294987</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 806 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          E_0        PX_0        PY_0        PZ_0         E_1       PX_1  \\\n",
       "0  130.172562  114.353065  -59.476486   18.188450   58.261734  56.383465   \n",
       "1  144.244232   44.394894   71.800972 -116.962021  176.547729  61.791088   \n",
       "2  164.264938   29.799849 -155.011093   45.458752  154.540726  34.444710   \n",
       "3  190.001175  -59.107372 -180.212479   11.411892   36.530270 -16.578011   \n",
       "4  197.350937 -190.924622  -28.691614  -40.889717   65.894730 -63.749012   \n",
       "\n",
       "         PY_1        PZ_1         E_2       PX_2  ...  E_199  PX_199  PY_199  \\\n",
       "0   -6.442986  -13.184177   43.596382  37.914547  ...    0.0     0.0     0.0   \n",
       "1   51.010948 -157.317657  134.554794  47.264194  ...    0.0     0.0     0.0   \n",
       "2 -140.643768   53.997475  102.357109  19.766624  ...    0.0     0.0     0.0   \n",
       "3  -32.546757    0.582238   28.827248 -10.152600  ...    0.0     0.0     0.0   \n",
       "4   -9.580021  -13.652922   53.116070 -52.294987  ...    0.0     0.0     0.0   \n",
       "\n",
       "   PZ_199       truthE     truthPX     truthPY     truthPZ  ttv  is_signal_new  \n",
       "0     0.0   496.599976  418.331940 -202.994308   43.685955    2              1  \n",
       "1     0.0  1061.465576  240.111099  465.869263 -906.836731    2              1  \n",
       "2     0.0     0.000000    0.000000    0.000000    0.000000    2              0  \n",
       "3     0.0     0.000000    0.000000    0.000000    0.000000    2              0  \n",
       "4     0.0     0.000000    0.000000    0.000000    0.000000    2              0  \n",
       "\n",
       "[5 rows x 806 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_events = load_dataset(\"dl4phys/top_tagging\", split=\"validation\")\n",
    "sample_df = raw_events.shuffle(seed=42).select(range(10_000)).to_pandas()\n",
    "sample_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "407e1e46-9e71-4751-a426-e088f6f082a6",
   "metadata": {},
   "source": [
    "Here we've used the `shuffle()` and `select()` methods to create a random sample, and then converted the result to a pandas `DataFrame`. It's convenient to reshape the events so that instead of having 800 columns, we have 200 columns, where each column groups the 4-vectors of a constituent in a single array. We can do this by casting our `DataFrame` to NumPy and applying the `reshape()` method as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac3cc2f1-e0b3-4893-a5f8-918a7d637de3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((10000, 200, 4), (10000,))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Reshape to size (num_jets, num_consituents, 4)\n",
    "events = sample_df.iloc[:, :800].values.reshape(-1, 200, 4)\n",
    "# Extract array of labels\n",
    "labels = sample_df[\"is_signal_new\"].values\n",
    "\n",
    "events.shape, labels.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa30dac-94bd-4872-8fb7-2ac0b1337c23",
   "metadata": {},
   "source": [
    "Now, each row in our `events` array corresponds to an array of shape `(num_constituents, 4)`. We can inspect one of these constituent 4-vectors as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "07eb663e-e6e4-40ef-a91e-cb2f9977575b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([130.17256 , 114.353065, -59.476486,  18.18845 ], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "events[0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "350009c1-2b7d-4793-909d-db60c46132b1",
   "metadata": {},
   "source": [
    "### Rotating to hadronic coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b19970d",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img alt=\"jet-images\" caption=\"Jet images\" src=\"images/jet-tagging.png\" id=\"jet-images\" width=600/>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3380fdd5-43d9-4858-a3f0-b87a03fac4a7",
   "metadata": {},
   "source": [
    "Each of the 4-vectors in `events` is currently stored in Cartesian coordinates $(E, p_x, p_y, p_z)$. However, it is convenient to rotate the basis to _hadronic coordinates_ $(p_T, y, \\phi, m)$, where:\n",
    "\n",
    "$$ p_T = \\sqrt{p_x^2 + p_y^2} \\,, \\quad y = \\mathrm{arctanh}\\, \\frac{p_z}{E} \\,, \\quad \\phi = \\arctan \\frac{p_y}{p_x}\\,, \\quad m = \\sqrt{E - p_x^2 - p_y^2 - p_z^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a79c5ba5-54b9-4182-9229-1acddd2eb694",
   "metadata": {},
   "source": [
    "Here $p_T$ denotes the transverse momentum, $y$ is the rapidity, $\\phi$ the azimuthal angle, and $m$ the mass. Although we could implement these formulas directly in NumPy, we'll use the [EnergyFlow](https://energyflow.network/) library instead, which provides many utility functions for precisely these cases. To convert our 4-vectors to hadronic coordinates, we use the `ptyphims_from_p4s()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4b40e710-da8a-4688-8505-28f56661c25b",
   "metadata": {},
   "outputs": [],
   "source": [
    "events_hadronic = ptyphims_from_p4s(events, phi_ref=\"hardest\", mass=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01882ffa-ad1b-4ffc-a629-542a4be8342f",
   "metadata": {},
   "source": [
    "Here, the `phi_ref` argument specifies which $\\phi$ value to use as reference within $\\pm \\pi$. The next step is to centre the collection of 4-vectors associated with each event. We can use the `center_ptyphims()` function to handle that for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "44834eb4-820a-4a03-aeeb-53d16bf5a056",
   "metadata": {},
   "outputs": [],
   "source": [
    "events_centered = [center_ptyphims(event) for event in events_hadronic]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d729a79-2aca-4ea2-9d0e-94512e6f2d55",
   "metadata": {},
   "source": [
    "The final step is to reflect and rotate our events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6a3e9d2c-c806-4226-9448-62a2ebe060d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "events_reflected_and_rotated = [\n",
    "    reflect_ptyphims(rotate_ptyphims(event)) for event in events_centered\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef6dc012-979b-4392-8cb9-f543a05eb0cf",
   "metadata": {},
   "source": [
    "The final step is to create a jet image from an array of 4-vectors.  Here we can use the `pixelate()` function from EnergyFlow and we'll follow the precription from the top tagging review to create images of $40\\times 40$ pixels. The end result is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "550f7da2-0bba-415e-b914-bdcfffa87081",
   "metadata": {},
   "outputs": [],
   "source": [
    "images = [\n",
    "    pixelate(\n",
    "        event,\n",
    "        npix=40,\n",
    "        img_width=0.8,\n",
    "        nb_chan=1,\n",
    "        norm=False,\n",
    "        charged_counts_only=False,\n",
    "    )\n",
    "    for event in events_reflected_and_rotated\n",
    "]\n",
    "\n",
    "images = np.array(images).reshape(len(images), 40, 40)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8630bcdb-af5f-4223-95ae-b89470c85048",
   "metadata": {},
   "source": [
    "Now that we have our images, let's split them into those corresponding to the top-quarks and QCD background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "69c93bbc-736a-4d14-8d59-a1034bdcf975",
   "metadata": {},
   "outputs": [],
   "source": [
    "top_images = images[labels == 1]\n",
    "qcd_images = images[labels == 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d87ee3-8c91-4344-9ec6-58e348a66dd9",
   "metadata": {},
   "source": [
    "And now we can plot a few examples to see how the images look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "003b6560-9c01-4758-a7d6-a3f5eb443992",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAEnCAYAAAC+IdmrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcU0lEQVR4nO3de5RlZXnn8e+vL3TLTRpBbAGDd+NKTOtqifGSGK/IJEucqFGTLEycITezcDSOqIlCJiaaeElmXGMGR4RExbtLYzSxVYxjLkCjzV0FEQLY0KIgDUjbl2f+2LvwdHNO16k6l6pd9f2sdVadevc+53l2dZ2nn3r3fs9JVSFJktRlKxY6AUmSpFHZ0EiSpM6zoZEkSZ1nQyNJkjrPhkaSJHWeDY0kSeo8GxrtJUklediYn/OlSb4yzuechiRPTXLDQuchaeFMoiZqMmxoFqm2Cbg0yV1Jbkryv5Pcd599HpHkI0luSfKDJJckeWWSlUmOa1+Id7S3m5N8OskzF+qYJC1tPfXmjiR7kvyw5/tfW+j8tLTZ0CxCSV4FvAV4NXBf4AnAccDnkqxu93kocD5wPfDTVXVf4AXARuCQnqc7rKoOBn4G2AR8IslLp3Mk45WGv7PSIlVVB8/cgP8Afrln7P0Lnd9cJFm10DlobvzPYZFJcihwBvAHVfWPVbWzqq4FXgg8BHhJu+sZwL9W1SuraitAVX2jql5SVbft+7xVdVNV/TVwOvCWWRqDE5Nc0878/OXMvkkemuSLSb7Xbnt/ksN6cj82yceTfLfd550DjvEvk3wlyX3b2aS3tc/37SQvb2eWVrX7finJm5L8C3AX8JAkT0xyYTsrdWGSJ/Y897VJntHz/elJ3tfen5m1OjnJf7QxX9+z732SnJ3k1iRXAI/fz89I0pCSrEnyV0m+097+KsmadttTk9yQ5HXta/La/c3mJHlwkn9Osj3JpiTv7HmN3+s0cW9NSHJ8kn9LcluSre1jD+jZt5L8fpKrgKv6xH5ykuuTPHUsPxiNlQ3N4vNEYC3w8d7BqroD+AzwrHboGcBH5/H8HwfuDzxyP/s8j2am53HAc4HfascD/DnwQOAngWNpGiSSrAQ+DVxHM5t0NPDB3idNsiLJu4HHAM+qqh8A/xV4DrChjXdSn3x+AziFZuZpO/APwP8E7ge8HfiHJPcb5uBbT6Y5/qcDb0jyk+34G4GHtrdnAyfP4TklDfZ6mpnmDTSzxccDf9Sz/QHAETR142TgzCSDatQHgIva/f8Hc3ud7gb+W/vYn6OpAb+3zz4nAT8LPLp3MMkJwLnAr1TVl+YQU1NiQ7P4HAHcUlW7+mzbChzZ3r9f+/1cfaf9evh+9nlLVX2/qv4D+CvgxQBVdXVVbaqqHVX1XZpm4hfaxxxP0+i8uqrurKq7q6r3QuDVNMXgcJpp6Lva8RcCf11VN1TVrcCb++RzdlVd3v5MngVcVVV/V1W7qupc4OvAL8/hZ3BGVf2wqi4GLqYpsDO5vKk99utpmiZJo/s14E+qaltbO86g+UOl1x+3teWfaf5oeeG+T5LkQTQzpzP7fhn4+2GTqKqLqurf29pxLfB/+HENm/HnbQ34Yc/YC9p9n1NVFwwbT9PlOcLF5xbgiCSr+jQ169vtAN9rv5+ro9uv39/PPtf33L+OplEhyVHAXwNPoZktWQHc2u53LHDdgEYM4GG0f5lV1Y96xh+4T7zrubfesQe2OfW6jh8f1zBu6rl/F3DwgFz2jSNpfvZ93d5TV1q3VtWd+9ne+zz99j12mCSSPILmD7GNwIE0/wdetM9u/WrQK4C/rarLhomjheEMzeLzb8AO4D/3DiY5mObUzJfaoc8DvzKP538esA34xn726S0OD+LHszp/BhTNRciHAr9OcxoKmiLwoP1cSHcl8JvAZ/eZSt4KHDMg9ozej4T/DvAT+2x/EHBje/9OmkI14wED8ulnK/c+dkmj2/d121tXANYlOWg/22dsHbDvjL1e/+2p8CN7tr+LZkb34W0Nex0/rmEzint7AXBSklP7bNMiYUOzyLTXlZwB/K8kJyRZneQ44MM0szMzKwXeCDyxvcD2AQBJHpbkfb0X6s5IclSSl7ePe21V7dlPGq9Osi7JscCpwIfa8UOAO4AfJDmaZhXWjAtois2bkxyUZG2SJ+1zbOfSFJDPp1mlRXtcpyY5us37NbP8iD4DPCLJS5KsSvKrNOe6P91u3wK8qP25bQSeP8vz9fow8Nr22I8B/mAOj5U02LnAHyU5MskRwBuA9+2zzxlJDkjyFOCXgI/s+yRVdR2wuWffJ7P36eZvAmuT/Kc0K0L/CFjTs/0Q4HbgjiSPAn53yPy/Q3O9zalJhn2MpsyGZhGqqr+g+Y//rTQXwX6b5q+OZ8xMtVbVt2guajsOuDzJD4CP0bzYt/c83W1J7gQuBU4EXlBVZ82SwidppmG30JzLfk87fgbNhbs/aMfvuXC5qnbTFJaH0SzXvAH41T7Hdg7wJ8AX20bt3cDngEuAr9E0LLtoLt7r97P5Hk2xexXNabf/DvxSVc2civtjmot6b23z/cAsx9rrDJrp62+3Of3dHB4rabA/palNl9DUoq+2YzNuonnNfofmj7bfqaqvD3iul9BctPt9mj/Q/nZmQ/sH4e8B/5dm1vZOmlo04w/bx2+nqT0fYkjtNYVPB05L8l+GfZymJ1X9Zte0mCT5TZom4Enti2rJSvIc4G+qat/TSpKWoHYJ9Puq6phZdh30+NOBh1XVr48xLXWQFwV3QFW9N8kumiXdS6qhSXIf4BdpZkSOovmL6xMLmpQkqXNsaDqiqpbq6Y/QnOr5EPBDmlNZb1jQjCRJneMpJ0mS1HleFCxJkjpvpIamXVb8jSRXJzltXElJUi9rjaTZzPuUU/uGRd8EnkmzLO5C4MVVdcWgxxyQNbWWgwZtltRh27n1lqo6cvY958Zao/nK6tV9x2vnzilnonG5mzv5Ue3Y980QgdEuCj4euLqqrgFI8kGaDzIcWGTWchA/m6ePEFLSYvX5+uikPirCWqN5WXVUv09PgF039nsTYnXB+fWFgdtGOeV0NHt/5sUNzO3zdCRpGNYaSbOa+LLtJKcApwCs3esjdiRpfKw10vI2ygzNjez9QX7H8OMPCLxHVZ1ZVRurauPqvT5SQ5KGYq2RNKtRGpoLgYcneXCSA4AXAZ8aT1qSdA9rjaRZzfuUU1Xtaj+9+Z+AlcBZVXX52DKTJKw1y9Gq9Q/oO17rDu07vvuKb/Yd9+Lf5WWka2iq6jM0n44sSRNjrZE0G98pWJIkdZ4NjSRJ6jwbGkmS1Hk2NJIkqfMm/sZ60nKTVf1fVrVr15Qzkbqp7t7Rd3z7I9f1Hb/Pug19x/MvW8aUkbrAGRpJktR5NjSSJKnzbGgkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z2Xb0pi5PFuLTdas6TteO/ovj15ou2+9te/4gZ84f07Ps+oBR/Ud33XTzXPOiRUr+4/v2T3359JEOEMjSZI6z4ZGkiR1ng2NJEnqPBsaSZLUeTY0kiSp81zlJElL3GJdzTRp81rNNIirmRY9Z2gkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6b6SLgpNcC2wHdgO7qmrjOJKSpF7WGkmzGccqp1+sqlvG8DyStD/WGvX1zbP697eP+K3NU85EC8lTTpIkqfNGbWgK+FySi5KcMo6EJKkPa42k/Rr1lNOTq+rGJPcHNiX5elV9uXeHtvicArCWA0cMJ2mZstZI2q+RZmiq6sb26zbgE8DxffY5s6o2VtXG1awZJZykZcpaI2k2825okhyU5JCZ+8CzgMvGlZgkgbVG0nBGOeV0FPCJJDPP84Gq+sexZLVcND+7e6uabh7S4mat0X65mkkwQkNTVdcAPzPGXCTpXqw1kobhsm1JktR5NjSSJKnzbGgkSVLn2dBIkqTOG8dnOWm+XM0kSdJYOEMjSZI6z4ZGkiR1ng2NJEnqPBsaSZLUeTY0kiSp82xoJElS59nQSJKkzrOhkSRJnWdDI0mSOs+GRpIkdZ4NjSRJ6jwbGkmS1Hk2NJIkqfNsaCRJUufZ0EiSpM6zoZEkSZ1nQyNJkjpv1oYmyVlJtiW5rGfs8CSbklzVfl032TQlLXXWGkmjGGaG5mzghH3GTgO+UFUPB77Qfi9Jozgba42keZq1oamqLwPf32f4ucA57f1zgJPGm5ak5cZaI2kU872G5qiq2trevwk4akz5SFIva42koYx8UXBVFVCDtic5JcnmJJt3smPUcJKWKWuNpP2Zb0Nzc5L1AO3XbYN2rKozq2pjVW1czZp5hpO0TFlrJA1lvg3Np4CT2/snA58cTzqStBdrjaShDLNs+1zg34BHJrkhycuANwPPTHIV8Iz2e0maN2uNpFGsmm2HqnrxgE1PH3MukpYxa42kUfhOwZIkqfNsaCRJUufZ0EiSpM6zoZEkSZ1nQyNJkjrPhkaSJHWeDY0kSeo8GxpJktR5NjSSJKnzbGgkSVLn2dBIkqTOm/WznCRJWu6yZk3f8dqxY8qZaBBnaCRJUufZ0EiSpM6zoZEkSZ1nQyNJkjrPhkaSJHWeq5yk2ST9x6umm4ekvqaxAsnVTIufMzSSJKnzbGgkSVLn2dBIkqTOs6GRJEmdN2tDk+SsJNuSXNYzdnqSG5NsaW8nTjZNSUudtUbSKIZZ5XQ28E7gb/cZf0dVvXXsGUmLjauZpuVsrDWaB1cgjW7F2rV9x/fcffeUM5m/WWdoqurLwPenkIukZcxaI2kUo1xD8/Ikl7TTxOvGlpEk7c1aI2lW821o3gU8FNgAbAXeNmjHJKck2Zxk806cFpQ0J9YaSUOZV0NTVTdX1e6q2gO8Gzh+P/ueWVUbq2rjavq/m6Mk9WOtkTSseTU0Sdb3fPs84LJB+0rSfFlrJA1r1lVOSc4FngockeQG4I3AU5NsAAq4FvjtyaWoScjqA/qO184fTTkTqWGt0bSsPPTQgdt23377FDNZPLq0mmmQWRuaqnpxn+H3TCAXScuYtUbSKHynYEmS1Hk2NJIkqfNsaCRJUufZ0EiSpM4b5rOctATNZzXT66/Z0nf8TQ/ZMFoykjRF01jJlFX9/3utXbvG+hj9mDM0kiSp82xoJElS59nQSJKkzrOhkSRJnWdDI0mSOs+GRpIkdZ7LtjU0l2dL0nDms9R6IZdnZ82avuO1Y8eUM5k/Z2gkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z4ZGkiR1nquctCStWLu27/ieu++eciaStPh1aTXTIM7QSJKkzrOhkSRJnWdDI0mSOs+GRpIkdd6sDU2SY5Ocl+SKJJcnObUdPzzJpiRXtV/XTT5dSUuVtUbSKIZZ5bQLeFVVfTXJIcBFSTYBLwW+UFVvTnIacBrwmsmlqq7Z85TH9h1f8f++NvnYrmbqImuNxmrlkUf2Hd/93e+OL8ahh/aPcfvtY4uh4cw6Q1NVW6vqq+397cCVwNHAc4Fz2t3OAU6aUI6SlgFrjaRRzOkamiTHAY8FzgeOqqqt7aabgKPGm5qk5cpaI2muhm5okhwMfAx4RVXtNZdWVQXUgMedkmRzks076f4b90iaLGuNpPkYqqFJspqmwLy/qj7eDt+cZH27fT2wrd9jq+rMqtpYVRtXs2YcOUtaoqw1kuZrmFVOAd4DXFlVb+/Z9Cng5Pb+ycAnx5+epOXCWiNpFMOscnoS8BvApUm2tGOvA94MfDjJy4DrgBdOJEN11qYPvbfv+I7a2Xf8pIf/wsDn2nPXXWPJSYuatUb7tfKw+/Yd333bD/qPj3E10yALuZppxUEH9R3fc+edU85kcZi1oamqrwAZsPnp401H0nJlrZE0Ct8pWJIkdZ4NjSRJ6jwbGkmS1Hk2NJIkqfOGWeUkAbBi7dq+44M+N+nZD9zQd/zuXz6+7/jauy6YV16SlodBq5lWHHhg3/FBqyPz+J8eGKMuvHTuiS2Q5bqaaRBnaCRJUufZ0EiSpM6zoZEkSZ1nQyNJkjrPhkaSJHWeDY0kSeo8l21raIOWZ8/V2r93ebY6JAM+XqpqunlMwMqHP2Tgtt1XXTPFTEYz1w+vXcil2auOPWbgtl3X3zDFTJYeZ2gkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z4ZGkiR1nqucNHVZ1f/XrnbtmnjsuX7AprQUVjMN0qWVTF3zoxMe33/DP1443USWEWdoJElS59nQSJKkzrOhkSRJnWdDI0mSOm/WhibJsUnOS3JFksuTnNqOn57kxiRb2tuJk09X0lJlrZE0imFWOe0CXlVVX01yCHBRkk3ttndU1Vsnl566YK6rlqaxmmkQVzMtatYajdXKdev6ju++9daJxz7A1UxTN2tDU1Vbga3t/e1JrgSOnnRikpYXa42kUczpGpokxwGPBc5vh16e5JIkZyXp3wpL0hxZayTN1dANTZKDgY8Br6iq24F3AQ8FNtD8VfW2AY87JcnmJJt3smP0jCUtadYaSfMxVEOTZDVNgXl/VX0coKpurqrdVbUHeDdwfL/HVtWZVbWxqjauZs248pa0BFlrJM3XMKucArwHuLKq3t4zvr5nt+cBl40/PUnLhbVG0iiGWeX0JOA3gEuTbGnHXge8OMkGoIBrgd+eQH7qgIVctaQlxVqjsZrGaiZNzopDDrnXWO4YPA8zzCqnrwDps+kzc0lMkvbHWiNpFL5TsCRJ6jwbGkmS1Hk2NJIkqfNsaCRJUucNs8pJkiQtQysOPLDv+J677pp47D3bt99rrHk7qv6coZEkSZ1nQyNJkjrPhkaSJHWeDY0kSeo8GxpJktR5NjSSJKnzXLYtSZL6GrQ8e+UR9+s7vvuW700ynf1yhkaSJHWeDY0kSeo8GxpJktR5NjSSJKnzbGgkSVLnucpJkqRZZM2avuO1Y8eUM1kcFnI10yDO0EiSpM6zoZEkSZ1nQyNJkjrPhkaSJHXerA1NkrVJLkhycZLLk5zRjj84yflJrk7yoSQHTD5dSUuVtUbSKIZZ5bQDeFpV3ZFkNfCVJJ8FXgm8o6o+mORvgJcB75pgrpKWNmvNlP3Td7YM3PbsB26YWh5dsPLII/qO77rhxonHXrF2bd/xPXff3f8BT3hM//F/v2RMGcHKdev6ju++9daxxZirWWdoqnFH++3q9lbA04CPtuPnACdNIkFJy4O1RtIohrqGJsnKJFuAbcAm4FvAbVW1q93lBuDoiWQoadmw1kiar6EamqraXVUbgGOA44FHDRsgySlJNifZvJPl+QZEkoZjrZE0X3Na5VRVtwHnAT8HHJZk5hqcY4C+JxKr6syq2lhVG1fT/50WJamXtUbSXA2zyunIJIe19+8DPBO4kqbYPL/d7WTgkxPKUdIyYK2RNIphVjmtB85JspKmAfpwVX06yRXAB5P8KfA14D0TzFPS0metmZBtL39i3/FnP3DKiXTYNFYzDTJwNdMgY1zNNMhCrmYaZNaGpqouAR7bZ/wamnPckjQya42kUfhOwZIkqfNsaCRJUufZ0EiSpM6zoZEkSZ03zConSVKH3f+d/zrxGFnd/zNDa+ePJh5bAmdoJEnSEmBDI0mSOs+GRpIkdZ4NjSRJ6jwbGkmS1Hk2NJIkqfNcti1JGpnLs0e34zmP7zu+5rMXTjmTxWHFIYfcayx3DJ6HcYZGkiR1ng2NJEnqPBsaSZLUeTY0kiSp82xoJElS57nKSZKkRWC5rmYaZM/27fcaq9ozcH9naCRJUufZ0EiSpM6zoZEkSZ1nQyNJkjpv1oYmydokFyS5OMnlSc5ox89O8u0kW9rbholnK2nJstZIGsUwq5x2AE+rqjuSrAa+kuSz7bZXV9VHJ5eepGXEWiNp3mZtaKqqgDvab1e3t5pkUpKWH2uNpFEMdQ1NkpVJtgDbgE1VdX676U1JLknyjiRrJpWkpOXBWiNpvoZqaKpqd1VtAI4Bjk/yU8BrgUcBjwcOB17T77FJTkmyOcnmnewYT9aSliRrjaT5mtMqp6q6DTgPOKGqtlZjB/Be4PgBjzmzqjZW1cbV+IeVpNlZayTN1TCrnI5Mclh7/z7AM4GvJ1nfjgU4CbhscmlKWuqsNZJGMcwqp/XAOUlW0jRAH66qTyf5YpIjgQBbgN+ZXJqSlgFrjaR7rDjwwHuN5YeD52GGWeV0CfDYPuNPm2NukjSQtUbSKHynYEmS1Hk2NJIkqfNsaCRJUufZ0EiSpM4bZpWTJEnSVO256657jVXtGbi/MzSSJKnzbGgkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z4ZGkiR1nsu2JUnSWPT7QEnovwR77LEnHkGSJGnCbGgkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z1VOkiRpLKaxmmkQZ2gkSVLn2dBIkqTOs6GRJEmdZ0MjSZI6z4ZGkiR1XqpqesGS7wLXtd8eAdwyteB7M7axjT1+P1FVR04x3kDWGmMbe8nGHlhnptrQ7BU42VxVG41tbGMvvdiLyXL9NzC2sZdD7F6ecpIkSZ1nQyNJkjpvIRuaM41tbGMv2diLyXL9NzC2sZdD7Hss2DU0kiRJ4+IpJ0mS1HkL0tAkOSHJN5JcneS0Kce+NsmlSbYk2TzhWGcl2Zbksp6xw5NsSnJV+3XdFGOfnuTG9ti3JDlxAnGPTXJekiuSXJ7k1HZ84se9n9jTOO61SS5IcnEb+4x2/MFJzm9/1z+U5IApxj47ybd7jnvDuGMvZtaZpVtn2jjWGmvN3qpqqjdgJfAt4CHAAcDFwKOnGP9a4Igpxfp54HHAZT1jfwGc1t4/DXjLFGOfDvzhhI95PfC49v4hwDeBR0/juPcTexrHHeDg9v5q4HzgCcCHgRe1438D/O4UY58NPH+Sx71Yb9aZpV1n2jjWGmvNXreFmKE5Hri6qq6pqh8BHwSeuwB5TFxVfRn4/j7DzwXOae+fA5w0xdgTV1Vbq+qr7f3twJXA0UzhuPcTe+KqcUf77er2VsDTgI+245M67kGxlzPrzBKuM21sa421Zi8L0dAcDVzf8/0NTOkXoVXA55JclOSUKcadcVRVbW3v3wQcNeX4L09ySTtVPJFp6BlJjgMeS9PFT/W494kNUzjuJCuTbAG2AZtoZghuq6pd7S4T+13fN3ZVzRz3m9rjfkeSNZOIvUhZZ5ZJnQFrDdYaYHleFPzkqnoc8Bzg95P8/EIlUs283TS723cBDwU2AFuBt00qUJKDgY8Br6iq23u3Tfq4+8SeynFX1e6q2gAcQzND8KhJxBkmdpKfAl7b5vB44HDgNdPKR9YZplBnwFqDteYeC9HQ3Agc2/P9Me3YVFTVje3XbcAnaH4ZpunmJOsB2q/bphW4qm5ufxn3AO9mQseeZDXNi/z9VfXxdngqx90v9rSOe0ZV3QacB/wccFiSVe2mif+u98Q+oZ0Wr6raAbyX6f+uLyTrzBKvM2CtsdbsbSEamguBh7dXZB8AvAj41DQCJzkoySEz94FnAZft/1Fj9yng5Pb+ycAnpxV45kXeeh4TOPYkAd4DXFlVb+/ZNPHjHhR7Ssd9ZJLD2vv3AZ5Jc179POD57W6TOu5+sb/eU9RDcz592r/rC8k6s4TrTBvHWmOt2du4ri6eyw04keaq8G8Br59i3IfQrHa4GLh80rGBc2mmHXfSnNN8GXA/4AvAVcDngcOnGPvvgEuBS2he9OsnEPfJNFO8lwBb2tuJ0zju/cSexnE/BvhaG+My4A09v3MXAFcDHwHWTDH2F9vjvgx4H+3qhOVys84s3TrTxrbWWGv2uvlOwZIkqfOW40XBkiRpibGhkSRJnWdDI0mSOs+GRpIkdZ4NjSRJ6jwbGkmS1Hk2NJIkqfNsaCRJUuf9f08PlPCnoCOMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAEnCAYAAAC+IdmrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAblUlEQVR4nO3de7SkVX3m8e/TTdOtQuQqQcCgqGMcJ7ZOh8RLMsYrMskSJ2rUxIWOM+RmFk6MI5pEIRMnmnhJZlxjBgcFo0HxtnSMSUTFuJxx0MY0yEUFFQRsaLkKIi00v/mj3oPVbVWfOqcu5+xT389atU6d/b5Ve7+nT/36Oft9d1WqCkmSpJatW+kBSJIkjctAI0mSmmegkSRJzTPQSJKk5hloJElS8ww0kiSpeQYa7SZJJXnohJ/zxUk+P8nnnIUkT0pyzUqPQ9LKmUZN1HQYaFapLgR8JckdSa5L8j+S3H+PfR6e5ANJbkhya5KLkvx+kvVJju5eiLd3t+uTfDzJ01bqmCStbX315vYk9yT5Qd/3v77S49PaZqBZhZK8Angj8Erg/sDPA0cDn0yyodvnGOB84GrgX1XV/YHnAluA/fue7oCq2g94NHAu8JEkL57NkUxWevydlVapqtpv4QZ8G/iVvrb3rvT4liLJPis9Bi2N/zmsMkl+AjgN+L2q+oeququqrgSeBzwEeGG362nA/62q36+q7QBV9bWqemFV3bLn81bVdVX1V8CpwBsXCQbHJ/lmN/PzFwv7JjkmyWeS3Nhte2+SA/rGflSSDyf5brfP24Yc418k+XyS+3ezSW/unu9bSV7WzSzt0+372SSvT/J/gDuAhyR5fJIvdbNSX0ry+L7nvjLJU/u+PzXJe7r7C7NWJyb5dtfnH/bte58kZya5OcmlwM/u5WckaURJNib5yyTf6W5/mWRjt+1JSa5J8pruNXnl3mZzkjw4yT8luS3JuUne1vca/7HTxP01IcmxSb6Q5JYk27vH7tu3byX53SSXA5cP6PuJSa5O8qSJ/GA0UQaa1efxwCbgw/2NVXU78Ang6V3TU4EPLuP5Pww8APgXe9nn2fRmeh4LPAv49117gD8DHgj8NHAUvYBEkvXAx4Gr6M0mHQG8r/9Jk6xL8g7gZ4CnV9WtwH8Engls7vo7YcB4XgScRG/m6Tbg74D/BhwMvAX4uyQHj3LwnSfSO/6nAK9N8tNd++uAY7rbM4ATl/Cckob7Q3ozzZvpzRYfC/xR3/afBA6hVzdOBE5PMqxG/S1wQbf/f2Fpr9NdwH/qHvs4ejXgd/bY5wTg54BH9jcmOQ44G/jVqvrsEvrUjBhoVp9DgBuq6u4B27YDh3b3D+6+X6rvdF8P2ss+b6yqm6rq28BfAi8AqKorqurcqtpZVd+lFyb+TfeYY+kFnVdW1fer6s6q6r8QeAO9YnAQvWnoO7r25wF/VVXXVNXNwBsGjOfMqrqk+5k8Hbi8qv6mqu6uqrOBrwK/soSfwWlV9YOquhC4kF6BXRjL67tjv5peaJI0vl8H/qSqdnS14zR6f6j0++OutvwTvT9anrfnkyR5EL2Z04V9Pwf871EHUVUXVNX/62rHlcD/5Ec1bMGfdTXgB31tz+32fWZVfXHU/jRbniNcfW4ADkmyz4BQc3i3HeDG7vulOqL7etNe9rm67/5V9IIKSQ4D/gr4BXqzJeuAm7v9jgKuGhLEAB5K95dZVf2wr/2Be/R3NT+uv+2B3Zj6XcWPjmsU1/XdvwPYb8hY9uxH0vLs+bq9t650bq6q7+9le//zDNr3qFEGkeTh9P4Q2wLcl97/gRfssdugGvRy4N1VdfEo/WhlOEOz+nwB2An8u/7GJPvROzXz2a7pU8CvLuP5nw3sAL62l336i8OD+NGszn8Fit5FyD8B/Aa901DQKwIP2suFdJcBLwH+fo+p5O3AkUP6XtD/kfDfAX5qj+0PAq7t7n+fXqFa8JNDxjPIdn782CWNb8/XbX9dATgwyf32sn3B9iH7Ltjt9d+dCj+0b/vb6c3oPqyrYa/hRzVsQfHjnguckOTkAdu0ShhoVpnuupLTgP+e5LgkG5IcDZxDb3ZmYaXA64DHdxfY/iRAkocmeU//hboLkhyW5GXd415dVffsZRivTHJgkqOAk4H3d+37A7cDtyY5gt4qrAVfpFds3pDkfkk2JXnCHsd2Nr0C8qn0VmnRHdfJSY7oxv2qRX5EnwAenuSFSfZJ8mv0znV/vNu+DXh+93PbAjxnkefrdw7w6u7YjwR+bwmPlTTc2cAfJTk0ySHAa4H37LHPaUn2TfILwC8DH9jzSarqKmBr375PZPfTzV8HNiX5t+mtCP0jYGPf9v2B7wG3J3kE8Nsjjv879K63OTnJqI/RjBloVqGq+nN6//G/id5FsN+i91fHUxemWqvqG/QuajsauCTJrcCH6L3Yb+t7uluSfB/4CnA88NyqeuciQ/govWnYbfTOZZ/RtZ9G78LdW7v2ey9crqpd9ArLQ+kt17wG+LUBx3YW8CfAZ7qg9g7gk8BFwD/TCyx307t4b9DP5kZ6xe4V9E67/Wfgl6tq4VTcH9O7qPfmbrx/u8ix9juN3vT1t7ox/c0SHitpuD+lV5suoleLvty1LbiO3mv2O/T+aPutqvrqkOd6Ib2Ldm+i9wfauxc2dH8Q/g7wv+jN2n6fXi1a8Afd42+jV3vez4i6awqfApyS5D+M+jjNTqoGza5pNUnyEnoh4Andi2rNSvJM4K+ras/TSpLWoG4J9Huq6shFdh32+FOBh1bVb0xwWGqQFwU3oKreleRueku611SgSXIf4JfozYgcRu8vro+s6KAkSc0x0DSiqtbq6Y/QO9XzfuAH9E5lvXZFRyRJao6nnCRJUvO8KFiSJDVvrEDTLSv+WpIrkpwyqUFJUj9rjaTFLPuUU/eGRV8HnkZvWdyXgBdU1aXDHrNvNtYm7jdss6SG3cbNN1TVoYvvuTTWGkkL7uT7/LB27vlmiMB4FwUfC1xRVd8ESPI+eh9kOLTIbOJ+/FyeMkaXklarT9UHp/VREdYaSQCcX58eum2cU05HsPtnXlzD0j5PR5JGYa2RtKipL9tOchJwEsCm3T5iR5Imx1ojzbdxZmiuZfcP8juSH31A4L2q6vSq2lJVWzbs9pEakjQSa42kRY0TaL4EPCzJg5PsCzwf+NhkhiVJ97LWSFrUsk85VdXd3ac3/yOwHnhnVV0ysZFJEtYaaUVl4IIiWIVvyjvWNTRV9Ql6n44sSVNjrZG0GN8pWJIkNc9AI0mSmmegkSRJzTPQSJKk5k39jfU0G+sPuP/A9l233DrjkUiS1oxVuJppGGdoJElS8ww0kiSpeQYaSZLUPAONJElqnoFGkiQ1z0AjSZKa57LtNWJiy7PXrR++7Z5dk+lDkqQJc4ZGkiQ1z0AjSZKaZ6CRJEnNM9BIkqTmGWgkSVLzXOWk3bmSSZLUIGdoJElS8ww0kiSpeQYaSZLUPAONJElq3lgXBSe5ErgN2AXcXVVbJjEoSepnrZG0mEmscvqlqrphAs8jDZcMbq+a7Ti0kqw1Gij/+l8ObK8LLpnxSLSSPOUkSZKaN26gKeCTSS5IctIkBiRJA1hrJO3VuKecnlhV1yZ5AHBukq9W1ef6d+iKz0kAm7jvmN1JmlPWGkl7NdYMTVVd233dAXwEOHbAPqdX1Zaq2rKBjeN0J2lOWWskLWbZgSbJ/ZLsv3AfeDpw8aQGJklgrZE0mnFOOR0GfCS91Sf7AH9bVf8wkVFJe3I10zyz1miv1sxqpqWu5nT1526WHWiq6pvAoyc4Fkn6MdYaSaNw2bYkSWqegUaSJDXPQCNJkppnoJEkSc2bxGc5SYN5Bb6kNe6mlzxuYPtB7/rC0p9sqbXRWrobZ2gkSVLzDDSSJKl5BhpJktQ8A40kSWqegUaSJDXPQCNJkprnsm1Nj0sKJa1xy1qeralwhkaSJDXPQCNJkppnoJEkSc0z0EiSpOYZaCRJUvNc5SRJGls2bhzYXjt3zngkmlfO0EiSpOYZaCRJUvMMNJIkqXkGGkmS1LxFA02SdybZkeTivraDkpyb5PLu64HTHaaktc5aI2kco8zQnAkct0fbKcCnq+phwKe771e3desH3yStFmeyFmrNnKqdOwfe1or1Bx448KbVY9FAU1WfA27ao/lZwFnd/bOAEyY7LEnzxlojaRzLvYbmsKra3t2/DjhsQuORpH7WGkkjGfui4KoqoIZtT3JSkq1Jtt7F2pl+lDRb1hpJe7PcQHN9ksMBuq87hu1YVadX1Zaq2rKBwe8kKUlDWGskjWS5geZjwInd/ROBj05mOJK0G2uNpJEs+llOSc4GngQckuQa4HXAG4BzkrwUuAp43jQHORH37FrpEUjaizVTa7QmfeKS8wa2P+OBm2c7EA21aKCpqhcM2fSUCY9F0hyz1kgah+8ULEmSmmegkSRJzTPQSJKk5hloJElS8xa9KFhtW3/A/Qe277rl1hmPRJLateTVTMM+K9AVt1PjDI0kSWqegUaSJDXPQCNJkppnoJEkSc0z0EiSpOYZaCRJUvNctr3GuTxbklaAy7NnzhkaSZLUPAONJElqnoFGkiQ1z0AjSZKaZ6CRJEnNc5XTvEoGt1fNdhySJE2AMzSSJKl5BhpJktQ8A40kSWqegUaSJDVv0UCT5J1JdiS5uK/t1CTXJtnW3Y6f7jAlrXXWGknjGGWV05nA24B379H+1qp608RHpNkYsppp/WEPGPqQXd+9cfCGGXxmyfpDDx3Yvuu7313aE7m6azU7E2uNpGVadIamqj4H3DSDsUiaY9YaSeMY5xqalyW5qJsmPnBiI5Kk3VlrJC1quYHm7cAxwGZgO/DmYTsmOSnJ1iRb72LnMruTNKesNZJGsqxAU1XXV9WuqroHeAdw7F72Pb2qtlTVlg1sXO44Jc0ha42kUS0r0CQ5vO/bZwMXD9tXkpbLWiNpVIuuckpyNvAk4JAk1wCvA56UZDNQwJXAb05viJqG9QcfNLB91/U7ZjyS0Sx1NdP6AwdfarHr5psnMRxNgbVG0jgWDTRV9YIBzWdMYSyS5pi1RtI4fKdgSZLUPAONJElqnoFGkiQ1z0AjSZKaN8pnOakFS/yMol03ru13mHc1kyTNF2doJElS8ww0kiSpeQYaSZLUPAONJElqnoFGkiQ1z0AjSZKa57LttWLI8uxZyMaNA9tr584Zj0SSNK+coZEkSc0z0EiSpOYZaCRJUvMMNJIkqXkGGkmS1DxXOWlsrmaS5oerGrVaOUMjSZKaZ6CRJEnNM9BIkqTmGWgkSVLzFg00SY5Kcl6SS5NckuTkrv2gJOcmubz7euD0hytprbLWSBrHKKuc7gZeUVVfTrI/cEGSc4EXA5+uqjckOQU4BXjV9Iaq1qx/+DED2+vb1w5sv3vLI4Y+17rPb5vEkLS6WWsa4GomrVaLztBU1faq+nJ3/zbgMuAI4FnAWd1uZwEnTGmMkuaAtUbSOJZ0DU2So4HHAOcDh1XV9m7TdcBhkx2apHllrZG0VCMHmiT7AR8CXl5V3+vfVlUF1JDHnZRka5Ktd+FUpaS9s9ZIWo6RAk2SDfQKzHur6sNd8/VJDu+2Hw7sGPTYqjq9qrZU1ZYNDH6HSUkCa42k5RtllVOAM4DLquotfZs+BpzY3T8R+OjkhydpXlhrJI1jlFVOTwBeBHwlybau7TXAG4BzkrwUuAp43lRGqFVv/cMeMrB919e/MbD9H7+zbWD7Mx545/A+Djl4cB833Lj3wakl1hpJy7ZooKmqzwMZsvkpkx2OpHllrZE0Dt8pWJIkNc9AI0mSmmegkSRJzTPQSJKk5o2yykkCYP0jHz6wvfZZWi5+xgM3L7nvia1mypBrTmvge7VJTcnGwe+/4+cvaR44QyNJkppnoJEkSc0z0EiSpOYZaCRJUvMMNJIkqXkGGkmS1DyXba9x6+5734Ht99xxx5Kfa9elXx93OCvP5dlaw1yerXnmDI0kSWqegUaSJDXPQCNJkppnoJEkSc0z0EiSpOa5ymmNW85qJkmSWuMMjSRJap6BRpIkNc9AI0mSmmegkSRJzVs00CQ5Ksl5SS5NckmSk7v2U5Ncm2Rbdzt++sOVtFZZaySNY5RVTncDr6iqLyfZH7ggybndtrdW1ZumNzxpNrJh36Hb6q4fznAkc81aM6YbX/q4ge0/cfVdA9s3fHLrNIcjzdSigaaqtgPbu/u3JbkMOGLaA5M0X6w1ksaxpGtokhwNPAY4v2t6WZKLkrwzyYGTHpyk+WStkbRUIweaJPsBHwJeXlXfA94OHANspvdX1ZuHPO6kJFuTbL0LP9pe0t5ZayQtx0iBJskGegXmvVX1YYCqur6qdlXVPcA7gGMHPbaqTq+qLVW1ZQMbJzVuSWuQtUbSco2yyinAGcBlVfWWvvbD+3Z7NnDx5IcnaV5YaySNY5RVTk8AXgR8Jcm2ru01wAuSbAYKuBL4zSmMT5oJVzKtCtaaPVz1J4NXLf3Ua78wsP3gMwa3D7Nu06ah2+65884lPdeakQxur5rtOLRko6xy+jww6F/4E5MfjqR5Za2RNA7fKViSJDXPQCNJkppnoJEkSc0z0EiSpOaNsspJcyRbHjV0W211taw0S8NWM03KclYyrdt//8HPddtt4w5ndXA1U7OcoZEkSc0z0EiSpOYZaCRJUvMMNJIkqXkGGkmS1DwDjSRJap7LtrUbl2ZL82Pdfe87dNs9d9wxuH2tLM/WmuMMjSRJap6BRpIkNc9AI0mSmmegkSRJzTPQSJKk5rnKSZLWuHWbNg1sH7aSSWqRMzSSJKl5BhpJktQ8A40kSWqegUaSJDVv0UCTZFOSLya5MMklSU7r2h+c5PwkVyR5f5J9pz9cSWuVtUbSOEaZodkJPLmqHg1sBo5L8vPAG4G3VtVDgZuBl05tlGpTMvgmDWatmZJ77rxz4G2S1m3aNPAmzcqigaZ6bu++3dDdCngy8MGu/SzghGkMUNJ8sNZIGsdI19AkWZ9kG7ADOBf4BnBLVd3d7XINcMRURihpblhrJC3XSIGmqnZV1WbgSOBY4BGjdpDkpCRbk2y9i53LG6WkuWCtkbRcS1rlVFW3AOcBjwMOSLLwTsNHAtcOeczpVbWlqrZsYOM4Y5U0J6w1kpZqlFVOhyY5oLt/H+BpwGX0is1zut1OBD46pTFKmgPWGknjGOWznA4Hzkqynl4AOqeqPp7kUuB9Sf4U+GfgjCmOUy2qWukRqC3WmoZNetWUtFSLBpqqugh4zID2b9I7xy1JY7PWSBqH7xQsSZKaZ6CRJEnNM9BIkqTmGWgkSVLzDDSSJKl5BhpJktQ8A40kSWqegUaSJDXPQCNJkppnoJEkSc0z0EiSpOaN8uGUkqRRrVs/uP2eXbMdxxzJhn0HttddP5zxSLSSnKGRJEnNM9BIkqTmGWgkSVLzDDSSJKl5BhpJktQ8VzlJ0iS5mmnmVuVqJle7zZwzNJIkqXkGGkmS1DwDjSRJap6BRpIkNW/RQJNkU5IvJrkwySVJTuvaz0zyrSTbutvmqY9W0pplrZE0jlFWOe0EnlxVtyfZAHw+yd93215ZVR+c3vAkzRFrjdaOtb6aaRWu4lo00FRVAbd3327objXNQUmaP9YaSeMY6RqaJOuTbAN2AOdW1fndptcnuSjJW5NsnNYgJc0Ha42k5Rop0FTVrqraDBwJHJvkUcCrgUcAPwscBLxq0GOTnJRka5Ktd7FzMqOWtCZZayQt15JWOVXVLcB5wHFVtb16dgLvAo4d8pjTq2pLVW3ZgH9YSVqctUbSUo2yyunQJAd09+8DPA34apLDu7YAJwAXT2+YktY6a42kcYyyyulw4Kwk6+kFoHOq6uNJPpPkUCDANuC3pjdMSXPAWiM1IhsGx4faubpXOV0EPGZA+5OnMiJJc8laI2kcvlOwJElqnoFGkiQ1z0AjSZKaZ6CRJEnNG2WVkyRJ0r1q5+p780pnaCRJUvMMNJIkqXkGGkmS1DwDjSRJap6BRpIkNc9AI0mSmmegkSRJzTPQSJKk5hloJElS8ww0kiSpeQYaSZLUPAONJElqnh9OKUmSpiob9h26re764UT6cIZGkiQ1z0AjSZKaZ6CRJEnNM9BIkqTmGWgkSVLzUlWz6yz5LnBV9+0hwA0z63x39m3f9j15P1VVh86wv6GsNfZt32u276F1ZqaBZreOk61VtcW+7du+117fq8m8/hvYt33PQ9/9POUkSZKaZ6CRJEnNW8lAc7p927d9r9m+V5N5/Tewb/ueh77vtWLX0EiSJE2Kp5wkSVLzViTQJDkuydeSXJHklBn3fWWSryTZlmTrlPt6Z5IdSS7uazsoyblJLu++HjjDvk9Ncm137NuSHD+Ffo9Kcl6SS5NckuTkrn3qx72Xvmdx3JuSfDHJhV3fp3XtD05yfve7/v4kwz+hbfJ9n5nkW33HvXnSfa9m1pm1W2e6fqw11prdVdVMb8B64BvAQ4B9gQuBR86w/yuBQ2bU1y8CjwUu7mv7c+CU7v4pwBtn2PepwB9M+ZgPBx7b3d8f+DrwyFkc9176nsVxB9ivu78BOB/4eeAc4Pld+18Dvz3Dvs8EnjPN416tN+vM2q4zXT/WGmvNbreVmKE5Friiqr5ZVT8E3gc8awXGMXVV9Tngpj2anwWc1d0/Czhhhn1PXVVtr6ovd/dvAy4DjmAGx72Xvqeuem7vvt3Q3Qp4MvDBrn1axz2s73lmnVnDdabr21pjrdnNSgSaI4Cr+76/hhn9InQK+GSSC5KcNMN+FxxWVdu7+9cBh824/5cluaibKp7KNPSCJEcDj6GX4md63Hv0DTM47iTrk2wDdgDn0pshuKWq7u52mdrv+p59V9XCcb++O+63Jtk4jb5XKevMnNQZsNZgrQHm86LgJ1bVY4FnAr+b5BdXaiDVm7ebZbp9O3AMsBnYDrx5Wh0l2Q/4EPDyqvpe/7ZpH/eAvmdy3FW1q6o2A0fSmyF4xDT6GaXvJI8CXt2N4WeBg4BXzWo8ss4wgzoD1hqsNfdaiUBzLXBU3/dHdm0zUVXXdl93AB+h98swS9cnORyg+7pjVh1X1fXdL+M9wDuY0rEn2UDvRf7eqvpw1zyT4x7U96yOe0FV3QKcBzwOOCDJPt2mqf+u9/V9XDctXlW1E3gXs/9dX0nWmTVeZ8BaY63Z3UoEmi8BD+uuyN4XeD7wsVl0nOR+SfZfuA88Hbh474+auI8BJ3b3TwQ+OquOF17knWczhWNPEuAM4LKqekvfpqkf97C+Z3TchyY5oLt/H+Bp9M6rnwc8p9ttWsc9qO+v9hX10DufPuvf9ZVknVnDdabrx1pjrdndpK4uXsoNOJ7eVeHfAP5whv0+hN5qhwuBS6bdN3A2vWnHu+id03wpcDDwaeBy4FPAQTPs+2+ArwAX0XvRHz6Ffp9Ib4r3ImBbdzt+Fse9l75ncdw/A/xz18fFwGv7fue+CFwBfADYOMO+P9Md98XAe+hWJ8zLzTqzdutM17e1xlqz2813CpYkSc2bx4uCJUnSGmOgkSRJzTPQSJKk5hloJElS8ww0kiSpeQYaSZLUPAONJElqnoFGkiQ17/8DKt1nokNWItQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    plt.figure(figsize=(8, 4))\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.title(\"QCD background\")\n",
    "    plt.imshow(qcd_images[i])\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.title(\"Top quark\")\n",
    "    plt.imshow(top_images[i])\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69ad82ef-8c05-4aa6-9335-316f2b9c21f5",
   "metadata": {},
   "source": [
    "Notice how the background events are much more clustered around the centre, while the top-quark jets have multiple \"prongs\". We can see this more clearly by averaging the pixel values across all events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e044c284-37ee-4934-83e9-408d1b6fe4f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAEnCAYAAAC61AmWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs5UlEQVR4nO3de5BtZ1nn8d+v+/TtXJKTkJAJSeQmpaJAwGO8ocNwcQCtIVgMEi9EBw1UyRTMMDNkmCqBURysElALSwlDSFQQYrilGLwEjIOIBE5iCCHRSbjEJJzkEM790n26ez/zx95HOmc/z0mvnN17d6/9/VSdOt1vr37Xu9Ze++m3V7/PehwRAgAAaIOJUQ8AAABgUJjYAACA1mBiAwAAWoOJDQAAaA0mNgAAoDWY2AAAgNZgYgMAAFqDic2Ysf03tvfanhn1WAbB9pNtX2d7v+2Dtv/a9g+dsM207TfZvtP2Ydtft32l7cf1vv43tud733/A9k22L2/LOQJGrU1xpxdL/mTU40CNic0Y6f0g/zFJIenfrUH/mwbd58Ps74mS/k7SlyQ9XtJjJH1U0vW2L1qx6bXqHu/PSjpd0tMk3STpOSu2eXVEbJN0rqTXSXqZpE/Y9hofBtBqbYs7WP+Y2IyXl0v6nKSrJF0qSbZnbO+z/X3HN7J9tu2jth/d+/ynbN/S2+6ztp+6Ytuv23697VslHba9qXe34yu9OyC3237xiu0nbb/N9oO2v2b71bbjeHCyfbrt99jeZfs+279he7I4njdJ+vuI+B8RsSciDkbE70n6E0m/1evvuZKeJ+lFEfGFiFiKiP0R8fsR8Z4TO4yIwxHxN+oG4B+W9JOP7FQD6GlN3LH9fElvkPQztg/Z/mKv/TG9O8d7bN9l+1dWfM+bbF9r+4O9sd1s+2kDPcN4CCY24+Xlkt7X+/dvbZ8TEQuSPizpkhXbvVTS/42I3bafLulKSa+U9ChJ75J03Qm3lC9RdwKwPSKWJH1F3d/QTpf0Zkl/Yvvc3ra/IukFki6U9AxJF58wxqskLUn6TklPl/QTkn65OJ7nSfqzpP0aST9me1bScyV9PiLuKfpIRcQ/S9rZOw4Aj1xr4k5E/IWk35T0wYjYGhHHJygfkHSvuneNXyLpN20/e8W3vkjdWHWmpPdL+qjtqfqU4VQwsRkTtp8p6bGSromIm9QNAj/b+/L71f3Ty3E/22uTpMskvSsiboyI5Yi4WtKCpJXrWH4vIu6JiKOSFBF/FhHfiIhORHxQ0p2Sjv9p6KWSfjci7o2IvZLeumKM50h6oaTX9u6c7Jb0jhPGttJZknYl7bskTaobRB5VbLMa3+j1AeARaGncOfEYL5D0o5JeHxHzEXGLpP+t7oTuuJsi4tqIWJT0dkmzJxwLBoiJzfi4VNJfRcSDvc/f32uTpBskbbb9g72/h18o6SO9rz1W0ut6t4P32d4n6QJ1fzM57iF3Q2y/fMUt5H2Svk/dSYh633dP8b2PlTQladeK732XpEcXx/SgumtiTnSuun/P/1bvX7bNapwnac8j/F4A7Yw7J3qMpD0RcXBF293qxo++/UVER9++u4M1wKKrMWB7Tt3fWCZt399rnpG03fbTIuKLtq9R99buA5I+vuJNeo+kt0TEW06yi38pEW/7sZLere7C3L+PiGXbt0g6vgh3l6TzV3zvBSs+vkfd38rO6t1afjiflPTvJb33hPaXSvpcRCzY/qSk19g+PyLuXUWfx4/jAknfr95aHQDNtDjuxAmff0PSmba3rRj/d0i6L9uf7YneWL6xin3hEeCOzXi4WNKypCer+1vRhZK+R9Lf6tu3S98v6Wck/Zy+fTtY6gaLV/V+q7LtLbZ/0va2Yl9b1H3jf1OSbP+Sur85HXeNuhON82xvl/T641+IiF2S/krS22yfZnvC9hNt/+tiX2+W9CO232L7TNvbbP9HSb8k6dd6fX5S0vWSPmL7+3uLDLfZfpXt/3Bih7Y39/b3MUmfl/SJYt8ATu5itTPuPCDpcb0Jinrr9z4r6X/ZnnV3kfMr1E1iOO77bf90b7Hya9WdSH2u6B+niInNeLhU0nsj4p8j4v7j/yS9U9LP2d4UETdKOqzu7dE/P/6NEbFT3YV375S0V9Jdkn6x2lFE3C7pbZL+Xt0A8BR1U7KPe7e6QeRWSf+g7sRhSd0AKHUD3rSk23v7u1bFn5Ii4k5Jz1Q3ffvrkvZJ+nVJL+5NaI57SW8/H5S0X9Jtknaoe8fnuHfaPtgb8+9I+pCk5/duGwNorpVxR99OWPiW7Zt7H18i6XHq3oX5iKQ3nhCDPqbuBG6vpF+Q9NO99TZYA4448a4aMDy2XyDpDyPisQPo63x1fwt6Y5bKDQDSYOPOKvb1JknfGRE/v9b7Qhd3bDBUtudsv7D3J6HzJL1R314weEp6a2heIOlc21sH0SeAjW8t4w7WHyY2GDaruzZmr7q3hO9Qbz3MIETElyLiNyLi0KD6BLDhrWncwfrCn6IAAEBrcMcGAAC0xrqY2Nh+vu1/crfGxuWjHg+AdiLWAO038j9FuVto7P+pW/fnXklfkHRJL30vNe2ZmNWWIY0QQBMHtffBiDh71OM4EbHmEWhQ3N4TDX9Pnij6doN+JvNto+q7aq9+DiYPe3C17XLDJ0MsV88CTMZY7LP++d3w5/oGXJEyr8M6FgvpC7oenjx8kaS7IuKrkmT7A+oWDCuDzay26Af9nCEND2uuQfAssVZs3fhkXHv3qMdQWF+xZhDXfdn3YG7Ge7KvwHVpYstcs75nZ/MvzBXtic5pm/P2uby+5PJc/iPPnTx+TB7pn3z4aP74mYnDR9P2clKyZ1++fTL5isV8EhTHikfhNHz8VmTHX/VRxdrqem6yfYM4fmN8qvzaevhT1Hl6aN2Oe/XQGhsAMAjEGmAMrIc7Nqti+zJ1K75qVvksHQBOFbEG2NjWwx2b+/TQgmTn66HFwyRJEXFFROyIiB1Tmhna4AC0BrEGGAPrYWLzBUlPsv1429OSXibpuhGPCUD7EGuAMTDyP0VFxJLtV0v6S0mTkq6MiC+PeFhYC2u5WBJ4GMQayUVWULp4VJKnih8RnSxdKP892ZuLRcVZH5I0v5A2L/+rR/X3vbicbCl1NuVjWZrNF0Mvz+bbT27pP/7Zb6abnmQhb3Fui/MSi0k/xeLhMstLxaLv5fx8NVpsPFH03bRe8BomfIx8YiNJEfEJdautAsCaIdYA7bce/hQFAAAwEExsAABAazCxAQAArcHEBgAAtMa6WDyMlhnEo7WBjajJY+LX8v1QZKhEkRTj6eli+/wbvKn/R0ccO5Z3XmTueCove6CqPdHZmo+7M51n7hw9K/+RN3U0P19Lc/2/+y+elu9zYrbo+4EDabuScyhJmp/va6pKW8RSVW8qV2W/Zcp9VplVZUcNfx4MAHdsAABAazCxAQAArcHEBgAAtAYTGwAA0BosHsbgjWKR8AgWqAF9suttHV2b1SLhUrXYdGL1vxOXi4QrVRmHpEzCxKFiwfLmfJ8zB/OFrwvb8oWyE0v9xz+xkPcxOV8s5J3Mz1VaOkGSZpLCq8UiYRcLtmO5YXmDtI+Gi4TXUazljg0AAGgNJjYAAKA1mNgAAIDWYGIDAABag4kNAABoDbKi0A7raEU+WqZJmYRT3fak48h/D3VSsqDJo/NPusvpBhlNnaKMQ5XRMzeX99Pgvbx8WpJBJEnFqT1yVp791Cl+Es7t6W9bOCvPLJv7RlF+YqHIftqcH398a29/Y1GWospOc4OsNakotVGU5dgIuGMDAABag4kNAABoDSY2AACgNZjYAACA1lgXi4dtf13SQUnLkpYiYsdoRwSgjYg1QPuti4lNz7+JiAdHPQgMwICyPUobeLU+1oVmsaZJxt0IsvOyDKgsU0qSYjHPUJoosp/KjKak5pSrGk/VOSm272yezbdvYGkuz35ansnPy9LmoqNk87lv5fFq/tF5ltOWfUfyro/M5/ucTMZe1G2qXs8qRpbZctn2VVxuGn9HUCuNP0UBAIDWWC8Tm5D0V7Zvsn3ZqAcDoLWINUDLrZc/RT0zIu6z/WhJ19v+x4j49MoNekHoMkmaVXXfEABOilgDtNy6uGMTEff1/t8t6SOSLkq2uSIidkTEjikVT5sEgJMg1gDtN/KJje0ttrcd/1jST0i6bbSjAtA2xBpgPKyHP0WdI+kj7q6c3iTp/RHxF6MdElataQbUWhnByntsOGsfa5rUiir7OPUsweisvq6UVGc/NR5LE0eOps3VHpdn+rOFlmfzuk1Vrajpg/lrcfTR+Te4099elX6aPJa3x1w+Ri8U37Ap+bFcZCJ5Mj9bEcXx5MlieXZVJ8/E0kTVSbH9CIx8YhMRX5X0tFGPA0C7EWuA8TDyP0UBAAAMChMbAADQGkxsAABAazCxAQAArTHyxcPY4LKMj0FlSjWpSUL2E9aDtXw/VJLMJWf1hk4iilpEnqqya5J9zuY1nmK5eB9vzmsr6ViedjRxtL99qji3kwv5j7ZDj8n3GZvy+DGzN+s733Z6X1FX68hC2q4qKypT1HiKtYx7ZaZpw1pRI4jN3LEBAACtwcQGAAC0BhMbAADQGkxsAABAa7B4GKszipIF1ePcs8VrlFTAejCAhcJl2YNiAamnkjBeLAZuus8m/cRiUWsge1y/pOqd6c35IuTsvTxxcL7YNi9eunl33t4pFlsvnNHfNvtgvstNh4vjnyjiWLV4eu/+vqaq5EW1SDyO5QuTq2so77xB/F1nuGMDAABag4kNAABoDSY2AACgNZjYAACA1mBiAwAAWoOsKJyatX5cfLrPZD6+AVbqY0xV2SWNty8ylJLMpTKDqiq1UGRFeXo63z7rv1NkP1UZPdNTed9FSQUnpRliNh9fVn5Bkibn8zgxtyc//k5yulwkFi3P5Oe2ejXjwMHiCw2yPtdSFVOrTNOmmanZ9gPKYuWODQAAaA0mNgAAoDWY2AAAgNZgYgMAAFpjaBMb21fa3m37thVtZ9q+3vadvf+TB1gDwOoRa4DxNsysqKskvVPSH61ou1zSpyLirbYv733++iGOCavVdCV8um0+j25aGyffmJpQ+BdXaT3Fmk6RzVS8d6JZmSdFkhU1MZfXIYqiblOa5XSyfRaZThnP5PWZGr9nk+PMMqUkSfN5raTZ3UfS9qlDeXZVZ7o/Zk0eyY99ssjEcpHlFVWdp6NF/ats2+J1aBQ7lcfg8jqcKDLrmmamrmHMHtodm4j4tKQ9JzS/SNLVvY+vlnTxsMYDoJ2INcB4G/Uam3MiYlfv4/slnTPKwQBoLWINMCZGPbH5FxERqivZy/Zltnfa3rmohSGODECbEGuAdhv1xOYB2+dKUu//3dWGEXFFROyIiB1TKv5mCwA5Yg0wJkY9sblO0qW9jy+V9LERjgVAexFrgDExtKwo238q6VmSzrJ9r6Q3SnqrpGtsv0LS3ZJeOqzxYP1ouoJ/ILVUyKJqrYHHmiY1bbL2IoukqtsUS0V2zaaitpL6+8kypbpjafY+iWN5dlFaQ6pTZMU03Gc5lqUkK6rIOKpMfisf48SR4s5c9tpX+5wqfpxWmWhVP1ntr4bxqsw0La6L6CT3OKqaZRugLt/QJjYRcUnxpecMawwA2o9YA4y3Uf8pCgAAYGCY2AAAgNZgYgMAAFqDiQ0AAGiNYdaKwkY2gJpQ9eYNMzWyFfzVSn2yn7AWmmbhJaoMlSpbqrzGs/dbkWk4saVhDalqjJv6f3SU9aOqmlgLecZVFNlCnk6ywpbzWBNHjuZjmZvN2zcVmWtH+x/QGPNFLacq42gqr0NVSl5/F3Gset0GkmnaVNNY2yTbsCHu2AAAgNZgYgMAAFqDiQ0AAGgNJjYAAKA1WDyMU9NgoXC5SLh6tHy1oLGJapEni4pxKgZw3VcLPKv2clFxtm3xeP+mi4TL92ayUDhbUHyyPly9N6vzdfhIf1vVd7VIuODqvCRj9FRe2qJ8PRsuNs5KU5QlFcqkiQYLzZVfo2VZjkFZwxjMHRsAANAaTGwAAEBrMLEBAACtwcQGAAC0BhMbAADQGmRF4aGqTIUqCyRbfd+0pEK1z6aZHcA6NIgsp8aSzJqTmiiyZYr3YJrpVL0vjxalBoqsGE/m+/S2bck+87IMZcZNVfbgwKF8+0xW2uEkfVfxrdMk66hpiYQq+6nKchtEBtQ6ykDljg0AAGgNJjYAAKA1mNgAAIDWYGIDAABaY2gTG9tX2t5t+7YVbW+yfZ/tW3r/Xjis8QBoJ2INMN6GmRV1laR3SvqjE9rfERG/PcRxYJCS1ffVyntPFiv15+bS9jjSXxumu8ukrkmnQdaWtK5W8GPgrtIgY02WYVJdV9m3V9lPTev5VO+f6en+xqre0vxCsc9i+2OL+eabkx8d1XEWWUReKMZS1ZzKMoOqbKGZoo9in2nGlaQ4lmRdHc1rP2X1s7qd5+e2vC6S7au+Y7FhvGpaW6pR300zt5LzMqD4O7Q7NhHxaUl7hrU/AOOJWAOMt/WwxubVtm/t3T4+Y9SDAdBaxBpgDIx6YvMHkp4o6UJJuyS9rdrQ9mW2d9reuajiFiYA5Ig1wJgY6cQmIh6IiOWI6Eh6t6SLTrLtFRGxIyJ2TGlmeIMEsOERa4DxMdKJje1zV3z6Ykm3VdsCwCNFrAHGx9Cyomz/qaRnSTrL9r2S3ijpWbYvlBSSvi7plcMaz9irsoKadpOt7C8yMtLsDUmeydsnipX6nUOHVzc4jKWBx5rkOixr7jSp6VPVVCveP5U0c2eqqGdUmNi2Ne97Mc+K0nJyTk7PM4uqGlLlmSprayXna8vmvO+lovbRTHEHrsg6ymJWWfuret2qa6LKXEvOV51pmrd3iuy3RtdnpcpcavozZQ0zUIc2sYmIS5Lm9wxr/wDGA7EGGG+jXjwMAAAwMExsAABAazCxAQAArcHEBgAAtMYwa0VhPRnUyvasi6rWy1TVXmRwFBkZaf9VLRVqSGENlNkl6XXVsFZUkdBT1lBqkgFVZXNVdY4mqnpGDWoLdZrVIaqyJFNlNlN+TsqaS6flWWHOMpSKuBQqsj6rsRwusjuX+y+ASLLQul/I28usvaUiyy3tpLrvUY1l/cRI7tgAAIDWYGIDAABag4kNAABoDSY2AACgNVg8PK4aLhIuHxeetLtYJOzZ2byPalHxIMo+VAs0gVPRyVf4lgvnM9Vj8qvrvnysfvL7abFtxQ1LMCiSfVYLXItz0vTdnZaOqBwrFslWZQxmitctWcirM07Lt31wX973/Hy+fbE413Nz/Y3FsUdxnOk1IUlRLCrOjrNMsKj6rla9Dx93bAAAQGswsQEAAK3BxAYAALQGExsAANAaTGwAAEBrkBWFh6pWvJeP105UGRmb8vbOtiQLQNLEQpEJkLVXGRlu8uh7ravHgmMdcP54+rKkQtZFkRVUPSbfp23JOyr2mZYJqPqusmWqx/7PL+T9ZFlUVeZk9d6ssiGrTJ/NSZxISh5IKmNNzM2k7Z3NRTmEqf5+lufyczW5Nc/6nNh7KO+7OC+x70AyjuIaKs5VrGUcKzIC1xPu2AAAgNZgYgMAAFqj0cTG9rTt/O8GADAgxBoAj9SqJza2XyNpl6S7bN9h+9VrNywA44pYA+BUPOzExvbv2r5U0mskfU9EnCfpxyU92favr3ZHti+wfYPt221/uRe8ZPtM29fbvrP3/xmP9GAAbFzEGgCDsJqsqBskPUPSWZI+a/uApFslfUnSq2y/PSL2rqKfJUmvi4ibbW+TdJPt6yX9oqRPRcRbbV8u6XJJr38Ex4LMIOotnaz7JGPE03mGQWdbnu3xzuvenbb/1Hv/W9r+hKvv6+/7gW/mAyyyn9ZRWRN82/qLNVHU0WmQJVhmPxUZSp1Dh9P2idOKGkVJHSHP5tk/Zb22KpNxe77P5dNW/1fCpa15FtGmfUXG1WL+5ozZ/h9Xnq+yovJzu+vHtud9Pye/rI78U//2m47kMXV5Ns9EmvnWtrT90Tfnxz+9u//cThw6km7rKiuqqC1VZuhlGXdrXWcv+9k0oGyuh313RsRHI+LXJH1O0oskPVfSVeoGjzMl/bXtr6yin10RcXPv44OS7pB0Xq/Pq3ubXS3p4sZHAWDDI9YAGIQmz7H5VUnXSLpF3d+gvkfSlyLiWbbzX9ELth8n6emSbpR0TkTs6n3pfknnNOkLQOsQawA8Yqu+nxoRd0r6QUnXSppV9xbxi3tfW3U9edtbJX1I0msj4iFPIoruU4XSe1G2L7O90/bOReW38ABsfMQaAKei0ZOHe0Hl//T+NWZ7St1A876I+HCv+QHb50bELtvnStpd7PsKSVdI0mk+k8fDAi1GrAHwSA3tAX22Lek9ku6IiLev+NJ1ki7tfXyppI8Na0wA2odYA4y3YdaK+lFJvyDpS7Zv6bW9QdJbJV1j+xWS7pb00iGOCSdwVe+las/q6CytPqtBkp44tTXv+3sPps3Lj+rPMpjYuy/fZ1VLRkVaVJVFRg2pjWR0sabKLmqwrbMsLElaKLKIkgwon55nM1W1ktTJM2CWzticti9u6890ik35e6czlbcvbsnjgZfz99r0vv6/Qi5tL7K/irfrgSfl5/aXn3Bz2v4jT7mzr+2Pv/kj6bY7778gbT981+lp++5n5GPfsqv/3J55cz5ubytiZ5FZV9WWyuL+mmeOrmFMHdrEJiI+I6nKPX7OsMYBoN2INcB4o1YUAABoDSY2AACgNZjYAACA1mBiAwAAWmOYWVEYhWrleZH9k9YMkeQqzSDLpujky+knv5VnOT3h2lem7du+lmeNTBx4INlnwxX2Va2fta6Pgo3FljclGUBV5lLaRbN6bVWttXL7zXnmUqrIitF0Xs/JS/n7YWZPf4bWsdOLcUf+XquynyaW8vaFs/qziDYdbpa6E5vz7Xds/mrafv9Sf0bT9qm8btN3nZU+Fkk7/znPUNv+1SJOLvSf887WvMbXxDf3pO2Vsj7Zsf6xZHUAJSmWikzTdZRRyh0bAADQGkxsAABAazCxAQAArcHEBgAAtAaLh7E6TRbVVgt5j86nzed+Jl90NnWoePz30f6Fi+XytGrcLBLGqkR6rZSlR7JFxTP5o/PLRZgTxWLbYoGvppIwXrzXtG1L2twpyp24eC8vF9tnpvfmBdk7s0VywHxRkiU75cXLsOlgUQS+ky/CvW7vM9L2p269p6/tP5396XTby+/9qbR9Zk9+nEfOzs/t6V9rsDD9tP7yMicTR47m/WQLhZvGyHVUdoY7NgAAoDWY2AAAgNZgYgMAAFqDiQ0AAGgNJjYAAKA1yIrCqUmyQOJYkZFQrLLf9rXDabsX8qyomE+yoqp9ViipgNWIvMxImRWVdVFkP03M5Rk6paKMQxzuf8S/N+WhPebysgdeKPqeyjN6OjP97VXVlcnD/e9XSVratjXveyp/b07v6c/0apKd1e08f90ePJZniy1G/3F+/NB3pds+Zdt9aftXfvBRafv+vzsnbd92TzLGKuOoKu1RXHPePJe2x/7+cjdVeR1N5NdEnYFKSQUAAIBHjIkNAABoDSY2AACgNZjYAACA1hjaxMb2BbZvsH277S/bfk2v/U2277N9S+/fC4c1JgDtQ6wBxtsws6KWJL0uIm62vU3STbav733tHRHx20McCwYklpOV8A1XwU/s78/qkCQv5JlOsdjfHsU+y5X9TVfwu0GmAkZtcLHGeQZUdV2lNXcKWXafpOa1orIMmKl824k9/dkvkrR03pn59vN5ds30A0km42SecbR0Wp79Nf3NvG5RJBlXktRJMrSiyKCaPyPP/po8LY8pF8ztTdt/dXt/rag/OnBWuu1EkRY2NZHHGlcloZLT6MVi42N55mglDuUZqE3qoZU1zrIYebL2NYyfQ5vYRMQuSbt6Hx+0fYek84a1fwDjgVgDjLeRrLGx/ThJT5d0Y6/p1bZvtX2l7TNGMSYA7UOsAcbP0Cc2trdK+pCk10bEAUl/IOmJki5U97estxXfd5ntnbZ3Lqq4hQsAPQOJNUGsATaaoU5sbE+pG2jeFxEflqSIeCAiliOiI+ndki7KvjciroiIHRGxY0ozwxs0gA1nYLHGxBpgoxlmVpQlvUfSHRHx9hXt567Y7MWSbhvWmAC0D7EGGG/DzIr6UUm/IOlLtm/ptb1B0iW2L5QUkr4u6ZVDHNP4KrN/is2rLJB0xX+e1VDWzDmYr9Sv6j+lmVhV9lNTI1jBj4EbXKwpakXV2yfZJZN5hk4UdX7caVivLBtf8V6r3oOb7t+X911kaHW2JTWHivfOZFGHauJYUc9o36G0femc0/v7Plpk6BQ6e/I7cFVG0xX7H9PXdtPBx6bbHlzMs7/2Hc3bN/WXvpIkdTYlWXibinsQVRZe0e65olZUki1VXZ8bIUYOMyvqM8p/bH5iWGMA0H7EGmC88eRhAADQGkxsAABAazCxAQAArcHEBgAAtMYws6KwnpQr25tlZESnf27sIsNAVe2RIvuplK3WL2s/NcwwAU6UZTpN55lO6bdX2U8NM1qiqAvkLVv6G2fy8Xkuz9CpxFT+I8KH+1N6nGUrnkSaWaUyMVOT+/trSy2etTXvo8hkmzgjjzUvPv2mtH1fZ3Nf242dJ6Tbfu1AXm9raWf+gOvZ+aK21KH+62XiUP6gyFgsroni+uzs3ZdvP9kfx8usqMpEcT03rcuXSWv1nWQoq+8ZAABgfWNiAwAAWoOJDQAAaA0mNgAAoDVYPDyuypIK1dK9Bl1XiyWjWDw8ny+Mc7FwMRaTx6i7mKNXC5lZVIzVSq6tcsF7su3E9FSz/VXvn5miIGcn2b5YyBt79+d9VGM8o7+MgVQszl3Kxx1b80XCEweO5Ns3OF/VguWjZ+eLpF/2vZ9N2x+zqX9hsiTdeviCvrZzZg6k23720OPT9rTqjKSpvHKEJrISFEXiRbVIOCuRcDJZmZpqcXssrWGixoDKNXDHBgAAtAYTGwAA0BpMbAAAQGswsQEAAK3BxAYAALQGWVF4qKbZUtnj5qtHxTfcZ/UI+Uar75uu1G+4+h5jILvGNxWZO8m2UV1TVWZV9f45mmfuaCL5/bR6v1aZVYv5WHw0z1hU8gh+dfL3mg8W2U9z+ViqfS6f1Z+h5SJDZ+7+/pIPkvT+W38gbX/pj38hbf/emXv72t511zPTbef358dzxt789XdxXXSm+19/Hype++x1OIkyi2ohOedFWYoyAzXLzpOaZdoOKP5yxwYAALQGExsAANAaTGwAAEBrMLEBAACtMbSJje1Z25+3/UXbX7b95l77423faPsu2x+0na9uAoBVINYA422YWVELkp4dEYdsT0n6jO0/l/SfJb0jIj5g+w8lvULSHwxxXFiNBtlSVa2ocjV9VvtJkhus+I9qBT9ZTuNozWNNWQ+tqOmTapo9mGU/Vf0s5e+pap9lBuKhoqDRVDInbJqBuFjss8jcmdzXP5bOlrwO1fLpea2oOJof/+0L56btH7j/or62qcki++tI3reL07J5d/4azexKalEVsTCO5tlf2lTU2cuyn6S0PlkdU9cw03Sj1YqKruNX5lTvX0h6tqRre+1XS7p4WGMC0D7EGmC8DXWNje1J27dI2i3peklfkbQvIo5PXe+VdN4wxwSgfYg1wPga6sQmIpYj4kJJ50u6SNJ3r/Z7bV9me6ftnYsqbqcBgIg1wDgbSVZUROyTdIOkH5a03fbxPwieL+m+4nuuiIgdEbFjSsXTMwFgBWINMH6GmRV1tu3tvY/nJD1P0h3qBp2X9Da7VNLHhjUmAO1DrAHG2zCzos6VdLXtSXUnVNdExMdt3y7pA7Z/Q9I/SHrPEMeEYapW00e+Ej6K5Ko0I6XKuBrQKntsKIONNdW1lcgySVwlP1VZJ5WqnyS7KIqEoyorylUNqSozMastVWTiVOfP1fZVZlmWuXX6lnTTznS+z4kiK+oNN744bY8j/WOcmM/7ntmXt2+5P89+Wi7GGDP9dcjKml1F9lvM59tX5zyNqZ0is658P1TxvcF1PqC4PLSJTUTcKunpSftX1f0bOACcMmINMN548jAAAGgNJjYAAKA1mNgAAIDWYGIDAABaY5hZUWijQdQBKftuUJOkzLgi+wmnKLu2qsyQbNsB1dyJJjXVirpSnipCfpWhVNUWmu7P3Cl1ivTG5YaZjMkxeSFP/9p0MO/7tDvzGlL7nppnS83d039eppNSTpK0/a58LFOHiuyiIjZNJDWxqnp6TXWOHFn1ttW1Ul2HpYkinS+7/gcUr7ljAwAAWoOJDQAAaA0mNgAAoDWY2AAAgNZg8TAGr+Ei4fRx3lKjR9k3RqkFrFZyHbp67H9S96C6vieKBbhRXYPVIuSsTELRR7kAuWl5h6y8QbHY1NPTxVjyxbau3pvJcXo+Ke0gaaLoY8v9RaxR/lqcdnf/+epM531PHssXg1flHaYO5GPXUv8YY34+37aKnUW75+bS9s7hZFFxVVKhqWrxeNNkkga4YwMAAFqDiQ0AAGgNJjYAAKA1mNgAAIDWYGIDAABag6woDF7TzKJBrI6v9kn2E05V8uj36DQoqVBk98Vys5IKaekESbHUn73iLFPqJLI+uh3l759s+4nqEfxHj+Z9F2OMJCtIkrKRxJY8y0fFuZ3bnWciTR3Kx7K4tb999ltFGYf9efmJmMpft4n9eXmDLFuszCybz/fp2Zlm25dZfkkfjbNYh1/uhjs2AACgNZjYAACA1mBiAwAAWoOJDQAAaI2hTWxsz9r+vO0v2v6y7Tf32q+y/TXbt/T+XTisMQFoH2INMN6GmRW1IOnZEXHI9pSkz9j+897X/mtEXDvEsWAjyDJMGvdB9tMYGlmsaZSNVF3fVRZVcS1ntZXKzJWiJlSVcaWJIqNnbnbVfTcVR/JsIZ+2rb9t/8G8j21b0vbJo3lG0+LW/Efh1rsP9zcuF+dwMT/nk/sOpe1lPa+srchmqlQ1wcpMvOy6bZr9VF3PI4jBQ5vYRPddefwVnur946cOgIEi1gDjbahrbGxP2r5F0m5J10fEjb0vvcX2rbbfYTtPwAeAVSLWAONrqBObiFiOiAslnS/pItvfJ+m/S/puST8g6UxJr8++1/Zltnfa3rmoZrflAIwXYg0wvkaSFRUR+yTdIOn5EbEruhYkvVfSRcX3XBEROyJix5T4RQvAwyPWAONnmFlRZ9ve3vt4TtLzJP2j7XN7bZZ0saTbhjUmAO1DrAHG2zCzos6VdLXtSXUnVNdExMdt/7Xts9UtBXKLpFcNcUxYD0ZRWwpttq5iTTTMFnKVWFVkqUSWpdI4o7CqZ1Vkxiw2eA8WmVWuMrSqmlPHkjpPE8XJ2pdnS1X1rDZ/dT7vJ6mh1JnN6zZNHCn6OJZnYkXR3uS1i4XiT6WdKkOpaE9e5qbXbWNZHB9QBtUws6JulfT0pP3ZwxoDgPYj1gDjjScPAwCA1mBiAwAAWoOJDQAAaA0mNgAAoDWGmRUFDAb1nzBqRXZJZElERW2dqq5ULOXZMt40lY8lydyxG4b2YiyusqIy1ftyscgKKrKlKp5Jtu/kNZG0KT9+7y/qNhXZUlrqP/7JKvupqsNUHedSMfbk2qpqP7k4zk6VLVXJXv/q3DatCTWCLFbu2AAAgNZgYgMAAFqDiQ0AAGgNJjYAAKA1WDwMAINSLBTOVOUKykXCTR6HX+20Gl+58DVf+NloUXGxeLZcEDubFx6NZLFttXg25osFvotJWQZJnpvLt08W1UYxPu07kLdXpQkmi/NyNB9j2nWxSNjFgt20/IaK16K63sprvOGi4rTvYqFxw4QR7tgAAIDWYGIDAABag4kNAABoDSY2AACgNZjYAACA1iArCgBOplFGRpIZ0iBTSqqzpSpOMpeiyMRxXjmhUcaVJEVRgiHdZyfv21UZg+r4s4yZatzFPrW5yH6qJK+99x/MNz1aZGIVxxnzRdmDZJ9Nr4ksg0yqr4tUdd12Grw+TQ2oXA53bAAAQGswsQEAAK3BxAYAALQGExsAANAaTGwAAEBrOAa0CnmYbH9T0t29T8+S9OAIhzMsHGe7tPk4HxsRZ496EIOwIta0+fVaieNslzYfZxlnNuTEZiXbOyNix6jHsdY4znYZl+Nsi3F5vTjOdhmX4zwRf4oCAACtwcQGAAC0RhsmNleMegBDwnG2y7gcZ1uMy+vFcbbLuBznQ2z4NTYAAADHteGODQAAgKQNPLGx/Xzb/2T7LtuXj3o8g2T7Stu7bd+2ou1M29fbvrP3/xmjHOOpsn2B7Rts3277y7Zf02tv23HO2v687S/2jvPNvfbH276xd/1+0Pb0qMeKXFtjzTjEGYlYM46xZkNObGxPSvp9SS+Q9GRJl9h+8mhHNVBXSXr+CW2XS/pURDxJ0qd6n29kS5JeFxFPlvRDkn619xq27TgXJD07Ip4m6UJJz7f9Q5J+S9I7IuI7Je2V9IrRDRGVlseaq9T+OCMRa8Yu1mzIiY2kiyTdFRFfjYhjkj4g6UUjHtPARMSnJe05oflFkq7ufXy1pIuHOaZBi4hdEXFz7+ODku6QdJ7ad5wREYd6n071/oWkZ0u6tte+4Y+zxVoba8YhzkjEGo1hrNmoE5vzJN2z4vN7e21tdk5E7Op9fL+kc0Y5mEGy/ThJT5d0o1p4nLYnbd8iabek6yV9RdK+iFjqbTIO1+9GNW6xpnXvv5WINa2/fiVt3InNWItuKlsr0tlsb5X0IUmvjYgDK7/WluOMiOWIuFDS+ereAfju0Y4IeHhtef8dR6wZHxt1YnOfpAtWfH5+r63NHrB9riT1/t894vGcMttT6gaa90XEh3vNrTvO4yJin6QbJP2wpO22N/W+NA7X70Y1brGmle8/Ys14xZqNOrH5gqQn9VZ7T0t6maTrRjymtXadpEt7H18q6WMjHMsps21J75F0R0S8fcWX2nacZ9ve3vt4TtLz1P0b/w2SXtLbbMMfZ4uNW6xp1ftPItZoDGPNhn1An+0XSvodSZOSroyIt4x2RINj+08lPUvdyqwPSHqjpI9KukbSd6hbbfilEXHiwr8Nw/YzJf2tpC9J6vSa36Du377bdJxPVXfB3qS6v0hcExH/0/YT1F2Ieqakf5D08xGxMLqRotLWWDMOcUYi1oxjrNmwExsAAIATbdQ/RQEAAPRhYgMAAFqDiQ0AAGgNJjYAAKA1mNgAAIDWYGIDAABag4kNAABoDSY2GBnbr7f9x7b/0vbdtv/LqMcEoH2INeOFiQ1G6SnqPgn0hZKeK+nlox0OgJYi1oyRTQ+/CbBmniLppyNi2faypA37OHMA6xqxZoxwxwYj0au2e3ZEfKXX9FR1a7kAwMAQa8YPExuMynepW3n2uAslfXE0QwHQYsSaMcPEBqPyFEm3rvj8QhFsAAwesWbMUN0b64LtuyQ9JSKOjnosANqLWNN+3LHByNneLukYgQbAWiLWjAfu2AAAgNbgjg0AAGgNJjYAAKA1mNgAAIDWYGIDAABag4kNAABoDSY2AACgNZjYAACA1mBiAwAAWuP/A3Ppw6/8R+3dAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(np.sum(qcd_images, axis=0))\n",
    "plt.title(\"Average QCD\")\n",
    "plt.xlabel(\"$\\eta$\")\n",
    "plt.ylabel(\"$\\phi$\")\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(np.sum(top_images, axis=0))\n",
    "plt.title(\"Average top\")\n",
    "plt.xlabel(\"$\\eta$\")\n",
    "plt.ylabel(\"$\\phi$\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "756c0ebe-715b-4bbf-b37c-59472ef5d965",
   "metadata": {},
   "source": [
    "Great, we now have a way to convert raw 4-vectors into jet images! For convenience, we've applied the above steps to the whole `top_tagging` dataset, so let's download it now from the [Hugging Face Hub](https://huggingface.co/datasets/dl4phys/top_tagging_images):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "92093aca-a879-4633-96cd-9a421b2b3606",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "351d1b9ce7b046498587b888d78c9a8c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/3 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "DatasetDict({\n",
       "    validation: Dataset({\n",
       "        features: ['image', 'label'],\n",
       "        num_rows: 403000\n",
       "    })\n",
       "    test: Dataset({\n",
       "        features: ['image', 'label'],\n",
       "        num_rows: 404000\n",
       "    })\n",
       "    train: Dataset({\n",
       "        features: ['image', 'label'],\n",
       "        num_rows: 1211000\n",
       "    })\n",
       "})"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "images_ds = load_dataset(\"dl4phys/top_tagging_images\")\n",
    "images_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5768bf5-a691-4244-a1a2-c461469c8c66",
   "metadata": {},
   "source": [
    "Here we now have three splits to work with, and we can access an individual image by indexing as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7c777e03-142f-4630-a682-5f63031a6cf4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD6CAYAAABnLjEDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAMIklEQVR4nO3df6jdd33H8efL26TpqiUG26xrstm5MilDM+mqMhlduo6s/6SCiIWNDAp1YwVlMsz8Rx0TOqZ2f2w4FLNmw1lL1SVINw1ZwAmSttYY00bX2FVMSBOdy5oOFpfkvT/ON3Kbe25yen7dc+/n+YDD+Z7P93vu9/1t7qvfcz73nO87VYWkle8VS12ApOkw7FIjDLvUCMMuNcKwS40w7FIjRgp7ki1JvpvkSJLt4ypK0vhl2L+zJ5kD/h24AzgKPA7cXVVPL/ac1bmy1nD1UPuTdHn/y//wkzqTfuuuGOHn3gocqapnAZI8BGwFFg37Gq7mzbl9hF1KupT9tXfRdaO8jL8B+MG8x0e7MUkzaJQz+0CS3AvcC7CGn5n07iQtYpQz+zFg47zHG7qxl6iqT1bVLVV1yyquHGF3kkYxStgfB25KcmOS1cC7gN3jKUvSuA39Mr6qzia5D/gyMAfsqKqnxlaZpLEa6T17VT0KPDqmWiRNkJ+gkxph2KVGGHapEYZdaoRhlxph2KVGGHapEYZdaoRhlxph2KVGGHapEYZdaoRhlxph2KVGGHapEYZdaoRhlxph2KVGjHRZqiTPAaeBc8DZqrplHEVJGr9xXDf+N6vqR2P4OZImyJfxUiNGDXsBX0nyja7zi6QZNerL+LdV1bEk1wF7knynqr46fwPbP0mzYaQze1Ud6+5PAl+k19n14m1s/yTNgKHDnuTqJK+6sAz8NnBoXIVJGq9RXsavB76Y5MLP+ceq+pexVCVp7Ebp9fYs8MYx1iJpgvzTm9QIwy41wrBLjTDsUiMMu9QIwy41wrBLjTDsUiMMu9QIwy41wrBLjTDsUiMMu9QIwy41wrBLjTDsUiMMu9QIwy41wrBLjbhs2JPsSHIyyaF5Y+uS7EnyTHf/6smWKWlUg5zZHwS2XDS2HdhbVTcBe7vHkmbYZcPedXj58UXDW4Gd3fJO4K7xliVp3Ia9lPT6qjreLT9P7xryfdn+SZoNI0/QVVXRa/C42HrbP0kzYNiwn0hyPUB3f3J8JUmahGHDvhvY1i1vA3aNpxxJkzLIn94+C3wd+OUkR5PcA9wP3JHkGeC3useSZthlJ+iq6u5FVt0+5lokTZCfoJMaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaMexlqbQCveHJLBg7+KZFL0KkZcYzu9QIwy41wrBLjTDsUiOGbf/0oSTHkhzobndOtkxJoxpkNv5B4K+Bv79o/IGq+ujYK9KSceZ9ZRu2/ZOkZWaU9+z3JTnYvcy3i6s044YN+yeA1wGbgOPAxxbbMMm9SZ5I8sT/cWbI3Uka1VBhr6oTVXWuqs4DnwJuvcS29nqTZsBQYb/Q563zduDQYttKmg2XnY3v2j/dBrwmyVHgg8BtSTbR6976HPDuyZUoaRyGbf/06QnUImmC/ASd1AjDLjXCsEuN8OIVurRXzC0cO39u+nVoZJ7ZpUYYdqkRhl1qhGGXGuEEnS7NybgVwzO71AjDLjXCsEuNMOxSIwy71Ahn4zU2c9dcs2Ds3AsvLEEl6sczu9QIwy41wrBLjRik/dPGJPuSPJ3kqSTv6cbXJdmT5Jnu3mvHSzNskAm6s8D7qurJJK8CvpFkD/D7wN6quj/JdmA78P7JlapJu2LjhgVjN+86tmBssTZRTsbNtkHaPx2vqie75dPAYeAGYCuws9tsJ3DXhGqUNAYv6z17ktcCvwrsB9ZX1fFu1fPA+vGWJmmcBg57klcCnwfeW1Uveb1WVUXvGvL9nmf7J2kGDBT2JKvoBf0zVfWFbvjEhc4w3f3Jfs+1/ZM0GwbpCBN6TSEOV9XH563aDWwD7u/ud02kQk3N6Tf93IKxv/zZLy0Y+42t9/Z9/lW7Hht4X37abvoGmY3/deD3gG8nOdCNfYBeyB9Ocg/wfeCdE6lQ0lgM0v7pa0AWWX37eMuRNCl+gk5qhGGXGmHYpUb4fXb9VL/Z9Ddf84cLxtbu+nrf58+tv27B2Lkf/mffbZ15nz7P7FIjDLvUCMMuNcKwS41YuRN0/fqKg+2MgFzZ/zsKdWbhF5XW/kP/ybh+zp/67z6DU/7vbT/5RXlmlxph2KVGGHapEYZdaoRhlxqxcmfjnYFdVL9Z91n+uS+L/+6L8swuNcKwS40w7FIjRmn/9KEkx5Ic6G53Tr5cScMapf0TwANV9dHJlSdpXAa54ORx4Hi3fDrJhfZPkpaRUdo/AdyX5GCSHXZxlWbbKO2fPgG8DthE78z/sUWeZ/snaQYM3f6pqk5U1bmqOg98Cri133Nt/yTNhkFm4/u2f7rQ563zduDQ+MuTNC6jtH+6O8kmet1bnwPePYH6xs+LG6hRo7R/enT85UiaFD9BJzXCsEuNMOxSI1bu99kX42ScGuWZXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaMcgFJ9ckeSzJt7r2Tx/uxm9Msj/JkSSfS7J68uVKGtYg32c/A2yuqhe7S0p/Lck/A39Mr/3TQ0n+FriH3rXktQTm1l/Xd/zciZNTrkSz6rJn9up5sXu4qrsVsBl4pBvfCdw1iQIljcegTSLmustInwT2AN8DTlXV2W6To9j/TZppA4W96/yyCdhAr/PL6wfdge2fpNnwsmbjq+oUsA94K7A2yYX3/BuAY4s8x/ZP0gwYZDb+2iRru+WrgDuAw/RC/45us23ArgnVKGkMBpmNvx7YmWSO3v8cHq6qLyV5GngoyZ8D36TXD05LxFl3Xc4g7Z8O0uvJfvH4syzSuVXS7PETdFIjDLvUCMMuNcKwS40w7FIjDLvUCMMuNcKwS40w7FIjDLvUCMMuNcKwS40w7FIjDLvUCMMuNcKwS40w7FIjDLvUiFHaPz2Y5D+SHOhumyZeraShjdL+CeBPquqRSzxX0owY5IKTBfRr/yRpGRmq/VNV7e9WfSTJwSQPJLEDhDTDhmr/lORXgD+l1wbq14B1wPv7Pdf2T9JsGLb905aqOt51eD0D/B2LXEPe9k/SbBi2/dN3klzfjYVeu+ZDkytT0qhGaf/0r0muBQIcAP5gcmVKGtUo7Z82T6QiSRPhJ+ikRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcaYdilRhh2qRGGXWqEYZcakV53pyntLPkh8P3u4WuAH01t59PjcS0/K+nYfqGqru23Yqphf8mOkyeq6pYl2fkEeVzLz0o+tvl8GS81wrBLjVjKsH9yCfc9SR7X8rOSj+2nluw9u6Tp8mW81Iiphz3JliTfTXIkyfZp73+ckuxIcjLJoXlj65LsSfJMd//qpaxxGEk2JtmX5OkkTyV5Tze+rI8tyZokjyX5VndcH+7Gb0yyv/ud/FyS1Utd6yRMNexdJ9i/AX4HuBm4O8nN06xhzB4Etlw0th3YW1U3AXu7x8vNWeB9VXUz8Bbgj7p/p+V+bGeAzVX1RmATsCXJW4C/AB6oql8C/gu4Z+lKnJxpn9lvBY5U1bNV9RPgIWDrlGsYm6r6KvDji4a3Aju75Z30etcvK1V1vKqe7JZPA4eBG1jmx1Y9L3YPV3W3AjYDj3Tjy+64BjXtsN8A/GDe46Pd2EqyvqqOd8vPA+uXsphRJXktvZbd+1kBx5ZkLskB4CSwB/gecKqqznabrMTfScAJuomq3p86lu2fO5K8Evg88N6qemH+uuV6bFV1rqo2ARvovdJ8/dJWND3TDvsxYOO8xxu6sZXkRJLrAbr7k0tcz1CSrKIX9M9U1Re64RVxbABVdQrYB7wVWJvkim7VSvydBKYf9seBm7rZz9XAu4DdU65h0nYD27rlbcCuJaxlKEkCfBo4XFUfn7dqWR9bkmuTrO2WrwLuoDcfsQ94R7fZsjuuQU39QzVJ7gT+CpgDdlTVR6ZawBgl+SxwG71vTZ0APgj8E/Aw8PP0vuH3zqq6eBJvpiV5G/BvwLeB893wB+i9b1+2x5bkDfQm4Obonegerqo/S/KL9CaL1wHfBH63qs4sXaWT4SfopEY4QSc1wrBLjTDsUiMMu9QIwy41wrBLjTDsUiMMu9SI/wcVd3naJg9pBAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "img = images_ds[\"train\"][0][\"image\"]\n",
    "plt.imshow(img);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26b567bf-881f-46c2-b906-8c981172c912",
   "metadata": {},
   "source": [
    "Okay, this looks like a jet image - let's now take a look at converting them to PyTorch tensors that are suitable for training a neural net!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e63935f-1d9a-4014-a829-e937dffffac8",
   "metadata": {},
   "source": [
    "## From images to tensors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d80fc93c-9254-4258-9c8c-2c445302f9f1",
   "metadata": {},
   "source": [
    "Now that we have a dataset of images, the next thing we need to do is convert them to PyTorch tensors and wrap them in a `DataLoader`. One way to do this is via `torchvision`'s `ToTensor` class, which converts a PIL image into a tensor as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "df6dd8e8-9d67-409d-bad1-ad9b3ed16914",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([1, 40, 40])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = ToTensor()(img)\n",
    "t.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3036062-e6f7-423e-9f1d-04dfcd9c07c1",
   "metadata": {},
   "source": [
    "Here we can see that we now have a tensor of shape `(num_channels, height, width)`, where `num_channels` is 1 because we're dealing with black and white images, and `height` and `width` denote the number of pixels. To apply `ToTensor` to multiple images, we can use a list comprehension and `torch.cat()` to concatenate all the tensors as a single tensor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "d3e25468-ab7f-4a05-8d7d-09efbe32bd21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([10, 40, 40])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = torch.cat([ToTensor()(img) for img in images_ds[\"train\"][:10][\"image\"]])\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "189ae367-45ba-4f5f-84dd-442144066fe9",
   "metadata": {},
   "source": [
    "Now this isn't quite what we want because CNNs in PyTorch expect rank-4 tensors with a shape `(minibatch, channels, height, width)`, where `channels` refers to the number of color channels (e.g. RGB) of the image. In our case, the images are black and white, so we just need to reshape the tensor to insert a new dimension in the second position:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "72aad838-5c44-4990-9cde-4e5e8ab00e74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([10, 1, 40, 40])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = torch.cat([ToTensor()(img) for img in images_ds[\"train\"][:10][\"image\"]])\n",
    "x.unsqueeze_(dim=1)\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "400612d5-1f28-41d5-b28d-916f9ac2520f",
   "metadata": {},
   "source": [
    "Great, we now have the image tensors in the right shape, so let write a simple `get_dataset()` helper function that applies this over a single training split:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "f1381d63-3d56-40ea-a794-6cd291967035",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dataset(dataset, num_examples=None):\n",
    "    if num_examples is not None:\n",
    "        dataset = dataset.shuffle(seed=42).select(range(num_examples))\n",
    "\n",
    "    x = torch.cat([ToTensor()(img) for img in dataset[\"image\"]]).unsqueeze(1)\n",
    "    y = torch.cat([torch.tensor(l).unsqueeze(0) for l in dataset[\"label\"]])\n",
    "\n",
    "    return TensorDataset(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3d5ec9-086e-4f72-90fe-e9c77082a997",
   "metadata": {},
   "source": [
    "Here we've also collected the labels in a single tensor and returned a `TensorDataset` object that we can iterate over:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "e2d1950a-c482-47f6-85d1-6b5686d636bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lower the training size if Colab RAM explodes 💣\n",
    "train_ds = get_dataset(images_ds[\"train\"], num_examples=350_000)\n",
    "valid_ds = get_dataset(images_ds[\"validation\"], num_examples=50_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82922865-3e13-4310-b175-856d63b5c98b",
   "metadata": {},
   "source": [
    "Now that we have the PyTorch datasets prepared, let's wrap them in dataloaders so we can generate minibatches during training:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "c134bd01-07bc-4e0e-976e-19444e549177",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_dl = DataLoader(train_ds, bs=128, shuffle=True)\n",
    "valid_dl = DataLoader(valid_ds, bs=128)\n",
    "dls = DataLoaders(train_dl, valid_dl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "598e674b-015c-4e78-9807-4fde32858eb3",
   "metadata": {},
   "source": [
    "## Fine-tuning a pretrained CNN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "178819c0-82f1-4761-867c-2dba249a318e",
   "metadata": {},
   "source": [
    "Now that we have prepared our data and created dataloaders, we can train a CNN! Instead of training a CNN from scratch, we'll use a _pretrained model_ that has been already been trained to classify over a million photos in the famous [ImageNet dataset](https://www.image-net.org/). These models are called \"pretrained\" because their weights have already been optimized by training on another dataset. The layers in a pretrained model are well suited for detecting edges, colour gradient and so on, which allows us to train a new model with much less labelled data, and much faster. \n",
    "\n",
    "In practice, we replace the last layer of the pretrained model with a randomly initialized one that has the appropriate shape for the task at hand. For instance, ImageNet involves 1,000 different classes, so if we want to re-use the weights of a model trained on this dataset for a new set of classes, we'll need to resize the last layer. This last layer is often referred to as the _head_ of the network, with every layer before it belonging to the _body_.\n",
    "\n",
    "This process of reusing the weights to solve a different task is called _transfer learning_, and now underpins much of the success that deep learning has achieved in computer vision, natural language processing, and audio. A cartoon of the process is shown below:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b34afde3-f713-48ce-b880-d579676aef16",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img alt=\"transfer\" caption=\"Transfer learning\" src=\"images/transfer-learning.png\" id=\"transfer-learning\" width=800/>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28f97df0-8ee1-4546-b80a-c1045a33560d",
   "metadata": {},
   "source": [
    "Anyway, enough jargon - let's train a model! As we did in lecture 1, we'll use the high-level API of fastai, which provides factory methods for various domains. Instead of a `tabular_learner()` we'll use the `vision_learner()`, which is designed for computer vision tasks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "1b47b804-1b6f-48d3-8d43-a7a3831810c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "learn = vision_learner(\n",
    "    dls,\n",
    "    resnet34,\n",
    "    metrics=[accuracy, RocAucBinary()],\n",
    "    n_in=1,\n",
    "    n_out=2,\n",
    "    normalize=False,\n",
    "    loss_func=F.cross_entropy,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bc8c472-2d65-4a21-902c-6677acf9b6de",
   "metadata": {},
   "source": [
    "Here we've provided the dataloaders and specified `resnet34` as the name of the pretrained model we wish to initialise with. ResNets are a popular type of CNN that we might cover in a future lesson, but for now it is enough to know that they are great for most computer vision tasks. Now, because we're using custom dataloaders instead of the factory methods provided by fastai, we need to specify that our model should run on a GPU instead of CPU (the default). We can do that by tapping into the model and selecting the desired device name with PyTorch's `to()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "bee77d68-8e49-44d3-86e5-c354819afe30",
   "metadata": {},
   "outputs": [],
   "source": [
    "learn.model.to(\"cuda\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5213ae1c-0e58-40d5-82ce-fa46cd570798",
   "metadata": {},
   "source": [
    "Now that our model is on the GPU, we can follow the same steps we've done in previous lectures. First we find a good learning rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "cd54e49f-a10c-4a91-9554-3734153c3ca8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "SuggestedLRs(valley=0.0010000000474974513)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAENCAYAAAD6/JlzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsXklEQVR4nO3deXiU5b3/8fc3G9kggRAWCfu+ChjBHRQVxbW2bqUurdZarZ5q67H9tafSVk9PT09r61b3pWpBS13RKtqCYBE1LLKD7IQ1QBJIJstkcv/+mAkGyGSBTJ6Z5PO6rlwkzzwzzydxzDf38ty3OecQERGpS5zXAUREJHqpSIiISFgqEiIiEpaKhIiIhKUiISIiYalIiIhIWBErEmb2rJntMbMVYR7PMLO3zewLM1tpZt+OVBYRETk2kWxJPA9cUM/jtwOrnHMnAhOB35tZUgTziIhIE0WsSDjn5gH76zsFaG9mBqSHzq2KVB4REWm6BA+v/QjwFrADaA9c7ZyrbuhJnTt3dn369IlwNBGR1mXRokV7nXPZTX2el0ViMrAUOAfoD3xgZvOdcweOPNHMbgFuAejVqxd5eXktmVNEJOaZ2ZZjeZ6Xs5u+DbzmgtYDm4AhdZ3onHvSOZfrnMvNzm5yIRQRkWPkZZHYCkwCMLOuwGBgo4d5RETkCBHrbjKz6QRnLXU2s3zgPiARwDn3OPBr4HkzWw4YcK9zbm+k8oiISNNFrEg4565t4PEdwPnNcS2/309+fj7l5eXN8XIxKTk5mZycHBITE72OIiKtiJcD180mPz+f9u3b06dPH4IzatsW5xz79u0jPz+fvn37eh1HRFqRVrEsR3l5OVlZWW2yQACYGVlZWW26JSUikdEqigTQZgtEjbb+/Yu0dh+s2s36PQdb/LqtpkjEkvT0dAA2b97MiBEjPE4jIrHg9pcX8/fF21v8um2zSCx7FR4cAdMyg/8ue9XrRCIiYVUFqqkMVJOaGN/i1257RWLZq/D2nVC8DXDBf9++87gKxU9+8hMeffTRQ19PmzaN+++/n0mTJjF27FhGjhzJm2++We9rBAIB7rnnHk4++WRGjRrFE088AcD111/PG2+8cei8qVOnNvhaItK6+PwBAFKSVCQi75+/An/Z4cf8ZcHjx+jqq6/m1Ve/KjKvvvoqN9xwA6+//jqLFy9mzpw5/OhHP8I5F/Y1nnnmGTIyMvj888/5/PPPeeqpp9i0aRM33XQTzz//PADFxcUsWLCAiy666JizikjsKasMFonUpJafkNoqpsA2SXF+0443wpgxY9izZw87duygoKCAjh070q1bN+666y7mzZtHXFwc27dvZ/fu3XTr1q3O15g9ezbLli1j5syZwTjFxXz55Zecf/753HbbbRQUFPD3v/+dr3/96yQktL3/bCJtme9QkWj5lkTb+22TkRPqaqrj+HG48sormTlzJrt27eLqq6/m5ZdfpqCggEWLFpGYmEifPn3qnaLqnOPhhx9m8uTJRz12/fXX89JLLzFjxgyee+6548opIrHHVxncRUHdTS1h0i8gMeXwY4kpwePH4eqrr2bGjBnMnDmTK6+8kuLiYrp06UJiYiJz5sxhy5b6F2CcPHkyf/7zn/H7/QCsW7eO0tJSAG688Ub++Mc/AjBs2LDjyikisUctiZY06qrgv//8VbCLKSMnWCBqjh+j4cOHc/DgQXr06EH37t2ZOnUql1xyCSNHjiQ3N5chQ+pc4PaQm2++mc2bNzN27Ficc2RnZx8asO7atStDhw7l8ssvP66MIhKbvCwSVt9gajTKzc11R+4nsXr1aoYOHepRosjz+XyMHDmSxYsXk5GREfa81v5zEGmr3luxk1tfWsy7d57JsBM6HNNrmNki51xuU5/X9rqbYsyHH37I0KFDueOOO+otECLSeqm7ScI699xzGxzPEJHWzcsioZaEiEiUq7lPQrObjkOsja00t7b+/Yu0Zj4Pb6ZrFUUiOTmZffv2tdlflDX7SSQnJ3sdRUQiwOevIikhjvi4ll/tuVWMSeTk5JCfn09BQYHXUTxTszOdiLQ+ZZUBT8YjoJUUicTERO3IJiKtlq8y4MkKsNBKuptERFqzssqAJ4PWoCIhIhL1fJVVngxag4qEiEjU86klISIi4ZT5vRu4VpEQEYlyPg9nN6lIiIhEubLKACmJGpMQEZE6BAeu1ZIQEZE6qLtJRETqFKh2VFRVa3aTiIgcrczv3TLhoCIhIhLVfJVVAKToZjoRETlSzV4SaWpJiIjIkbzclQ5UJEREopq6m0REJCy1JEREJKyaIpGi/SRERORIZWpJiIhIOF91N2lMQkREjvDVwLVaEiIicoRW291kZs+a2R4zW1HPORPNbKmZrTSzjyKVRUQkVvn8ARLjjcR4b/6mj+RVnwcuCPegmWUCjwGXOueGA1dGMIuISEwK7iXhTSsCIlgknHPzgP31nPJN4DXn3NbQ+XsilUVEJFYF95LwZtAavB2TGAR0NLO5ZrbIzK4Pd6KZ3WJmeWaWV1BQ0IIRRUS85eVeEuBtkUgATgIuAiYD/2Vmg+o60Tn3pHMu1zmXm52d3ZIZRUQ8VVYZ8GxmEwR/UXslH9jnnCsFSs1sHnAisM7DTCIiUaUttyTeBM4wswQzSwXGA6s9zCMiEnV8/oBni/tBBFsSZjYdmAh0NrN84D4gEcA597hzbrWZvQcsA6qBp51zYafLioi0RWWVVXTvkOzZ9SNWJJxz1zbinN8Bv4tUBhGRWNeWu5tERKQBXg9cq0iIiEQxtSRERKRO1dWOMo8HrlUkRESiVHmVt4v7gYqEiEjU8nrrUlCREBGJWmUeb10KKhIiIlHL613pQEVCRCRqlYZ2pVN3k4iIHOVQd5OKhIiIHKmmuylN3U0iInIkX6i7SS0JERE5SpmmwIqISDi6T0JERMIq82vgWkREwvBVVhEfZyTFe/erWkVCRCRK+SoDpCbGY2aeZVCREBGJUl7vJQEqEiIiUcvrvSRARUJEJGr5Kr3dSwJUJEREolaZv0otCRERqZu6m0REJKyyyoCne0mAioSISNRSS0JERMLSwLWIiIRVVqmBaxERqYNzDp9f3U0iIlKHiqpqnPN2cT9QkRARiUqHlgnX7CYRETlSza50qRq4FhGRI9XsSqfuJhEROUppFOxKByoSIiJRqaa7SS0JERE5StmhloTGJERE5Ag+dTeJiEg4hwauNQVWRESO9NUUWBUJERE5gs+vMQkREQmjrDKAGSQnevtrOmJXN7NnzWyPma1o4LyTzazKzL4RqSwiIrHGVxkgNTEeM/M0RyRL1PPABfWdYGbxwG+B2RHMISISc4p8ftLaedvVBBEsEs65ecD+Bk67A/g7sCdSOUREYtGqnQcY3K291zG8G5Mwsx7A14A/e5VBRCQalVUGWLf7ICfmZHodxdOB6z8C9zrnqhs60cxuMbM8M8srKCiIfDIREQ+t2nmAQLVjZE6G11HwssMrF5gRGpTpDEwxsyrn3BtHnuicexJ4EiA3N9e1ZEgRkZa2PL8IICpaEp4VCedc35rPzex5YFZdBUJEpK1Zll9Mdvt2dO3QzusokSsSZjYdmAh0NrN84D4gEcA593ikrisiEuuWbS9mVI8Mz6e/QiOLhJmlAWXOuWozGwQMAf7hnPOHe45z7trGhnDO3djYc0VEWrOSiio2FJRwyagTvI4CNH7geh6QHJqRNBu4juB9ECIi0oxWbC/GORgVBYPW0PgiYc45H3AF8Jhz7kpgeORiiYi0TcvziwGiYmYTNKFImNmpwFTgndAxb5cmFBFphZZtL6ZHZgqd070ftIbGF4kfAj8FXnfOrTSzfsCciKUSEWmjluUXMbJHdLQioJED1865j4CPAMwsDtjrnLszksFERNqaYp+fLft8XJXb0+sohzSqJWFmfzWzDqFZTiuAVWZ2T2SjiYi0Lcu3B8cjouEmuhqN7W4a5pw7AFwO/APoS3CGk4iINJMvQndaR1N3U2OLRKKZJRIsEm+F7o/Q8hgiIs1oeX4xvbNSyUhN9DrKIY0tEk8Am4E0YJ6Z9QYORCqUiEhbtHx7MaOiqKsJGlkknHMPOed6OOemuKAtwNkRziYi0mbsLalge1EZo6KoqwkaP3CdYWZ/qFmu28x+T7BVISIizWDDnhIABkXBRkO1Nba76VngIHBV6OMA8FykQkXCsvwi7n51KQfLwy43JSLimUJf8HdTVlqSx0kO19gi0d85d59zbmPo45dAv0gGa26FPj+vLd5+6JZ3EZFoUlxWCUBmFA1aQ+OLRJmZnVHzhZmdDpRFJlJkjA4NBi3ZVuRpDhGRutS0JDqmRldLorH7SdwK/MXMakZUCoEbIhMpMjJSE+mXncaSrUVeRxEROUqhr5LEeCM1KbqWxWvs7KYvnHMnAqOAUc65McA5EU0WAWN6dmTptkKc0y0eIhJdikr9ZKYmRcVGQ7U1trsJAOfcgdCd1wB3RyBPRI3ulcnekkryC2Oqp0xE2oCisko6Rtl4BDSxSBwhuspdI4zpmQloXEJEok+hz09mSnSNR8DxFYmY67MZ0q09yYlxLNW4hIhEmSJfZdTNbIIGBq7N7CB1FwMDUiKSKIIS4uMY1SOTJdsKvY4iInKYIp+fjj1jrCXhnGvvnOtQx0d751xjZ0ZFldG9Mlm5/QAVVQGvo4iIAOCco8jnJzMt+loSx9PdFJPG9MykMlDN6p0HvY4iIgKArzJAZaC61Y1JxKTRvTIBWLI1fJfTwo37uO3lRQSqY27YRURiUFFZzY10akl4rntGCt06JLO0nhlOT83byLvLd7GhoKTlgolIm1VYWrMkh1oSUWF0z8ywd14fKPcz/8u9AFrnSURaRFFoSY5onN3UJovEmF6ZbN3vY19JxVGPfbhqN5WBauCr/WZFRCKpKLS4X7St2wRttEiMDt1UV1eX0zvLdtIjM4WTendUkRCRFvHV4n5qSUSFkTkZxMfZUV1OxWXBrqYLR3RjZI8MVu04oMFrEYm4otCYRDTtbV2jTRaJ1KQEBndtz8fr9x622F9NV9NFo7ozskcGZf6ABq9FJOIKfX5Sk+JplxBdK8BCGy0SANeO78XSbUU88/GmQ8feXR7sahrdM5NROcFV0TV4LSKRFlzcL/rGI6ANF4lvje/FecO68tv31rAsv4jiMj/zvizgwhHdMDP6ZaeTmhSvcQkRibginz8qZzZBGy4SZsbvvjGK7PR23DF9CW8u3Y4/4LhoVHcA4uOMYd07sEJFQkQirDBKF/eDNlwkIHjjykPXjiG/sIxfvr3qUFdTjRE9MlipwWsRibBinz8qb6SDNl4kAHL7dOLu8wYRqHaHuppqaPBaRFpCoS86NxyCxu9x3ardOqE/7ZMTuGB4t8OOj6w1eD2oa3svoolIK1dd7Sgui84Nh0AtCSA4/nD9qX3o0iH5sOP9s9NJSax/8Hrz3lK+0E53InKMDpZXUe2ic0kOUEuiXvFxxvAT6h68LvcHeGzuBv48dz1mxic/OYes9HbNct39pZX8/I3lbNhTSmq7eNLbJdApLYm7zxtE76y0Y37dyqpqHpmzng7JCQw7oQPDu2eQEG/MW1fA7FW7mbt2Dyf2zOTBq0bTMS38XzXl/gB/W5TPzqIyqqodlVXVpCbF861TenNCZsztRSXiqUJf9C7JASoSDRrRI4NXPt9GoNoRHxccr/hkwz5+9vpyNu4t5dyhXfhw9R6mf7aVH5wz8Livtzy/mFtfWkRBSQVnDcymoipAaUUVS7cWsXDjPl6+eTwDuhxb19evZ63ixYVbDjuWEGdUVTs6piZyav8sPly1h0sf/Zinrs9lSLcOh53rnGP2qt088M5qtu73kRBnJMbHkRBvlFUGeObjTXznjL58f2J/OiRH519FItHmUJGIwg2HQEWiQSN7ZPD8gs1sKCjhhMwUHnhnNdM/20rvrFRevGkcZw7M5rpnPuXFhVu45az+JCUcew/e3/K28bM3VtA5LYmZt57KqJzMQ4+t3XWQqU9/ytVPLOQvN41j+AkZTXrt1xbnhzL245az+rFyxwFW7ijmQFkVEwdnk9u7IwnxcSzZWsitLy3iiscW8JsrRjKgSzr7SyvZX1rJ3/Ly+Xj9XgZ1Teflm8dz+oDOh14/v9DH72ev489zNzDjs618f2J/rh3Xi/YqFiL1qtlLIiNKxySs9rIUzfrCZs8CFwN7nHMj6nh8KnAvwf2yDwLfd8590dDr5ubmury8vOaOG9a63Qc5/8F5TB3fi7lrC9hRXMZ3z+zH3ecNIjkxeAv9v9bs5jvP5/Gna0Zz2egeTb5GoNrx61mreH7BZk4fkMVD14yps+tq095Spj61kJKKKl74zjjG9OrYqNdfvfMAX3vs34zumclLN40nIb7+QrbnQDnfe2nRUWtbdUhO4O7zBvGtU3qHfY0V24v57XtrmP/lXtq3S+CacT258fS+JMYbW/f52LzPR1WgmnOGdDlqDEikLXp9ST53vfIF//rRBPplp0fsOma2yDmX2+TnRbBInAWUAH8JUyROA1Y75wrN7EJgmnNufEOv29JFIlDtGHHf+5T5A/TrnMbvrhzFSb07HXZOdbVj0h8+okNKIm/efvphx7cXldGzU2rY1y+pqOLO6Uv415o93HRGX3564ZB6f4nnF/qY+vSnFBys4InrTuLMgdn15i8u83PpIx9T7g8w644zyW7fuHGTiqoA76/cTWKc0Sktiaz0JHpkppKS1Li1ZZblF/HU/E28u3xnnfeZmMG4Pp24eFR3Jg/vpoIhbdazH2/iV7NWseS/zqt3LPB4RV2RADCzPsCsuorEEed1BFY45xr8M7yliwTAQ//8knJ/gDsnDTzUejjSCws2c99bK3ntttMY26sjZZUB/mPGEmav2s09kwdz+9kDjnrOzuIyvvN8Hut2H2TapcO57pTejcqz50A51z/7GRsKSvj9VaO59MQTDnu83B/g8837+WhtAR+s3s32wjJe+d4pRxW3lpBf6OPNpTtIS4qnd+c0endKxR9w/GPFTmYt28n6PcF7UE7smcl5Q7swaWhXhnRrf9j9KjWKy/wszy9m2fYilm0rZm9JBZeN6cHXx/YgNUk9pxKb/jB7LQ/PWc/6B6YcGveMhFgvEj8Ghjjnbg7z+C3ALQC9evU6acuWLXWd5qnSiipO+c0/mTAom2mXDuemF/JYll/EmJ6ZLN5axA/OHsCPzh+EmeGcY9aynfxq1irKKgM8OnUsEwbV3yI4UnGZn+++kMfnW/Yz7ZLhnDmwMx+tK+CjdQUs3LiPcn81SfFxjOvbietP7c35R9wDEi3W7T7IB6t2M3vV7kNTiTunt+PU/lmc1j+LjJREPtu0n0837WfNrgPUvF17Z6WSkhjPml0H6ZCcwDXjevHt0/vQPUOzqyS2/OLNFby5dAdf3Hd+RK8Ts0XCzM4GHgPOcM7ta+g1vWhJNNb9s1bx3ILNdM9IZm9JBX+6ZgznDu3Kz99YzvTPtvGd0/vyjZNymPb2Sj7btJ/hJ3Tg91edeNQsosYq9we4c3qwtVKjX+c0zhqUzYRB2Yzv1ymm/sLec6CcuWsLWLBhLws27GPPweDOgcmJcZzUuyMn9+nE2F4dGZWTQWZqEs45Fm8t4tl/b+K9FbtITYzn/q+NOKZxIRGv3DF9Ccvzi5h7z9kRvU5MFgkzGwW8DlzonFvXmNeM5iKxbb+PCb+bQ6e0JJ6+4eRD60A55/jVrFU89+/NQHD3qR9PHsw1J/c67uZlVaCalxZuIT4+jgkDs+mVFX78I5Y459hQUEpJRRXDundocNbY1n0+7np1KYu2FPL1sTn88rLhpLcLFsgiXyWGReWGLiLXPfMpB8ureKPWeGYkHGuR8OzPTDPrBbwGXNfYAhHtenZK5ZXvnUpOx5TDuj3MjF9cPIyuHZLZe7CCH5wzoNkW80qIj+PG0/s2y2tFEzNjQJfGz/TolZXKK7ecwkP/Ws8j//qSTzfto3N6OzbvK6XI5ycpIY7ffn0kXxuTE8HUIk1X5POTlR6d018hgkXCzKYDE4HOZpYP3AckAjjnHgd+AWQBj4UGKauOpcpFm5P71D04bGbcOqF/C6dpWxLi47j7vEGc3j+L3763hpSkeKaM7E7frDQ+XL2bu175gtU7D3LvBUMiOkAo0hSFvsom/UHU0iJWJJxz1zbw+M1AnQPVIsdjfL8sXrvt8Kb7jaf34f5Zq3hy3kZW7zzAn64ZQ6cITjcUaaxo3nAItMCftBGJ8XH88rIR/M8VI1m4cR+n/uaf/HDGEhZs2Eu19gsRj/gD1ZRUVEXtCrCgZTmkjblmXC/G9u7ISwu38PqS7byxdAe9s1K5+Yy+XJnbM+x9MCKRUOQLLskRres2gVoS0gYN6tqeX102gs9/di4PXn0indKS+K83V3Lm/87hiY82UFJR5XVEaSOKQov7ReuudKCWhLRhyYnxfG1MDpeP7sEnG/fx2JwN/OYfa3hkznquHdeLG07rQw8tfS4RVLO4X7TuSgcqEiKYGaf178xp/TuzdFsRT83fyDMfb+KZjzdxwYhu/Oi8QRFdeE3arsLSUEtCYxIisWF0z0we/eZYtheV8ZcFm/nrp1uZu2YPv/n6qKPWyBI5XjVjEprdJBJjemSm8NMpQ3n/rrMY0r0Dd05fws/fWE65P4Bzjn0lFeRt3s+iLfvZH/prUKSpvtpwSC0JkZh0QmYKM245hf97fy1PzNvIB6t2U1YZ4ED54YPbGSmJDOiSzj2TB3NKvyyP0kqsKSrzkxBnpDVyCX4vqEiINCAxPo6fThnKuL6dmP7ZNrpnJNO3cxp9O6fhcGwsKGXT3lLmf7mX6575lP/+2kiuzO3pdWyJAUW+SjJTk+pcGj9aqEiINNKkoV2ZNLTrUcfPGRL8t7jMz+0vL+aemcvYtLeUH58/mDgt/yH1KCz1R/XMJlCREGk2GSmJPPftk/nFmyt5bO4Glm4r4rT+WQzoks6ALun065yuoiGHKSqrjOpBa1CREGlWifFx/PfXRjCwSzrPfLyJBRu+2iJldM9MHrpmTKtZzl2OX5HPX+/2xtFARUKkmZkZ3zmjL985oy+lFVVsKChhydYi/m/2WqY8NJ/7Lx/B5WO0MVJb55xj235f1E90UJEQiaC0dgmMyslkVE4mk4Z24YczlvLDV5bywardjMrJICkhjqSEOIZ0a+/JHuTinT0HKyitDNAvO83rKPVSkRBpITkdU5lxyyk8OmcDj85ZzzvLdx72+MWjuvPzi4bRLSPZo4TSkjYUlADQr3N0382vIiHSghLi4/iPcwdy+9n9qaiqprKqmvKqAK9+ns+jc9czZ80e7jpvEDee1oeEeN3r2pptLCgFiPqWhN6FIh5IiI8jrV0CHdOS6J6Rwn+cO5AP75rA+H5Z3P/Oaq59aiG7D5R7HVMiaGNBKcmJcXTrEN0tRxUJkSjRKyuVZ27I5U/XjGbljgNc9NB8FmzY63UsiZCNe0voGwPTolUkRKKImXHZ6B68efvpZKQk8q2nP+VPH35JcWhJaWk9Nu0tjfquJlCREIlKA7u2580fnMGUkd158MN1nHz/h9zylzzeWbaTcn/A63hynCqqAmzb76N/5+gvEhq4FolS6e0SePjaMXz3zH68uXQHby/bwexVu+mRmcJ9lwzj/MA8+OevoDgfMnJg0i9g1FVex5ZG2LrPR7UjJvYpUZEQiWJmxok9MzmxZyY/u2go878s4DfvrmHWy39iQrtnaOcqgicWb4O37wx+rkIR9TaEZjb1jYGWhLqbRGJEfJwxcXAXZt15Bg+0f/2rAlHDXxZsWUjU27g3dI9EDIxJqCUhEmMS4+NIrNhV94PF+S0bRo7JpoJSstu3o31ydC/uB2pJiMSmjJw6D1e115pQsWDj3lL6xUBXE6hIiMSmSb+AxJTDDpW5JO4pvIzz/vAR339pEX+YvZa1uw56FFDqs7GgJCa6mkDdTSKxqWZwutbspuKT/5NuJeMo2VPCml0HeX/lLh6bu4FbJ/TnB+cMIDkxerfIbEsKSysp9Pmjfs2mGioSIrFq1FWHzWTqBtxb6+HC0kruf2c1j8xZz7vLd/LfV4yM+mWp24JYGrQGdTeJtFod05L4/VUn8uJN46gMVHPNkwu55S95rNpxwOtobdpXC/vFRktCRUKklTtzYDaz7zqLu84dxCcb9zHlofnc9vIi1u3WeIUXNu4tJSHO6NkxpeGTo4CKhEgbkJqUwH+cO5CP//Mc7jhnAB+tLWDyH+dx5/QlrN9T4nW8NmVjQQm9slJjZin42EgpIs0iIzWRH50/mPn3nsOtE/rz4erdnP/gR9z1ytJDm+BIZG0sKI2ZQWvQwLVIm9QpLYl7LxjCzWf05cl5G3nhk828sXQ7U0Z057az+zP8hIyjnuOc49/r9/FK3jacc3RISaR9cgI9O6YyZWR3OqUlefCdxJZAtWPLPh/nDOnidZRGU5EQacOy0tvx0ylD+e5Z/Xj24028+MkW3lm+kzMHdubU/lmMzslkRE4GizYX8tC/vmTJ1iKy0pLokJLIwXI/B8qrqKyq5pdvr+ScIV34xkk9mTg4m8QY6UppadsLy6gMVMfMzCZQkRARoHN6O/7zgiF8b0J//rJgM68t2c7/vrf2sHN6ZKZw/+UjuDI3h3YJX91zsWbXAWbm5fPG0u28v3I3mamJXDiiGxePOoFT+mURH+Wb6rSkDaHpr31jqLvJnHNeZ2iS3Nxcl5eX53UMkVavyFfJsvxiluUX0S0jhctGn1BvC8EfqOajtQW8vWwHH6zaja8yQFZaEqcN6Mwp/TpxSr8s+nVOw6ztFo2n52/k/ndWk/fzc+mc3q5Fr21mi5xzuU19nloSIlKnzNQkzhqUzVmDsht1fmJ8HOcO68q5w7pSVhlg7to9vL9yF59s3MfbX+wAIKdjCleM6cEVY3PoEyNrFzWnJduKOCEjucULxPFQkRCRZpeSFM+FI7tz4cjuOOfYvM/HJxv28Y8VO3l4znoe+td6TurdkUtGdWfyiG50z4iNewaOh3OOvM37Gdc3tu56j1iRMLNngYuBPc65EXU8bsCfgCmAD7jRObc4UnlExBtmRt/OafTtnMY3x/diZ3EZbyzZwetL8pn29iqmvb2K0T0zuXhUd6aO701KUutcYyq/sIzdByrI7d3R6yhNEskpCM8DF9Tz+IXAwNDHLcCfI5hFRKJE94wUvj+xP7PvmsCHd0/gnsmDCVQ77n9nNWf/31xmLsqnujq2xkobY9GWQgBy+6hIAOCcmwfsr+eUy4C/uKCFQKaZdY9UHhGJPgO6pHP72QN4+44zePV7p9K1Qzt+/LcvuPjhj3l/5S78gWqvIzabvC37SW+XwJBuHbyO0iRejkn0ALbV+jo/dGynN3FExEvj+nbi9dtO5+1lO/jf99byvRcXkZWWxKWjT+AbJ+XUeYNfLMnbXMiYXpkxNyU4JgauzewWgl1S9OrVy+M0IhIpcXHGZaN7MGVkd+atK+Dvi/N5eeFWnvv3ZiYN6cJPpwxlQJfYucegRnGZn7W7D3LhiNjrLPGySGwHetb6Oid07CjOuSeBJyF4n0Tko4mIlxLj45g0tCuThnalyFfJXz/bymNzNjD5j/P41vhe/PDcQXSMoWVAlmwtxLnYG48Ab4vEW8APzGwGMB4ods6pq0lEDpOZmsRtEwdwVW5PHvxgHS8u3MKMz7dx9uAuXDSqO5OGdiE1Kbo7RfI2FxIfZ4zumel1lCaL5BTY6cBEoLOZ5QP3AYkAzrnHgXcJTn9dT3AK7LcjlUVEYl/n9HY88LWR3HBaH15euIV3V+zivZW7SE6M44bT+nD3eYMOWy4kmuRt2c+w7h1IaxfdxawuWpZDRGJSoNrx+eb9vPr5Nl5bsp0h3drzx2tGR93sIX+gmpHT3ueak3sx7dLhnuU41mU5tFSjiMSk+DjjlH5Z/OHq0Tx7Yy57Syq49OF/89S8jVRF0dTZlTsOUO6vjsnxCFCREJFW4JwhXXn/h2cxcXA2D7y7mgv/NJ9/rdlNNPSU5G0O3i6W27uTx0mOjYqEiLQKWenteOK6k3j8W2PxB6r5zvN5fPOpT1mxvdjTXIu2FJLTMYVuGcme5jhWKhIi0mqYGReM6M7suyYw7ZJhrNl1gEsf+Zhpb63kYLm/xfMEqh15Wwpjbr2m2lQkRKTVSUqI48bT+zL3nrOZOr43L3yymUm//4hZy3a0WBfU/C8LuOih+RQcrGDC4MYttx6NVCREpNXKSEnk15eP4I3bTqdLh3b84K9LeOCd1REtFOv3lHDDs59x3TOfUVpZxaPfHMvlo3tE7HqRFnuTdkVEmujEnpm8efsZ/PLtlTz98SY6pQdv0Gtu63Yf5KonPqG62vHzi4Zy3am9o/bejcZSkRCRNiE+zph2yXCKy/z873tr6ZSaxDXjmm8tuK37fHzr6U9Jio/jb7efSu+s1rHznoqEiLQZcXHG775xIkU+P//v9eWkJMUzcXAX2rdLIC7OOFjuZ+m2IhZtKWRDQSlXjO3B2YO7NPi6u4rL+ebTC/EHqnnle62nQIDuuBaRNshXWcXUpz9lydYiAOIM2icncqDcj3NgFhzPKPL5mTg4m59fNOyo1WfL/QHW7ylh3e6DPDpnPbsPVPDX745nVE5my39DjXCsd1yrSIhIm1RSUcUHq3axr6SS4jI/RT4/WelJnNS7I6N7ZtIuIZ4XFmzmoX9+SZk/wIRB2VRUVVNUVklhqZ+dxWXUbKDXPjmBp67P5ZR+0bt/tYqEiEgE7C2p4MEP1vHppv10SE4gMzWJzJREenZKZVDX9gzulk7vrDQS46N7suixFgmNSYiI1KNm9dm2KrpLn4iIeEpFQkREwlKREBGRsFQkREQkLBUJEREJS0VCRETCUpEQEZGwVCRERCSsmLvj2swKgC1ABlCzL2FDn9f82xnYewyXrf2aTXn8yOP1fX1k1trHjiV3S2au/bkXP+uGjkVj5rqOt9b3R3Nlri9XUzKHyxnu89bynu7tnGv67kfOuZj8AJ5s7Oe1/s073ms15fEjj9f39ZFZjzd3S2b2+mfd0LFozNyW3h/Nlbmh3I3NHAvvD68z1/6I5e6mt5vwee1jx3utpjx+5PH6vq4r6/HkbsnMtT/34mfd0LFozFzX8db6/miuzA09v7GZjzwWje8PrzMfEnPdTcfDzPLcMSxw5bVYzK3MLScWcytzy2iOzLHckjgWT3od4BjFYm5lbjmxmFuZW8ZxZ25TLQkREWmattaSEBGRJlCREBGRsFQkREQkLBWJEDM708weN7OnzWyB13kaw8zizOwBM3vYzG7wOk9jmdlEM5sf+nlP9DpPY5lZmpnlmdnFXmdpDDMbGvoZzzSz73udp7HM7HIze8rMXjGz873O0xhm1s/MnjGzmV5nqU/oPfxC6Oc7tTHPaRVFwsyeNbM9ZrbiiOMXmNlaM1tvZj+p7zWcc/Odc7cCs4AXIpk3lO24MwOXATmAH8iPVNbamim3A0qAZFogdzNlBrgXeDUyKQ/XTO/p1aH39FXA6ZHMWytfc+R+wzn3XeBW4OpI5g1la47MG51zN0U2ad2amP8KYGbo53tpoy5wvHfjRcMHcBYwFlhR61g8sAHoByQBXwDDgJEEC0Htjy61nvcq0D4WMgM/Ab4Xeu7MWPlZA3Gh53UFXo6RzOcB1wA3AhfHQubQcy4F/gF8M1beH7We93tgbIxlbpH/D48j/0+B0aFz/tqY10+gFXDOzTOzPkccHgesd85tBDCzGcBlzrnfAHV2F5hZL6DYOXcwknmheTKbWT5QGfoyEMG4hzTXzzqkEGgXkaC1NNPPeiKQRvB/tDIze9c5Vx3NmUOv8xbwlpm9A/w1UnlrXa85ftYG/A/wD+fc4ghHbu73dItrSn6CLfccYCmN7ElqFUUijB7Atlpf5wPjG3jOTcBzEUvUsKZmfg142MzOBOZFMlgDmpTbzK4AJgOZwCMRTRZekzI7534GYGY3AnsjWSDq0dSf80SC3QvtgHcjGawBTX1f3wGcC2SY2QDn3OORDBdGU3/WWcADwBgz+2momHgpXP6HgEfM7CIauXRHay4STeacu8/rDE3hnPMRLGwxxTn3GsECF3Occ897naGxnHNzgbkex2gy59xDBH+ZxQzn3D6CYyhRzTlXCny7Kc9pFQPXYWwHetb6Oid0LJrFYmaIzdzK3HJiMXcsZq6t2fK35iLxOTDQzPqaWRLBQce3PM7UkFjMDLGZW5lbTizmjsXMtTVf/pYeiY/Q6P50YCdfTQW9KXR8CrCO4Cj/z7zOGeuZYzW3Mit3a8vckvm1wJ+IiITVmrubRETkOKlIiIhIWCoSIiISloqEiIiEpSIhIiJhqUiIiEhYKhLSKphZSQtfr1n2HLHg3hrFZrbUzNaY2f814jmXm9mw5ri+SENUJETqYGb1rmvmnDutGS833zk3GhgDXGxmDe39cDnB1WhFIk5FQlotM+tvZu+Z2SIL7oQ3JHT8EjP71MyWmNmHZtY1dHyamb1oZv8GXgx9/ayZzTWzjWZ2Z63XLgn9OzH0+MxQS+Dl0FLXmNmU0LFFZvaQmc2qL69zrozgEs49Qs//rpl9bmZfmNnfzSzVzE4juEfE70Ktj/7hvk+R5qAiIa3Zk8AdzrmTgB8Dj4WOfwyc4pwbA8wA/rPWc4YB5zrnrg19PYTgsubjgPvMLLGO64wBfhh6bj/gdDNLBp4ALgxdP7uhsGbWERjIV8u+v+acO9k5dyKwmuByCwsIrsFzj3NutHNuQz3fp8hx01Lh0iqZWTpwGvC30B/28NUGRznAK2bWneCuXZtqPfWt0F/0Nd5xzlUAFWa2h+BuekduufqZcy4/dN2lQB+C27NudM7VvPZ04JYwcc80sy8IFog/Oud2hY6PMLP7Ce67kQ6838TvU+S4qUhIaxUHFIX6+o/0MPAH59xboY15ptV6rPSIcytqfR6g7v9nGnNOfeY75y42s77AQjN71Tm3FHgeuNw590Vos6OJdTy3vu9T5Lipu0laJefcAWCTmV0JwS0xzezE0MMZfLW2/g0RirAW6FdrW8mrG3pCqNXxP8C9oUPtgZ2hLq6ptU49GHqsoe9T5LipSEhrkWpm+bU+7ib4i/WmUFfOSoJ7/EKw5fA3M1sE7I1EmFCX1W3Ae6HrHASKG/HUx4GzQsXlv4BPgX8Da2qdMwO4JzTw3p/w36fIcdNS4SIRYmbpzrmS0GynR4EvnXMPep1LpCnUkhCJnO+GBrJXEuziesLbOCJNp5aEiIiEpZaEiIiEpSIhIiJhqUiIiEhYKhIiIhKWioSIiISlIiEiImH9fxcU/v1EIJnHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "learn.lr_find()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ea05e63-e468-4fe8-9679-f1f10eba36ea",
   "metadata": {},
   "source": [
    "And finally we can train the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed45964a-7729-4023-bd58-06474c817f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "learn.fine_tune(3, 3e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be6db71d-7527-4dd6-bf7c-a1f6e65e7c8e",
   "metadata": {},
   "source": [
    "Note that this time we didn't call `fit()` or `fit_one_cycle()` as we've done in previous lectures. That's because we're starting with a pretrained model whose weights we \"fine-tuned\" to the new dataset of jet images. This process involves some tricks which is why it is called `fine_tune()`. Anyway, we've now got a model that scores pretty well on the validation set. Let's compute the same metrics we did in lecture 1, but now on the test set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f7737f-487b-43d0-ac6c-ed2df45e8587",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_metrics(learn, test_ds):\n",
    "    test_dl = learn.dls.test_dl(test_items=test_ds)\n",
    "    preds, targs = learn.get_preds(dl=test_dl)\n",
    "    fpr, tpr, _ = roc_curve(y_true=targs, y_score=preds[:, 1])\n",
    "    acc_test = accuracy_score(targs, preds.argmax(dim=-1))\n",
    "    auc_test = auc(fpr, tpr)\n",
    "    background_eff = interp1d(tpr, fpr)\n",
    "    background_eff_at_30 = background_eff(0.3)\n",
    "\n",
    "    print(f\"Accuracy: {acc_test:.4f}\")\n",
    "    print(f\"AUC: {auc_test:.4f}\")\n",
    "    print(\n",
    "        f\"Backround rejection at signal efficiency 0.3: {1/background_eff_at_30:0.3f}\"\n",
    "    )\n",
    "    return fpr, tpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bf680b5-6fcd-4c91-b27f-f9a8687a1891",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_ds = get_dataset(images_ds[\"test\"])\n",
    "_, _ = compute_metrics(learn, test_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7682d011-6728-4be9-a6e8-ccaad8c5b6b7",
   "metadata": {},
   "source": [
    "Comparing our results against those in the top tagging review, shows that our CNN is getting quite competitive results with the state-of-the-art!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b866eb4a-a3a1-4dd1-909d-33231b351cc3",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img alt=\"transfer\" caption=\"Transfer learning\" src=\"images/top-tagging-scores.png\" id=\"transfer-learning\" width=800/>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f3cb02a-45b4-451a-8732-4d43425e987b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}