<jupyter_start><jupyter_text>Particle Net ReplicationOriginal Authors:* <NAME>* <NAME>Original Code:https://github.com/hqucms/ParticleNetSource Research Paper:https://arxiv.org/pdf/1902.08570.pdf Authors of this Notebook:* <NAME>* <NAME>* <NAME>  Download Data and Packages  Installing LibrariesPython libraries being installed:* Awkward: it's used with numpy to manipulate operations on numpy arrays without needing to come back to regular python.* Uproot-methods: it's used as helper methods to uproot, used to manipulate root files.<jupyter_code>!pip install awkward>=0.12.0
!pip install uproot-methods>=0.7.0
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt
import tensorflow as tf
import awkward
import uproot_methods
import logging
import urllib.request
import os

logging.basicConfig(level=logging.DEBUG, format='[%(asctime)s] %(levelname)s: %(message)s')
logging.getLogger('matplotlib.font_manager').disabled = True<jupyter_output><empty_output><jupyter_text>Getting the Data<jupyter_code>from google.colab import drive
drive.mount('/content/gdrive/')
PROJECT_DIR = '/content/gdrive/My Drive/Scratch Space/particle-net'
# urllib.request.urlretrieve(
#     'https://zenodo.org/record/2603256/files/train.h5?download=1',
#     os.path.join(PROJECT_DIR, 'original/train.h5')
# )
# print('Downloaded train.h5')

# urllib.request.urlretrieve(
#     'https://zenodo.org/record/2603256/files/val.h5?download=1',
#     os.path.join(PROJECT_DIR, 'original/val.h5')
# )
# print('Downloaded val.h5')

# urllib.request.urlretrieve(
#     'https://zenodo.org/record/2603256/files/test.h5?download=1',
#     os.path.join(PROJECT_DIR, 'original/test.h5')
# )
# print('Downloaded test.h5')

os.listdir(os.path.join(PROJECT_DIR, 'original'))<jupyter_output><empty_output><jupyter_text>What the Dataset looks like[Top Quark Tagging Reference Dataset](https://zenodo.org/record/2603256.XsaZZWgzZPY)A set of MC simulated training/testing events for the evaluation of top quark tagging architectures.In total 1.2M training events, 400k validation events and 400k test events. Use “train” for training, “val” for validation during the training and “test” for final testing and reporting results.* $14 TeV$, hadronic tops for signal, qcd diets background, Delphes ATLAS detector card with Pythia8* No MPI/pile-up included* Clustering of  particle-flow entries (produced by Delphes E-flow) into anti-kT 0.8 jets in the $pT$ range $[550,650] GeV$* All top jets are matched to a parton-level top within $\Delta R = 0.8$, and to all top decay partons within $0.8$* Jets are required to have $|\eta| < 2$* The leading 200 jet constituent four-momenta are stored, with zero-padding for jets with fewer than 200* Constituents are sorted by $pT$, with the highest $pT$ one first* The truth top four-momentum is stored as truth_px etc.* A flag (1 for top, 0 for QCD) is kept for each jet. It is called is_signal_new* The variable "ttv" (= test/train/validation) is kept for each jet. It indicates to which dataset the jet belongs. It is redundant as the different sets are already distributed as different files.  Understanding the Data  The Tables<jupyter_code>df = pd.read_hdf(os.path.join(PROJECT_DIR, 'original/val.h5'), key='table')
df.head()<jupyter_output><empty_output><jupyter_text>So each row is an event. For each row:There are Columnns for all the particles 4-coordinate values in Energy-Momenta terms labelled `PX_i`, `PY_i`, `PZ_i`, and `E_i`.There is a 4-coordinate value called Truth, with labels `TruthPX`, `TruthPY`, `TruthPZ`, `TruthE`, I presume that it's the true coordinate value of the jet.`is_signal_new` is the label associated with the jet, whether it's a Top or a QCD jet.There is a train/test/validation identifier called `ttv`.  Number of Particles<jupyter_code>def _col_list(prefix, max_particles=200):
    return ['%s_%d'%(prefix, i) for i in range(max_particles)]
plt.figure(figsize=(25, 5))
g = sns.countplot(np.count_nonzero(df[_col_list('E')].values, axis=1))
g.set_xticklabels(g.get_xticklabels(), rotation=90)
plt.title('Number of Particles in a Jet')
plt.show()<jupyter_output><empty_output><jupyter_text>Preparing the DataThis time doing it the way it's done in the paper. $px$, $py$, $pz$ and $E$ are all in cartesian coordinates.We convert it to a Cylinderical coordinates, z-axis along the directions of the collision. The transverse axis should be perpendicular to this collision axis, we convert to $\theta$ and $\phi$ coordinates.Root Documentation: https://root.cern.ch/doc/master/classTLorentzVector.html<jupyter_code>def _transform(df, start=0, stop=-1):
    r"""
    Takes a DataFrame and converts it into a Awkward array representation
    with features relevant to our model.

    Parameters
    ----------
    df: Pandas DataFrame
        The DataFrame with all the momenta-energy coordinates for all the particles
    start: int
        First element of the DataFrame
    stop: int
        Last element of the DataFrame
    
    Returns
    -------
    v: OrderedDict
        A Ordered Dictionary with all properties of interest

    Notes
    -----
    Here the function is just computing 4 quantities of interest:
    * Eta value relative to the jet
    * Phi value relative to the jet
    * Transverse Momentum of the Particle (log of it) 
    * Energy of the Particle (log of it) 
    """
    from collections import OrderedDict
    v = OrderedDict()

    df = df.iloc[start:stop]
    # We take the values in the dataframe for all particles of a single event in each row
    # px, py, pz, e are in separate arrays
    _px = df[_col_list('PX')].values
    _py = df[_col_list('PY')].values
    _pz = df[_col_list('PZ')].values
    _e = df[_col_list('E')].values

    # We filter out the non-0 non-negative energy particles
    mask = _e > 0
    n_particles = np.sum(mask, axis=1) # Number of particles for each event where energy is greater than 0

    # _p[mask] filters out the >0 energy particles, and flattens them, so that they can be recollected for each event from counts array.
    px = awkward.JaggedArray.fromcounts(n_particles, _px[mask])
    py = awkward.JaggedArray.fromcounts(n_particles, _py[mask])
    pz = awkward.JaggedArray.fromcounts(n_particles, _pz[mask])
    energy = awkward.JaggedArray.fromcounts(n_particles, _e[mask])
    # These are jagged arrays with each row for 1 event, and all particles in the row

    p4 = uproot_methods.TLorentzVectorArray.from_cartesian(px, py, pz, energy)
    jet_p4 = p4.sum() # Sum of Lorentz Vectors for all particles in the jet

    # Getting the Labels
    _label = df['is_signal_new'].values # the target labels, QCD or Top
    v['label'] = np.stack((_label, 1-_label), axis=-1) # Making it categorical [Top, QCD]
    # Transformed features relative to the Jet and log features
    v['part_pt_log'] = np.log(p4.pt)
    v['part_e_log'] = np.log(energy)
    # Flip particle ETA if Jet Eta is negative
    # Particle's phi relative to the Jet
    _jet_etasign = np.sign(jet_p4.eta)
    _jet_etasign[_jet_etasign==0] = 1
    v['part_etarel'] = (p4.eta - jet_p4.eta) * _jet_etasign
    v['part_phirel'] = p4.delta_phi(jet_p4)

    del px, py, pz, energy, _px, _py, _pz, _e, p4, jet_p4, df
    return v
def convert(source, destdir, basename, step=None, limit=None):
    r"""
    Converts the DataFrame into an Awkward array and performs the read-write
    operations for the same. Also performs Batching of the file into smaller
    Awkward files.

    Parameters
    ----------
    source: str
        The location to the H5 file with the dataframe
    destdir: str
        The location we need to write to
    basename: str
        Prefix for all the output file names
    step: int
        Number of rows per awkward file, None for all rows in 1 file
    limit: int
        Number of rows to read.
    """

    if limit is not None:
        df = pd.read_hdf(source, key='table', start=0, stop=limit)
        logging.info('Restricting to the first %s events:' % str(df.shape[0]))
    else:
        df = pd.read_hdf(source, key='table')

    if step is None:
        step = df.shape[0]
    logging.info('Total events: %s' % str(df.shape[0]))

    idx = 0
    # Generate files as batches based on step size, only 1 batch is default
    for start in range(0, df.shape[0], step):
        if not os.path.exists(destdir):
            os.makedirs(destdir)
        output = os.path.join(destdir, '%s_%d.awkd'%(basename, idx))
        logging.info(output)

        if os.path.exists(output):
            logging.warning('... file already exist: continue ...')
            continue
        v = _transform(df, start=start, stop=start+step)
        awkward.save(output, v, mode='x')
        idx += 1

    del df, output
# convert(os.path.join(os.path.join(PROJECT_DIR, 'original'), 'train.h5'), 
#         destdir=os.path.join(PROJECT_DIR, 'converted'),
#         basename='train_file',
#         limit=500000)
# convert(os.path.join(os.path.join(PROJECT_DIR, 'original'), 'val.h5'), 
#         destdir=os.path.join(PROJECT_DIR, 'converted'), 
#         basename='val_file')
# convert(os.path.join(os.path.join(PROJECT_DIR, 'original'), 'test.h5'), 
#         destdir=os.path.join(PROJECT_DIR, 'converted'), 
#         basename='test_file')<jupyter_output><empty_output><jupyter_text>Making the DataSet<jupyter_code>def pad_array(jagged_array, max_len, value=0., dtype='float32'):
    r"""
    Pads a Jagged array with value to make them equal to max_len along each column

    Parameters
    ----------
    jagged_array: awkward array
        The array to be padded
    max_len: int
        Number of columns in each row (truncate or extend to)
    value: int or float
        The value to pad with
    dtype: string or type
        Final Data type of array

    Returns
    -------
    rectangular_array: np.array
        Padded version of the input Jagged Array
    """
    rectangluar_array = np.full(shape=(len(jagged_array), max_len), fill_value=value, dtype=dtype)
    for idx, jagged_element in enumerate(jagged_array):
        if len(jagged_element) != 0:
            trunc = jagged_element[:max_len].astype(dtype)
            rectangluar_array[idx, :len(trunc)] = trunc
    return rectangluar_array
class Dataset(object):

    def __init__(self, filepath, value_cols = None, label_cols='label', pad_len=100, data_format='channel_first'):
        self.filepath = filepath
        self.value_cols = value_cols if value_cols is not None else {
            'points': ['part_etarel', 'part_phirel'],
            'features': ['part_pt_log', 'part_e_log', 'part_etarel', 'part_phirel'],
            'mask': ['part_pt_log']
        }
        self.label_cols = label_cols
        self.pad_len = pad_len
        assert data_format in ('channel_first', 'channel_last')
        self.stack_axis = 1 if data_format=='channel_first' else -1
        # Here we make the arrays which will keep out data and load the first batch
        self._values = {}
        self._label = None
        self._load()

    def _load(self):
        logging.info('Start loading file %s' % self.filepath)
        counts = None
        with awkward.load(self.filepath) as a:
            # Load output labels from the awkward array
            self._label = a[self.label_cols]

            for k in self.value_cols:
                cols = self.value_cols[k]
                assert isinstance(cols, (list, tuple))
                arrs = []
                for col in cols:
                    if counts is None:
                        counts = a[col].counts
                    else:
                        assert np.array_equal(counts, a[col].counts)
                    arrs.append(pad_array(a[col], self.pad_len))
                self._values[k] = np.stack(arrs, axis=self.stack_axis)
        logging.info('Finished loading file %s' % self.filepath)


    def __len__(self):
        return len(self._label)

    def __getitem__(self, key):
        return self._label if key == self.label_cols else self._values[key]
    
    @property
    def X(self):
        return self._values
    
    @property
    def y(self):
        return self._label

    def shuffle(self, seed=None):
        # Get a random permutation
        if seed is not None: np.random.seed(seed)
        shuffle_indices = np.random.permutation(self.__len__())
        # Reorder the table
        for k in self._values:
            self._values[k] = self._values[k][shuffle_indices]
        self._label = self._label[shuffle_indices]<jupyter_output><empty_output><jupyter_text>The dataset exposes `_values` and `_labels`* `_values` is a dict with the following items, they will be of dimentions (Examples, MaxLen, Number of Columns),:    * `features`: 4 columns    * `points`: 2 columns    * `masks`: 1 column* `_labels` is the categorical output targets    * Top Jet Column: one for each jet    * QCD Jet Column: one for each jet<jupyter_code>train_dataset = Dataset(os.path.join(PROJECT_DIR, 'converted/train_file_0.awkd'), data_format='channel_last')
val_dataset = Dataset(os.path.join(PROJECT_DIR, 'converted/val_file_0.awkd'), data_format='channel_last')<jupyter_output><empty_output><jupyter_text>Visualizing the Dataset<jupyter_code>fig, ax = plt.subplots(3, 3, figsize=(15, 15))
ax = np.reshape(ax, (-1))
for i in range(len(ax)):
    event_id = np.random.choice(np.reshape(np.argwhere(val_dataset['label'][:, 0] == 1), (-1)), len(ax))
    ax[i].scatter(
        val_dataset['points'][event_id][:, 0], 
        val_dataset['points'][event_id][:, 1], 
        color='k', s=3)
fig.suptitle('Top Jets')
plt.show()
fig, ax = plt.subplots(3, 3, figsize=(15, 15))
ax = np.reshape(ax, (-1))
for i in range(len(ax)):
    event_id = np.random.choice(np.reshape(np.argwhere(val_dataset['label'][:, 1] == 1), (-1)), len(ax))
    ax[i].scatter(
        val_dataset['points'][event_id][:, 0], 
        val_dataset['points'][event_id][:, 1], 
        color='k', s=3)
fig.suptitle('QCD Jets')
plt.show()<jupyter_output><empty_output><jupyter_text>Making the Model  Helper Functions for Nearest Neighbors<jupyter_code>def batch_distance_matrix_general(A, B):
    """
    Gives the Euclidian distance between all possible pairs of particles of a
    single Jet
    Distance is computed between the channel vectors
    All the separate Jets are batch processed

    Parameters
    ----------
    A, B: TensorFlow Tensors of shape (N, P_A, C and N, P_B, C)
    
    Returns
    -------
    D: TensorFlow Tensors of shape (N, P_A, P_B)
        The distance between the two 

    Notes
    -----
    The Distance Formula used here is:
        D = A ** 2 + B ** 2 - 2 * A * B
        (x1 - x2) ** 2 + (y1 - y2) ** 2
    """
    with tf.name_scope('dmat'):
        # For all N rows: [(P_A, C) -> (P_A)], get A**2
        r_A = tf.reduce_sum(A * A, axis=2, keepdims=True)
        # For all N rows: [(P_B, C) -> (P_B)], get B**2
        r_B = tf.reduce_sum(B * B, axis=2, keepdims=True)
        # For all N rows: [(P_A, C) * (C, P_B) -> (P_A, P_B)]
        m = tf.matmul(A, tf.transpose(B, perm=(0, 2, 1)))
        D = r_A - 2 * m + tf.transpose(r_B, perm=(0, 2, 1))
        return D
def knn(num_points, k, topk_indices, features):
    """
    To select the features for the K-nearest neighbors provided

    Prarameters
    -----------
    num_points: int
        number of particles per jet
    K: int
        the top K neighbors we want
    topk_indices: Tensorflow Tensor of shape (N, P, K)
        the list of those K neighbors
    features: Tensorflow Tensor of shape (N, P, C)
        the features layer from the previous layer

    Returns
    -------
    features: TensorFlow Tensor of shape (N, P, K, C)
        C length feature vector of the closest K for all P particles in N jets
    """
    with tf.name_scope('knn'):
        queries_shape = tf.shape(features) # Stores the literal shape tuple (N, P, C)
        batch_size = queries_shape[0] # Stores the literal values N

        # N array stacked of size (P, K, 1) with all values equal to the Batch ID, overall shape (N, P, K, 1)
        batch_indices = tf.tile(tf.reshape(tf.range(batch_size), (-1, 1, 1, 1)), (1, num_points, k, 1))
        # Top Indices get expanded to the shape (N, P, K, 1)
        # The concatenated results is (N, P, K) pairs with (jet_index (batch_index), closest_neighbor), net shape (N, P, K, 2)
        indices = tf.concat([batch_indices, tf.expand_dims(topk_indices, axis=3)], axis=3)
        # We chose all the C features of the top K neigbors for each particle, done for the entire batch of N
        # Output shape = (N, P, K, C)
        return tf.gather_nd(features, indices)<jupyter_output><empty_output><jupyter_text>Defining the Layers ![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABoMAAAM4CAYAAAAHzZsOAAAgAElEQVR4AeydB1RU59b33/W97zWJGkDEhthLNE1jTMxNT0w3zVjTk5t6jenRWLHX2HtBVBDpvfc+hYHplGFghq6JsXcs/2/tZxhEIwqK1D1rzRplzpzy2/s855z9f/be/3Px4kXwmxmwD7APsA+wD7APsA+wD7APsA+wD7APsA+wD7APsA+wD7APsA+wD7APsA+wD7APtE4f+J/Lly+D38yAfYB9gH2AfYB9gH2AfYB9gH2AfYB9gH2AfYB9gH2AfYB9gH2AfYB9gH2AfYB9oHX6AItBLI<KEY>ZBlYcMr+IMqghMSN+n9HIEhiIjI7F+/XqRSq5Pnz7Yv38/Ll68WO+uPul5KsMrgThdXV2YmJiwO6jGZhveESvACrACrIDKKSAFCQQTZJ5UL3Hy5Mm4dOlSifdFcInjSwp7mUgD0rJ97evji86dO4tFTOSCl/0bxxIVtqv5xBVYgfLu05UrV6JXr16gdNkl71EauxW4sap+6hXMwUOHDsWYMWOQmZFZpr/B/a3qV41KtJ9hkEp0MzeyIgXu378vamOEh4fLhSuIgsjyBoP0At1xzNsF2sUuIA9oezvjqG9ZECRgTYArjnq7QiewlGsowA1Hxec8JG6iym5X5C6qbxBEx2cYVL8wSApZCFRQSsf8/HwcPnwY/fv3R7du3bB9+3bxnjw4/KTnqsiv9M/QihUrcOXKFfEluaJxlN9nBVgBVoAVYAVYgaopIIVBly9frtoHeWuFU8DPz0/AoJMnT4r/NxWuAXzCrIAKKEA1agkGUbYYfii/Ap9++inGjh0rsmAof2u5haxAWQUYBpXVhN9RIQW8vb2xadMm4TaQl6Ct3MEgOYIx9Q2EGAbJBwySvVcpjRy5gvT09PDOO++ge/fuoH+2b9y4wU6he5XvL1lNpT/T6tUdO3bAxcUFjx8/VqGZgZvKCrACrAArwArUrgIMg2pXX3naO8MgeeoNPhdWoHwFGAaVr4uyvsswSFl7lttVWQUYBlVWKd5O6RQg++/WrVvh6OgoV8Xou/fsgXm7tkAnoMhBwzCmTF2i+oJCi4/sQyctLZHSQxow59fygUNtpIl7kdaUQo4cLAcPHhSrL6m2UFBQEOj9F32O/1Z+/0l1cXNzw7Zt23glq9LNgNwgVoAVYAVYgfpUgGFQfapft8dmGFS3evPRWIFXUYBh0KuoprifYRikuH3HZ14zCjAMqhkdeS8KqEBmZiYWLlwo0kpJA5/y8Pr5F1/glzl/QNvPTW4gSH3BF3k77rhFf+GDjz7C5StXGDC8xHVS1zBIeu+SUyg5ORkzZ86EhoYGFi1ahNzcXLlJAyk9T0V5JcC2fPlyJCYm4unTpwo40vMpswKsACvACrAC8qcAwyD565PaOiOGQbWlLO+XFag5BRgG1ZyWirAnhkGK0Et8jrWpAMOg2lSX9y3XClhbW+PIkSOg4LE8BWb19fXRum0brDx+FLrsDpIbILbB8BQ0O7TDocOHRL0aebpm5PFc6gsGkRZUM4jSxJGr5cMPP8SgQYNAKSHZJfRiF1BF1xGl4DM3N8fDhw/lekznk2MFWAFWgBVgBRRFAYZBitJT1T9PhkHV15D3wArUtgIMg2pbYfnaP8Mg+eoPPpu6V4BhUN1rzkeUAwUoRdyGDRvg6ekpVyniKBhLcGrZ8uV4o0VzDPrkI3z60/f4bNQP/KxHDd77/FM0bd4cc/6ci0uXL8sVPKwogF/f79cnDJK2nWreFBYWYvHixdDU1MSuXbtw7do17r+XuLqk+klfAwICRG010pMfrAArwAqwAqwAK1B9BRgGVV9DRdkDwyBF6Sk+T1VWgGGQavU+wyDV6m9ubVkFGAaV1YTfUQEFKED8119/oaCgQC4Dwzdv3kRgYCA2btyIRYsXi1RXlO5K0Z5z585Fo8aNMX3GDCxUwPOX6r123Vp4eXkxSKgCRJAHGCSFGXQ/WVhYoHv37hg1ahRSUlI4bVwV+pJSxS1ZsgTp6el49uyZCswQ3ERWgBVgBVgBVqB2FWAYVLv6ytPeGQbJU2/wubAC5SvAMKh8XZT1XYZBytqz3K7KKsAwqLJK8XZKpYCHh4dwCVCQWBowlrdXWoVPLiE6R0V9Ul2mpk2bIjo6WqTtUtR20HlTf8jbNSLP5yNPMIh0ohRxVPfmu+++Q58+feDj48Np46oAhA4ePAgnJyc8evRIqeYCbgwrwAqwAqwAK1AfCjAMqg/V6+eYDIPqR3c+KitQFQUYBlVFLcXflmGQ4vcht6B6CjAMqp5+/GkFVeDo0aOwsrLiYHAVgsGvAh7y8vLQrFkzJCQksBOjlrV+lf6pzc/IGwyitlItoatXr4q0cR07dhRuIXmrGVabfVKdfTs4OIgaa/fv31fQUZ9PmxVgBVgBVoAVkB8FGAbJT1/U9pkwDKpthXn/rED1FWAYVH0NFWkPDIMUqbf4XGtDAYZBtaEq71OuFaA0R0uXLkVERAQDiloGFAyD7qmsm0geYZAUhpDTa9++faKO0KFDh4RrTfo3fi3/mo2Li8PKlSsFQJfrAZ5PjhVgBVgBVoAVUAAFGAYpQCfV0CkyDKohIXk3rEAtKsAwqBbFlcNdMwySw07hU6pTBRgG1ancfDB5UODSpUuYP38+Ll68qLKB+roKeDMMKj+wXlf61+dx5BkGkS7kCDIzM0P79u2xYsUK4RiqT73k/djXr18XjqqMjAyuGyQPExmfAyvACrACrIBCK8AwSKG7r0onzzCoSnLxxqxAvSjAMKheZK+3gzIMqjfp+cByogDDIDnpCD6NulOAHEEbNmwQwWB5D8Aq+vkxDGIYRLV65PU6pnPz9vZGz549MW/ePFy5ckVuz1UeNNy2bRsCAwPx5MmTuhuw+UisACvACrACrICCK/Dvv/9CR0cHAQEBePjwoWiNLAx6/PgxzM3N4ezsLL6HKHhzVf70aQHNgQMHkJycjKdPn6I0DKKUxeRQT01NFX9XecFYAFagHhSg/wFpYSDdjzRGl4ZBVPNYT08Ply9froez40PWpALUv/r6+vD09BSZgWjfsjCI/re1t7cHpUWnBaP8YAVUQQGGQarQy9zGEgrY2tqKf8jkOUgtD4HfmjgHhkEMg+T9Prtz5w6CgoIEEFqwYAEDoRekjjx16hQsLS2LA1klBlb+hRVgBVgBVoAVYAUqVMDd3R1du3YVjmRHR0d069YNEydOFPX4tLS0xO8XLlwQQckKd8J/UAgFKPC4fv16tGvXDj/99BOOHz8OqlW5f/9+TJ8+XaQpHj9+PAoKChSiPXySrIAyKkD/ow4bNgzNmzfHzp07MWvWLPTq1QuGhoZ455130KRJE5w5c6YYHiijBqrUJlqMMWDAALRu3RomJiZ49913MXr0aFDKeFqc0alTJ4SEhHAGDFW6KFS8rQyDVPwCUMXm6+rqwsrKChQErgngwfsoCTxIV5pUabVbbm4umjVrhoSEBGRmZopVcrRajjUrqZky6HH37l3QvRUfHy/uLWmaOLoGDh48CErPKK/tpGuWviBSYIbqidEKMXk91/o8L1otRSub79+/r4pTB7eZFWAFWAFWgBV4ZQXIIdKjRw/85z//QaNGjcQrBSGlP1PKWvqOzA/lUCArK0v092uvvYamTZviv//9r3h9/fXX8X//939wdXUFOcL4wQqwAvWnwO7duwWcbdCggYBCNB7Tk+7Xxo0bIzIysv5Ojo9cowrQHDx48GAx/jZs2BA0Nr/xxhuin2lenjt3LvLz82v0mLwzVkCeFWAYJM+9w+dWKwps375dWEQZBtUekCAYRDBg6tSpoC9XU6ZMQYsWLbB161bcvHmTA+0vcF/UZ7C/usemVTYaGhr49ddfQV+yqP/pOiDHjbwDFhoPfHx8QKtzKU3AtWvX+DotdZ36+/tj165dQpdaGZx5p6wAK8AKsAKsgBIr4OTkJNwiBAMo+ERP+rlNmzZi4dSzZ8+UuPWq1TRyB61evRotW7YUfUz9LH2OGTOGg46qdTlwa+VUAYpLvP/++wL+yI7JdK9Sijj6/5AfyqOAND089a9sf7dq1UosDOU5WHn6mlvycgUYBr1cI95CyRSglXe0yoOcDNUNfvPnywdK5AIhR5B09Zv0NT09nXUvFWBXpmuIvlAT/KGVNvQlS9rvlHNZEe43+sLv5eUlgBBBY3axlby/4+LiRGCD/zFSskmRm8MKsAKsACtQJwpQXYIffvhBrESWggF6pVqm9B2KH8qlALmDKO2UbF+T44BdQcrVz9waxVaAMlu0b9+++D4lSECAiGp+8UO5FCB30KRJk4QjSHZcXrJkCS5evKhcjeXWsAIvUYBh0EsE4j8rlwI0AcyfPx8pKSkKEZxWZFCwZ88eYbemiZZcIvSP7o0bNxjAKTEMouuVijMSEJL2+7x58+TeFSR7nxHooHRobdu2hZGRETvZZK5XSvW4cOFCLqypXNMit4YVYAVYAVagDhVwcXEp4Q6iujKJiYlcK6gO+6CuDiV1B0m/F9N347Fjx7IrqK46gI/DClRCgVu3buHjjz8uXsxI9+mJEyfYFVQJ7RRxE19f32JIT+CPspoEBgZyrSBF7Ew+52opwDCoWvLxhxVNASoUSMFpKtgpGwDmn0s6AGpCD3IHUSE+mmTpn6CMjAwGcDKB9ZrQWB73QStbBw4cKOz2lHM5JiZG4fqdxomjR4+K69fT0xP0uzxqXdfndPnyZSxatAiFhYUctFK0yY/PlxVgBVgBVkAuFJC6g+g7EgUd2RUkF91Saych6w5iV1Ctycw7ZgWqpQClhJO6g9gVVC0p5f7DUncQZbGhOZhdQXLfZXyCtaQAw6BaEpZ3K58KkP2TYBCnf6p5+FNeYJrcQeQK+ueff9gVpAIgSHoNGBgYiCK5c+bMwZUrVxQSpNAqsaVLl6Jfv36IjY0Vq8Ok7VPVV3JNLV68WIBdzqksn3McnxUrwAqwAqyA/Csg6w6i7xjkIOGHcipAfUsLaah26o8//oi8vDzlbCi3ihVQYAXo/z5azMiuIAXuxCqcuqw7iGoG8/+1VRCPN1UaBRgGKU1XckMqowA5gggGkXtBVQO6Nd1uqgVz+8YN3Lh4Edeys3AlNQWXEhNwMS4W6UGBGPvZp0jw8kRhbAwuJcTjckqK2O5GYaH4nCLUkqlpzZRlf3du3sTNS5dwPTsbV9JScSkxUfR7TkQ4Jnw1HOGO51EQE42Lot+TcTUrEzcKC3D7+nWFcAtdu3YNVOR36NChoJWdfK3eEzAoLS2NvzRXZsLhbVgBVoAVYAWqrQAF0x88eCAWZdCiBGV4UtrkkSNHisVSBAeUoU3choqvzfj4eLG4yNbWFvTdUlm0ou/FtMqeH6qnwOPHj0UsRVmuZcoCcfDgQXz//feitjT9rixt43aUHZtpDh4/fjyWL18O+r9Wmfr70aNHvMBE9YbkV2oxw6BXko0/pKgK5Obm4q+//mIYVEmXyt3bt3EtKwv5EeFId3FCkrERog8fQOiWjQhYvgQ+f86G++QJcJvwG1zH0/NXuE74De6Tx8Nj2iR4zpgC50m/w2P6ZPE7bUt/F9sVbes2eQK8585CwNJFCN20AVGH9iPR6CzSnc4jLzQEVzMyQNBBWQCKIrTj7p07uJ6Tg4LICGS6ueCCmQmijxxC6KZ/ELBkIXxm/1G238f/CvcJv8GD+nvyeDiNGwN3+nnS73Cf+PvzPpf2+6Tf4TVzGvwXzEfIujWI3LcHCfqnkX7eAXkhwbiSniYX/U5jxocffoiZM2cqrMupJq85stJTzTVeQaWosyCfNyvACrAC8q8ABZijo6Nx6tQprFmzBn///beo+Ul1P5XlOXr0aEyfPl0sUlOWNnE7yr8+aSHiqFGjMGvWLKW5fqV9TbUkt2/fDgJdlIGDvx/K//j6qmdI4Nra2hrbtm0TNUSl14CyvM6YMUMsApw9e7bS3afK0kc12Q6q3zZ16lSlnIPXrVuHM2fOiHqEBIf4wQqUpwDDoPJU4feUVgEK7NI/lOwMKpsm7mZBAXKDg5BsbiYC8wFLFwqoQ+CGgvZBSxYgYsNaxO/ZgZRjh5F55iTyTI1xydYaV12ccNPLA3cC/HA3OBD3QoIqfNLfabtbXh645uqMS3Y2yD9ngiz9U0jRPoKEvbsQtXE9gpYthves6QI0uU36Hf6L/0bE7h24YGqMHH8/3MjLxb27dxkSVRLsVQQEbl68KKBbipUFIvfvFbCHYB71u+eUCfCfPhkh0yYhcupExE6bhKTpk5HyxxRkzJyK7NnTkTdnBgrm/oHCuX/g4p8zX/ik7Wh7+lzmzKlI/WMqkmZMRvz0SYiaOhFh0yYhYPoUeE+dCLfxv8Jt4jj4zZ+DsG2bkWh4Ftk+XsJVRrCqovbUxvuhoaHo0KGDCErRSqLaOIai7HPZsmVITk7mf/aVdpbkhrECrAArUL8K0IKDnTt3isVblJqUVmubmJjA0tJSBCIpGKkMTysrK9BTGdrCbXj5NamMfU335OnTp7F582Zxvy5YsAAWFhaglFv8UB4FKG5iZGQkaqsQ2CQYpK+vL/pa2e59ZbxPla2Paqo9yjgHU5uMjY1x4MABkZ6U7lf6mWp3s4tTecbkmmoJw6CaUpL3oxAKMAySQCCy9VM6t7Tz9gIA+MyfK6CL9+w/ELZmJRL37RZwhkDNbT8fAXbuhwbjfhg9Q/AgPFTyjAjDg4gwPKRnZHjVn0Wfp30U7zMsRBxDHC8kCHf8fXH5vD2yz55B0oG9CF+3Gj5zZ0rOd+4sROzZiRQbK1xOSkRdQwJFCeBLz5P6/WpmBjKcHRF9cD98/5L0u+fkCQiaPlkAH4I9BHryZeDOpT9nQjznzcKlouflebNQ3ad0X+K16BgX/5wlgBKBo6xZ03BhxmTETJuEkOmT4EOQaMJvAk6Gb92Mbc+jHQAAIABJREFUZItzuBgfhzu3btUqoCF7vY6ODrS0tEBgiH6XaqpqrwyDFGKq45NkBVgBVkDhFKB0cN7e3qB5hlK3ODk5obCwUCzgohQu/GQN+BqQv2uAwA/BgqSkJGhra4OAENWMpdTsXAtL4YbhMidMsRNyfpEjg5yaBOup9jL1O9+P8nc/cp9wn9A1QPdnfn4+7OzssGLFCvG9Kjw8HOwSKjPEqfQbDINUuvtVr/GqDINuXb2KbB8fke6L0rKR6yJwyUIk7NmJHCMD3HB3LYY+D6TA51Uhz6uAoRd8phgWESgiKBUSJJxIuSaGAlyFLF8iUpJ5zZ6BqP17kenuBnK8sHPonqjPkxvoj1jtI/D5cxbcJoyD37RJiJg6UbhyyKUjcfTMKgY+1YU8Nfl5CSiSnBudJzmLyJkUPW0igqZNgvvEcfCcPhkRO7cJyHUjP79W+p3+0Z0yZQo+//xzUfxX1SCQtL0Mg1Rv3uQWswKsACtQFwr4+/tj6dKlIpCcnp4uAo3SuYdfyzr6WRPWRJ6uAVpwRt+V3d3dxX28e/duXL58uS6GDj5GLSlAaf8IBNF3/4CAAFHvivpZnq47PhceB/kaqPgaIDBEGT3oPqbvV5R+l+p98YMVIAUYBvF1oFIKqBoMunXtGrK8PBC+Y5uo8eIzZyZitm9BjtFZ3PL1FmBFCn5e2d3zAojzSm6hKuxPuImKANEdPx/kmRkjftd2+P45R8Ahqm2U4eqMm5cvq9QX19s3biDHzxeRu3fCY+pEeE0ej7CpE3FhxhQBUwiqSJ0+NQlu6mpfUicRtYMcRJRujtLM+ZJzaOI4BK9dhTQHO9woKKhRMJSTk4O3335b/FN09epVlbqmpF+0GQap1JTJjWUFWAFWoE4UoFoUq1evxo4dO5CdnQ0OOFYc3JHOx/zKGsnjNUDueQ8PDxF4PHfunHDT18kgwgepUQWePHkisiJQrdDg4GBOsV/NtOzyeK/yOanGHELfpyhN3JYtW7Bp0yZQPJRdmzU6XCrszhgGKWzX8Ym/igKqAoMuJSQIJ4jH9MkgAETuH6rtczc4oDjNW22DmvrYPzmIKI0d1Syi1HKUVs53/hwBRKIO7kNBVCTu3L6ttEH8yykpiDuhC68/psJr8gTh/kmbOVXU81Fk+PMyyCR1DxEcIqcT1TYiMOQ+6XeE79yGvJDgGkklR1+mKIVN27ZtYW9vLyzYqvZFmmHQq8w8/BlWgBVgBViBihSgoISenp6oR5GQkKDSqVhV7TsFt1c5g5G0Gv3kyZOgml+pqalcZ7KiwU+O3w8JCRHppWxtbXHlyhWl/d+ZxyDlHIO4X0v2K8UwyN23atUq0D1N0J4frADDIL4GVEoBZYZBVC8nx88PQetWwW3yeERuWIcCczPcDSoCQHKS8q3OIBHVISoCQxdtLBG9ZSPcp0xAwIqlIo3cnZs3leKLLU3u+aEhCN2wDm6TfkfQ1IkijRo5ZpQZAL0IEBEcKpw7E5mzpgkg5jFxHPwXzEeagz1uX7tWrX6nL09r1qzB4MGDxeplVfuyyTBIpaZMbiwrwAqwArWuALmCKKc9FSjnoGPJAI6qfcfg9ipP/6elpYnAo6mpKQcea30UrdkDPHv2DEeOHMGGDRtAKTvpf02+N5Xn3uS+VM2+JEh/+PBh4Q6iekL8YAUYBvE1oFIKKCMMkkIg/2WL4TVjqnDD3PD4f/beAyqu7MrXf/PezJv3nze227Fttz0z9sybv8M4je3pViRHAcoJRBAoC5AIyjnnHFoJCQnlnBAK5JwEIkgi51g5oVZ3S+976xaUKKJQQE3B1Vq1ClXdOnef37m1z6nz3b33g+YUcAMNAHWRYs4QMaSOjabwwF5iZ3iTMN/PpKGQHgJlpJOyZCFRAvzzctNHxehTwM2dQXewZCC91zhnBrWzfMif6k6sh6s+Uk4PhRSKt/5h09jYyJ///Gd9SpuBli5OhEEDasoUOysqICogKtDrCjx48ID58+eTl5cnRgWJqYjeem0mbnD2rQ1O4XeKsPEopCYSIK/4z3QUEGoFrVy5knPnziH85hG/W33ruyWOhzgeb3sNCCk8Fy5cSHZ2NkIqSPHfwFZAhEEDe/wHXO/7GwySFDwlddVyor09Kdy7G41QBygjDVOs//OhIoaEOkPapASKDx0gZvpUkpcspP7RIwSo9rYT64f+nKyslIwNa4nycCPbc4q+DpAQBTSQIM+b9lWIFhKipR57e+ihUOJ8P336ON1bpg2MiIjQp4sT0sYJd9p86GvgmzqfCIMG3LQpdlhUQFRAVKBXFTh58iTLly9HiBD6puY28bzi5pp4Dbz/a+DatWv6VHHC72+xRkWvutH32riwUSxEayYmJoq1gkRAL87L/egaKCgoYMWKFQj7GAKwF/8NbAVEGDSwx3/A9b6/wCCtUsmTsJNEek4hZ+NaVDFR+pRoHwqo9IfzCFBIgGf52zbr08flHT2Euo/nRBbAReHli/oaSOmeblTN8tGngntTMDKQjzdAoRyhppDbJLJ370RZV0fTG6ZAENLFzZ49GwsLC+qEz/ejhWJ3fRFh0ICbNsUOiwqICogK9KoCQvTA5s2bqa+vHzBzaXfzrPje+4cSoqbfjKbCXehBQUEUFhaKdYN61Yu+38bj4+MJDg4mJydHjNYcIL/vRB/5zfjID627sGchpH+8dOmSvvbx+/UcYmumpoAIg0xtxER730mB/gCD5BUVpCxdTILfbOqvXtKngxMjgTJ5W0AlRFI13rpO8nxfkhYEIikseGMw8CEmcgFYpK1ZSaynm74mkJAObiBDnXftu5A+rnKmN8mebsT7zqI+J+eNo8MqKyv51a9+xbZt21C8Q9q5D3H9vK9ziDDonaYg8cOiAqICogKiAu0U2LNnj34eFWHQwNiMel/rEbGdvn+9CNHzAgx6+vQpL168aPfNF//bVxUQxk2AQfn5+SIMEmGQeJNGP7oGhBSQa9as4cKFC/qov77qg0S7PowCIgz6MDqLZ+kjCpg6DKp7lE3MrGlkr1mJOj4WIbrlbSGI+LlWgCTUFBJSx+VsWk/09KnUpCS/MRjozR+lQjpAAVikeLpRLUQDiTWB3hsIq589jWwvN33dpYroSN4kbZwQXn3q1Cl+8Ytf8OTJkwFRYFWEQX1kMhPNEBUQFRAV6CcKiDCo72/q9+YaV2y7/45/XFycCINM0E8bwyDht474He2/31FxbAfW2IowyAQdci+aLMKgXhRXbLrvKWDKMKg26yHRPp4U7N2NLiVJhEBZrTDnfYGtpvRUSg4dJHqqB9WJCX0CCDU+eULsTB+yPN30NW/eNSJG/PyMDiBJqLf01NudqCmTKL0XgU6t7vEPH7VajZWVFf7+/voCuf19US3CoL43r4kWiQqICogKmLICIgwaWJtR/X2dJPav9XoWYZBpemYRBrVew+L3WdSiP10DIgwyTZ/cW1aLMKi3lBXb7ZMKmCoMEiJDon28KDqwl6a0FBEE9QIIMgAlIdqqLOSIHgjVZWd9o0BIXl5OrBAJ5uVGvZgWrgPEeZ9gS4i2KvTxIMptIlUCCNRoegyEhJzoH3/8MUlJSf0+nYIIg/rk1CYaJSogKiAqYLIKiDBI3GzrT5ttYl9ar2cRBpmmWxZhUOs1LH6fRS360zUgwiDT9Mm9ZbUIg3pLWbHdPqmAKcIgtVxGQoA/+Vs3iSCoFyGQAQYJzwIQerp7B/F+s1HW1HwjNYS0ajUpSxaS7imCoPcJfbprSwBC+VOnEOPtgaykuMfjrtVqmTJlCuPHj6exsbHHEMkUF5ciDOqTU5tolKiAqICogMkq0N9hkLw4lfsRt7hx40bXj5vhpBQ2otR0vfGmlRSQdD+CmEfVyFW9n7pJ21hISvgZDu/azPq169i0/QAnbyZT1KBEo+vazm92baOl6lE098JvttX65k1uhd8nIbuEBpUGnVENCHlxGg/uRpJeIulW/7fql7aRp0kPuBuXQ41c1ea8b9Wekd2m8HkRBpmmWx4oMEhekk7U3dttfUV7P30znKSn9cjVXftcrbSI1Kh7xGRXIlVqe/13oLaxmLS75zi6Zwsb1q1j47b9hF5P4GmtHHVf9c3aanLjH3DnZlvffPPmLcLvx5NVVIeijcZySjOiuR+ZSlGDHM179n3CfJoSdY+4Rx9mzPqKvxZhkGn65N6yWoRBvaWs2G6fVMAUYVD+8WMkB/ihSUoQI4I+EAwSgNCztBRSFwSQ+/kBtCrVe1nYHTt2TF9bRoAHr1sUFF66QJzXFGpmT+vViJju4MhAfK9hznTSPFzJ3LgOjVz+2nEyjGNeXh4//vGPuXbtGkLqOMPr/e1ZhEF9cmoTjRIVEBUQFTBZBfo3DNJSfmkxY+zMGDRoUNePwdYEnn1CtazrDUdNwWmCxjjgtSOOiobeXGcoKbq/n8DJ9pgNGcRQCzucnJ2wtxzG4MFDsBwfzLG4EiQfAEi9+RpKQewOD5yshnSq9eDBw3H0WselrGpkevCmpezSEsY5jmP59SfUKLrW/81taaJJncuJeWNxnL6XpNIG1O95Q/OtbPqANogwyDTd8sCAQVoqrq3C1cmiU1/xyl8PtmLeySzKGrv2DZqi8yyd7IzXlnsU1ih78TegkpKowyzyGIHF0MEMMbdlhOCbrYYzRPDNYwP4PPIp9YrX7zN8cN+hTGDfzDHYDulsHhzM4OEOeKw6S2q5pNlPasu5tnoKI8cv4UJ2OdL3DLnU+adZONmFaTuieFLzfvZ5Primb+HLRRhkmj65t6wWYVBvKSu22ycVMDUYJESlRAqRIVcv88XDjNfDoIwr7DD/Nyatuo0yo+vj1ftd+dMnP+ZPE3dQlpLGF68gSyynR/2S/3TdjSItjZ4eZxxV8zzjJvusfonD/LOo0ru2oc1nXp0/k2eRa3H9P2ZsOR+P6sYSRv9yKJvOxqPJfP81grqywfC65PZNojxckRYV9jhKpLuFwPnz5/n+97/P5MmTEeBBV1BILZUSPX2qPm2ZEK3yzlBmUQhfVGej3eSP1Pc9tGds0+qjPHtSxYvqWNSLfN/B1plI12xFuXUJEr+Z79DOu/evapYPUZPH8yZpAoUCq4sXL2bIkCHU1tb24g+Bb/ZuXBEG9cmpTTRKVEBUQFTAZBXo3zBIQ/mFBThbDsN9WwQJGbn69Z+wBmz7eExpvQqNtus5XvP0JPOcrZiyNYbyXoNBWsrubsbbyZxhDtPYcCaOx1UNSOVy5A2lpF7cgLejOWbOwZzJrGoBKl3b3N2auHfeUxCzxRU785EsPh1FSnaLzrnZZCSEc3ylFw5mQ7H1O0FGhRRNk4bS8wsYaTOSJVefUC3vesP3rexV5xAyxwWbqbtIKKkXYZDJeqmBZfiAgUFXljHezgy3jdeJTu3GN9cpUXfnmwvPsnCcPW4bIiio7i0YpKXiwU5mj7FimL03a0IjySmvRyqTI28sI/3qFmaOtMLcaT6hySU0qvuSX26iSRnPbm9nrFyCCYlIIKtlDszNziTxzknWTnPGYpgNcw4nUlSvoUlbzuXlE3EctYBzWWXvHwblnSJwvANeWx/wuFqEQQPLw4m9NSggwiCDEuLzgFDA1GBQwYXzpC4M6nl6uLgNTPjhvxNwJLZbeKLa4cIn/+tv+R9//0vm7LuHJtMAbWI4avER/+SyFXlqGj09zgBPmp9jSduxnLNh99G9arfnIOfZvaU4fvcPrDgdhybmDBdXbCExMoWmhz1vo609b/85AcBlLFnI41OhCGnb3uqHoNFdGyqVim9/+9v87d/+Ld/61reYNm0a+fn5HaBQacQd4r09qHtfUUFLTvFcVYRuy/uGQTORX8rnZV0qmn2rkfm/C8SZjSK8hBfZZ5AFzvlGYZAA3zI8XcnZuwuNQtHjcRf8yy9+8QtOnTqFMNbver30xc+LMGhATJViJ0UFRAVEBT6YAgMDBg1n+pEcKhs1CDePdPowWi82z/9qpLXllJbXIdNoUXcDgzTSWioqa5Eota9JRaZF2VBFeUWtPp2Rcbo04ZzamntscndguKU3e6IKqFUI9ho2FXVoNXLyzwQzxsqMMWvDKaxTtj2fRk5teRnlNRKU2teBFRWSqnLKKhtRaHRt2+mghcGG1z0bYNA41t0rpFahfaW1VqtB3RjDDk8HzK0DOJNfhUzXPQzSyGspLy+nRqJE+0qHTmxQSakuL6W8TobGuN8iDCIoKIinT5/y4sWLD+ZTxBO9mwIDDQZ570+jqK4b39zhu69GVldBmeBH1RrU3cAgjayOyspaGhWa7n1IkxZlYzUVFTVI2vjdZn+jrY1kxzQXLCw92X4nl0pZR9/85MJSJtlZMHrFVXIr5W19qkZOXUU55TWCv33dPKFCWl1BeWUDcvXrju3EH3bmvw0waOxKbuRUINca5kEtWo0aSfw+Pcyy9jtBZkkjutfAII28jgrBNzcq2vrc9udWSampKKO8VopK0zonqUUYpN+reDdPIX7a1BUQYZCpj6Bo/xspYGowKHXFMooPHdDXsOkJ4NCETOE/PrYjNDIZnRBtkxFD8aUzPLqbyDMjmCJAnp/+439g+9cf8J3feBARnUKTPjqnGQb93AgG9eS4trbFkbVvPVfPRephkOTyAc4eukz11f0c9PMmeNo8Qk/dNQJQKRQdW8PG6V4s9l9D1Jlg7D5qgUGxF7m5YRdp0SnN9mdGk7FjEau8PVnst5zw63FGkCiJgqNr2Dzdi4VzlnD1YpQRjEqlMmwT22d5EegzlwP7ryLJyDCKiOoaGJWFHCF58YI3ggLdbeKfPHlSD4T+23/7b/zDP/xDp1Aoc8sGHk11R0hZ1mVU0LwlKEIvort7B93ZfcgXzjY6dg6yHSFow++gu3YcxeZzPFcXt8Ig/wUoTgifvYX2+CZk63ag2rsKiQHmdNu2IfJmJtJ1e9Cm1fGyLh1d2A5k82YhmTsX2S7h3OHoboahWh+IxNcIEi1ej+r8VXT3I9BdOopiua/+fenGfWizJLysTEITshnp/ADk+w6h2hDQGim0aC3Ko9uRzZ+NZP5ylEd2Il+yEtW5q2j2L0Iq2P/K9nB051p0MURD+QWhOHYO3d276G6EoVw7D4nhPeNop7kz9FFZ8TO8EaK0uhtP4/eEzZ1Vq1YxdOhQ6urqevw54zb6+t8iDHqjKUc8WFRAVEBUQFTgNQoMHBiUS5W0dTOq6/leS13GeTbOHout+TCGDBmGzbh57D25GW87S9yMIoO0jdlc2ebHBDtzhg0dgpmjO4uPXiVkxXS8l4VRVm1Yw2ipe3iJ7fMm42g5nKFDhmDu6Mb8HVd5WGWoM6Gl5MJCRlsNZdSq2zytUbTdTGzZZNPVx3BoyTx815wju1qIrhHu+i4m6ugKpo+xwWzYEIYMMcfBNYBdN3OpkbfU6NGWcmPTXHwWHyPi6n6WTRuNjdlQhgwZjt3khYQklCJRaaiK2M7cqVPw/zyR8kbjG2s0PLm4mtme09hyu5C6TlO6tcKg9feKqGufyk6dyecznLC0C+Lc465hkLI4mmMrpzPG1oxhQ4YwxNwB1/m7uJFTg9xoQ7FJUciDQ0vxHmmN2dAhDBlqzgj3BXx+v6A5VVOnMEhKzqUN+HpNYc62W+TWyNG238DsJ/8X08S9xvn10bcHGgzy2Z9GcX1PUqtpqcu6zDb/iThYDGPokGFYj/FlZ8gmZjjbtokM0jbmcmN3IK6OlgzX+2Y3Fhy8yNHVc5i+7DhPyxpafidqqcu+zp5gd5ytzPS+2cxxMn5bLpJWJkGlB1FaSq+uYJL9cFyWXSG7XNa5b26IJ2RVEH6rwkgva0lLqSwhNnQts8fbYzFM8Ldm2E/yZ/vVh1RI1c3taMsI3xHAjMWHuHnpAKtnjcNOP/cMx3ZiEIejCqhXaqi+v4fA6R747YumsNb4RgANBVc34u89nY3X8qjsbJ57BYNWcTO3EkU7wKbOOob/WBvs5oXysLRrGKQsjePkujlMsLfQz3lDzOyZ5L+NK5kVSIxrDimKiA5ZxYwxtpgLc9JQcxzdAtl3J49qmYZOYZA0j2vbAvBxn83mqw+plGr63W95MU1cH3W635BZIgz6hoQXT/vNKGBqMCjK26PnKeKy0ihY8Ck//Is/2Ump6B4cZaPlr/nriMUkRqV2hEHfMmfP8QDsPv4un00/RkNaOs+zOoFBPTiuDQzKuMq2z37I4NlhKNNSyQ/6jB//8rf85d/+wOhRY3D99J/53g8+ZfO5eLSZaTxd58Cvv/8Jf7VwxtPmj/z6lz/jx//<KEY>//<KEY>//<KEY>//<KEY>//<KEY>//<KEY>//<KEY>)<jupyter_code>def edge_conv(points, features, num_points, K, channels, with_bn: bool = True, 
              activation='relu', pooling='average', name='edgeconv'):
    r"""
    Implements Edge Convolution

    Inputs
    ------
    points: (N, P, C_p)
    features: (N, P, C_0)

    Parameters
    ----------
    K: int 
        Number of neighbors
    in_channels: int 
        Number of input channels
    channels: tuple 
        Tuple of output channels
    pooling: str
        Pooling method ('max' or 'average')
    with_bn: bool
        With Batch Normalization, if not then we use bias on Conv layers

    Returns
    -------
        transformed points: (N, P, C_out), C_out = channels[-1]

    Warning
    -------
    This function assumes that the Data Format is Channel Last, do not change
    format in the DataSet class.
    """

    with tf.name_scope('edgeconv'):

        # Compute the K-nearest neighbors to each particle in the same jet (on eta-phi or latent space)
        D = batch_distance_matrix_general(points, points)  # Get the distance of all particle pairs: (N, P, P)
        _, indices = tf.nn.top_k(-D, k = K + 1, name='%s/topk' % (name))  # Take the particle indices which are closest to us: (N, P, K+1)
        indices = indices[:, :, 1:]  # Remove self from the list of nearest neigbors: (N, P, K)

        fts = features
        knn_fts = knn(num_points, K, indices, fts)  # Gives the filtered Feature Vectors: (N, P, K, C)
        knn_fts_center = tf.tile(tf.expand_dims(fts, axis=2), (1, 1, K, 1))  # My Features repeated K times: (N, P, K, C)

        # (My Feature, Neighbor - My Feature) for all (N, P, K): (N, P, K, 2*C)
        knn_fts = tf.concat([knn_fts_center, tf.subtract(knn_fts, knn_fts_center)], axis=-1)
        x = knn_fts # We shall call this our X (input feature vector for the Edge Conv model)

        for idx, channel in enumerate(channels):
            # Features after expanding is of shape (N, P, K, 2*C) and the convolution gets applied
            # on a P * K sized image going from 2 * C input channels to the output channels
            x = tf.keras.layers.Conv2D(
                channel,
                kernel_size=(1, 1), 
                strides=1,
                data_format='channels_last', 
                use_bias = False if with_bn else True, 
                kernel_initializer='glorot_normal',
                name='%s/conv%d' % (name, idx)
            )(x)
            # Batch normalize and apply Activation
            if with_bn:
                x = tf.keras.layers.BatchNormalization(name='%s/bn%d' % (name, idx))(x)
            if activation:
                x = tf.keras.layers.Activation(activation, name='%s/act%d' % (name, idx))(x)

        if pooling == 'max':
            fts = tf.reduce_max(x, axis=2)  # (N, P, C')
        else:
            fts = tf.reduce_mean(x, axis=2)  # (N, P, C')

        # Shortcut - Important: It's a convolution to correct dimentions, not a straight addition
        # Features after expanding is of shape (N, P, 1, C_0) and the convolution gets applied
        # on a P * 1 sized image going to the output channels
        sc = tf.keras.layers.Conv2D(
            channels[-1],
            kernel_size = (1, 1), 
            strides = 1, 
            data_format = 'channels_last',
            use_bias = False if with_bn else True, 
            kernel_initializer = 'glorot_normal', 
            name='%s/sc_conv' % name
        )(tf.expand_dims(features, axis=2))

        # Batch normalize and apply Activation, output with Skip connection
        if with_bn:
            sc = tf.keras.layers.BatchNormalization(name='%s/sc_bn' % name)(sc)
        sc = tf.squeeze(sc, axis=2)
        if activation:
            return tf.keras.layers.Activation(activation, name='%s/sc_act' % name)(sc + fts)  # (N, P, C')
        else:
            return sc + fts
def _particle_net_base(points, features=None, mask=None, setting=None, name='particle_net'):
    r"""
    Parameters
    ----------
    points: tf.keras.layers.Input, (N, P, C_coord)
        keras input layer for the point
    features: tf.keras.layers.Input, (N, P, C_features), optional
        keras input layer for the features
    mask: tf.keras.layers.Input, (N, P, 1), optinal:
        keras input layer for the mask
    setting: _DotDict()
        stores some Hyper parameters in a free object
    name: str
        name of the model to be returned (prefix of all parts of this model)

    Notes
    -----
    These are the features included in each of the input facets:
    'points': ['part_etarel', 'part_phirel'],
    'features': ['part_pt_log', 'part_e_log', 'part_etarel', 'part_phirel'],
    'mask': ['part_pt_log']

    Points exist on the Eta-phi plane, we can choose to augment the information
    with the log of energy and transverse momentum
    """

    with tf.name_scope(name):

        # If explicit features are not passed, then the features are the points
        if features is None:
            features = points

        # Makes a boolean mask out of the input mask, 1 if it's valid (in the mask), 999 if it's not.
        # Represents a padding particle or Transverse momentum is 1
        # TODO: Find out if an why Transverse Momentum = 1 needs to be checked.
        if mask is not None:
            mask = tf.cast(tf.not_equal(mask, 0), dtype='float32')  # make valid positions to 1
            coord_shift = tf.multiply(999., tf.cast(tf.equal(mask, 0), dtype='float32'))  # make non-valid positions to 999
        else:
            raise NotImplementedError('Implement an allow all mask')

        # CHANGES: fts = tf.squeeze(tf.keras.layers.BatchNormalization(name='%s_fts_bn' % name)(tf.expand_dims(features, axis=2)), axis=2)
        fts = tf.keras.layers.BatchNormalization(name='%s/fts_bn'%name)(features)

        # Extracting CONV Parameters (K, (C1, C2, C3)) from the settings we passed in
        for layer_idx, layer_param in enumerate(setting.conv_params):
            K, channels = layer_param
            # Masked out points get a (+999, +999) - coordinate shift on the 
            # eta-phi plane, after each layer
            pts = tf.add(coord_shift, points) if layer_idx == 0 else tf.add(coord_shift, fts)
            # Now we repeat edge-conv over all layers starting with (Points, Features) and
            # then the next layers are (Previous Layer masked, Previous Layer)
            fts = edge_conv(pts, fts, setting.num_points, K, channels, with_bn=True, activation='relu',
                            pooling=setting.conv_pooling, name='%s/%s_%d' % (name, 'edge_conv', layer_idx))

        # Filter out the masked out particles to 0
        if mask is not None:
            fts = tf.multiply(fts, mask)

        # Output shape: (N, C) - Takes the average of all the particles in the jet
        pool = tf.reduce_mean(fts, axis=1)  

        # Extracting CONV Parameters (Units, Dropout) from the settings we passed in
        if setting.fc_params is not None:
            x = pool
            for layer_idx, layer_param in enumerate(setting.fc_params):
                units, drop_rate = layer_param
                # Add the Dense and Dropout layers
                x = tf.keras.layers.Dense(units, activation='relu')(x)
                if drop_rate is not None and drop_rate > 0:
                    x = tf.keras.layers.Dropout(drop_rate)(x)
            # Final Classification layer
            out = tf.keras.layers.Dense(setting.num_class, activation='softmax')(x)
            return out  # (N, num_classes)
        else:
            return pool
class _DotDict:
    pass<jupyter_output><empty_output><jupyter_text>Putting together the Final Models<jupyter_code>def get_particle_net(num_classes, input_shapes, model_name='particlenet-lite'):
    r"""
    ParticleNet-Lite model

    Parameters
    ----------
    num_classes: int
        Number of output classes
    input_shapes: dict
        The shapes of each input (`points`, `features`, `mask`)

    Example
    -------
    get_particle_net_lite(2, {'points': (100, 2), 'features': (100, 4), 'mask': (100, 1)})
    """
    assert model_name in ['particlenet-lite', 'particlenet-full']

    # Creates the Setting dictionary for configuring up the Network 
    setting = _DotDict()
    setting.num_class = num_classes
    if model_name == 'particlenet-lite':
        setting.conv_params = [
            (12, (32, 32, 32)),
            (12, (64, 64, 64)),
        ] # conv_params: list of tuple in the format (K, (C1, C2, C3))
    elif model_name == 'particlenet-full':
        setting.conv_params = [
            (16, (64, 64, 64)),
            (16, (128, 128, 128)),
            (16, (256, 256, 256)),
        ]
    else:
        raise NotImplementedError('Unknown Model Name')

    setting.conv_pooling = 'average' # conv_pooling: 'average' or 'max'
    if model_name == 'particlenet-lite':
        setting.fc_params = [(128, 0.1)] # fc_params: list of tuples in the format (C, drop_rate)
    elif model_name == 'particlenet-full':
        setting.fc_params = [(256, 0.1)]
    else:
        raise NotImplementedError('Unknown Model Name')

    setting.num_points = input_shapes['points'][0]

    # Make the Keras input layers and send them and the settings to the net generator
    points = tf.keras.Input(name='points', shape=input_shapes['points'])
    features = tf.keras.Input(name='features', shape=input_shapes['features']) if 'features' in input_shapes else None
    mask = tf.keras.Input(name='mask', shape=input_shapes['mask']) if 'mask' in input_shapes else None
    outputs = _particle_net_base(points, features, mask, setting, name='particle_net')

    # Returns the final Keras model
    return tf.keras.Model(inputs=[points, features, mask], outputs=outputs, name='ParticleNet')<jupyter_output><empty_output><jupyter_text>Running the Model  Setup for the Model<jupyter_code>model_name = 'particlenet-lite'

num_classes = train_dataset.y.shape[1]
input_shapes = {k:train_dataset[k].shape[1:] for k in train_dataset.X}
model = get_particle_net(num_classes, input_shapes, model_name=model_name)
batch_size = 1024 if 'lite' in model_name else 384
epochs = 30

def lr_schedule(epoch):
    return 1e-3 if epoch <= 10 else 1e-4 if epoch <= 20 else 1e-5
model.compile(loss='categorical_crossentropy',
              optimizer=tf.keras.optimizers.Adam(learning_rate=lr_schedule(0)),
              metrics=['accuracy'])
model.summary()
tf.keras.utils.plot_model(model, 'model.png')<jupyter_output><empty_output><jupyter_text>Training the Model<jupyter_code>train_dataset.shuffle()
model.fit(
    train_dataset.X, train_dataset.y,
    batch_size=batch_size,
    epochs=10,
    validation_data=(val_dataset.X, val_dataset.y),
    shuffle=True,
    callbacks=[
        tf.keras.callbacks.ModelCheckpoint(
            os.path.join(PROJECT_DIR, 'checkpoints/best.h5'), 
            monitor='val_loss', verbose=True, save_best_only=True,
            save_weights_only=False, mode='auto', save_freq='epoch'
        ),
        tf.keras.callbacks.LearningRateScheduler(lr_schedule),
        tf.keras.callbacks.ProgbarLogger('steps'),
        tf.keras.callbacks.TensorBoard('logs/', histogram_freq=1)
    ]
)
model.save(os.path.join(PROJECT_DIR, 'model-full.h5'))<jupyter_output><empty_output><jupyter_text>Performance Evaluation  Training Stats Plot<jupyter_code>fig, ax = plt.subplots(1, 2, figsize=(10, 5))

ax[0].plot(model.history.history['loss'])
ax[0].plot(model.history.history['val_loss'])
ax[0].legend(['Training', 'Validation'])
ax[0].set_xlabel('Epoch')
ax[0].set_ylabel('Loss')
ax[1].set_title('Loss')


ax[1].plot(model.history.history['accuracy'])
ax[1].plot(model.history.history['val_accuracy'])
ax[1].legend(['Training', 'Validation'])
ax[1].set_xlabel('Epoch')
ax[1].set_ylabel('Loss')
ax[1].set_title('Accuracy')

fig.suptitle('Training Curves')
plt.show()<jupyter_output><empty_output><jupyter_text>Running TensorBoard<jupyter_code>import shutil
shutil.rmtree(os.path.join(PROJECT_DIR, 'logs'))
shutil.copytree('logs', os.path.join(PROJECT_DIR, 'logs'))
os.listdir(os.path.join(PROJECT_DIR, 'logs'))
# Load the TensorBoard notebook extension
%load_ext tensorboard
%tensorboard --logdir logs<jupyter_output><empty_output><jupyter_text>Visualizing the Network  Who is in my Neighborhood<jupyter_code>topk_layers = [ layer.output[1] for layer in model.layers if 'topk' in layer.name ]
aux_model = tf.keras.Model(inputs=model.inputs, outputs=model.outputs + topk_layers)
# Plotting part
fig, ax = plt.subplots(3, 2, figsize=(10, 15))

for j, axis in enumerate(ax):

    event_id = np.random.randint(len(val_dataset['points']))
    trace_id = np.random.randint(20)

    event = {
        'points': train_dataset.X['points'][event_id : event_id+1], 
        'features': train_dataset.X['features'][event_id : event_id+1], 
        'mask': train_dataset.X['mask'][event_id : event_id+1]
    }
    knn_group = np.concatenate(aux_model.predict(event)[1:], axis=0)
    mask = tf.cast(tf.not_equal(event['mask'], 0), dtype='float32') 
    n_particles = int(tf.reduce_sum(mask).numpy())

    print('Class of Sample', j, ':', np.argmax(model.predict(event)))

    for i in range(len(axis)):
        neighborhood = knn_group[i, trace_id]
        colors = ['k' if particle not in neighborhood else 'r' if particle != trace_id else 'g' for particle in range(n_particles)]
        sizes = [3 if particle not in neighborhood else 20 if particle != trace_id else 40 for particle in range(n_particles)]
        axis[i].scatter(
            event['points'][0, :n_particles, 0], 
            event['points'][0, :n_particles, 1], 
            color=colors, s=sizes)
        axis[i].set_title('Sample %d - Layer %d'%(j + 1, i))
fig.suptitle('Neigborhood Plots')
plt.show()
mask = tf.cast(tf.not_equal(event['mask'], 0), dtype='float32')  # make valid positions to 1
n_particles = int(tf.reduce_sum(mask).numpy())
print('Number of Particles:', n_particles)
coord_shift = tf.multiply(999., tf.cast(tf.equal(mask, 0), dtype='float32'))  # make non-valid positions to 999
pts = tf.add(coord_shift, event['points'])[:, 0:n_particles, :]

D = batch_distance_matrix_general(pts, pts).numpy()  # Get the distance of all particle pairs: (N, P, P)
sns.heatmap(D[0, :, :])
plt.show()
# _, indices = tf.nn.top_k(-D, k = K + 1, name='%s/topk' % (name))  # Take the particle indices which are closest to us: (N, P, K+1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09338764601922496, 'fraction_numerical': 0.058041294883400145, 'mean_word_length': 4.815868638756581, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 86, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '595435', 'license': 'MIT', 'n_tokens_mistral': 22873, 'n_tokens_neox': 21254, 'n_words': 2978, 'path': 'code/ParticleNet - Initial Replication.ipynb', 'repo_name': 'Blizzard57/jet-tagging', 'size': '432395'}