exploration_iteration.md

Data Exploration and Iteration

The first script we will run will be in the scripts directory, called data_exploration.py. From your repo directory you can run

python scripts/data_exploration.py

This will run this script tutorial-style, going through different lines print lines. We'll also go through the lines here. If you want to skip the tutorial style, you can add the option -s to skip the prompts and print out everything at once.

First, we will open the .h5 file

    f=h5py.File(filepath,"r")

.h5 files are very performant when reading from disk, and so are widely used in ML. The code will open the file and print out the different keys, which are variables labelled with a name.

Next, we'll look at the number of events. We print out the shape of the 'labels' variable

print(f"Number of events: {f['labels'].shape}")

Next, investigate the shape of the PMT variables, which is the detector variables we will eventually use to train our networks. We print out the shape of the 'event_data' variable

print(f"Shape of the barrel-only data: {f['event_data'].shape}")

We have 900k simulated scattering events here! labels are 0, 1, 2 for $\gamma$,electrons and muons respectively. The 'event_data' contains only the barrel portion of the tank which has been 'unrolled'. The first dimension (900k) enumerates over the events, the second two dimensions (16,40) enumerate over the row and column in the module 'grid'. Finally last two dimensions enumerate over the PMT within the module (again there are 19 in each mPMT module) Note: the first 19 entries correspond to charge collected on a given mPMT and last 19 correspond to the time.

We can now take a look at some details about our data. Note that the object returned by the subscript looks like an array - we can subscript it and even do fancy indexing:

#f['event_data'][[42,1984],:,:,:]

In fact the object is not a numpy array -it is a hdf5 Dataset object - the data itself lives on disk until we request it

We print out the type of the 'event_data' variable

print(f"The data type is: {type(f['event_data'])}")

Next, we see the size of the dataset will make it difficult to load all at once into memory on many systems by checkign the size of the PMT data

print("Size of the bulk of the data is {:.1f} GB".format( (f['event_data'].size * 4 / (1024**3)) ))

One important feature of the dataset it is uncompressed and contiguous or 'unchunked', as we can see when we look at the chunks and compression

print("dataset chunks: {} compression: {}".format(f['event_data'].chunks,f['event_data'].compression))

The dataset has been prepared as contiguous and uncompressed so that we are not obliged to load it all in memory but we can access it very fast. BUT it will take more spave on disk. In the next section we will see an example of how to deal with datasets with these sizes.

Using Pytorch Dataset object

Let's import and create a Dataset object - you are welcome to look at the source code in utils/data_handling.py. We've added a more detailed walk-through of data handling in the comments of data_exploration that you can read at your leisure. Since this tutorial is more about networks than data handling, we leave it as optional.

The data handling utility has a class named WCH5Dataset that we now instantiate, it should have the same length as the file we were looking at earlier

dset=WCH5Dataset("/fast_scratch_1/TRISEP_data/NUPRISM.h5",val_split=0.1,test_split=0.1)
print(f"Length of dataset object: {len(dset)}")

Next, let's get some random event and label from the training dataset, and plots an event display of this single event using a plotting script

event, label, energy=dset[dset.train_indices[1984]]
print("Label {} and energy: {} (MeV) ".format(label,energy))
#Make some event displays, save them to disk
event_displays(event, label, plot_path='plots/data_exploration/')

The default place to store the plots you made are in plots/data_exploration, so you can check them there. You may have to make the data_exploration directory using mkdir, and will have to do the same for plot-making in subsequent tutorials.

To look at plots, you can navigate to them in the explorer if you are using an IDE such as VSCode. If you are using terminal, and have X11 forwarding enabled, you can use a number of plot display commands. One example is feh, which you would use by:

feh plots/data_exploration/all_mpmt_charge.png

One last option would be to use scp to download your plots to your local computer to view them there.

There should be 3 different event display plots, do you understand what they all show?

Making plots of your data

Always try to learn as much as possible about the dataset before throwing ML at it. Let's quickly histogram the charges. We won't load the full dataset but taking few thousand should be fine (since we have 12k PMTs). We'll prepare the dataset, then finally use a plotting function (plot_pmt_var)[https://github.com/felix-cormier/HK_ML_tutorial/blob/trisep_dev/scripts/data_exploration.py#L214]

#Using plotting functions
plot_pmt_var(charge_data_to_plot, charge_labels_to_plot, charge_colors_to_plot, bins = charge_bins, xlabel = 'PMT Energy (photo electrons)', plot_path='plots/data_exploration/all_mpmt_charge.png', do_log=True)

Can you do the same for time variables, rather than charge? You can work directly in the script (and use option -s to skip tutorial mode). Hint: charge is indices 0-19 (one for each mPMT), then time is the indices after that...

Let's also plot the total energy in the event and also the true particle energy. Again we'll prepare the dataset, then finally use a plotting function (plot_pmt_var)[https://github.com/felix-cormier/HK_ML_tutorial/blob/trisep_dev/scripts/data_exploration.py#L239]

plot_pmt_var(sum_charge_data_to_plot, sum_charge_labels_to_plot, sum_charge_colors_to_plot, bins = sum_charge_bins,xlabel = 'PMT Energy Sum (photo electrons)', plot_path='plots/data_exploration/sum_mpmt_charge.png', do_log=False)

Muons seem to look very different - but is that expected?

This is the end of the data exploration script. Are there other plots you could make? Either event-level or PMT-level? You can work by running the script with the -s option to skip the tutorial step-by-step instructions.

Data Iteration

This is a short script to walk you through how to use PyTorch DataLoader objects to iterate through the data. It is not necessary for the rest of the tutorial, but might be good practice if you wanted to try your hand at any pre-processing while running training scripts.

To the run the script do

python scripts/data_iteration.py

Here we load DataLoader from PyTorch

from torch.utils.data import DataLoader

then define train, validation and testing DataLoaders from our dataset

train_iter=DataLoader(dset,batch_size=64,shuffle=False,sampler=SubsetRandomSampler(dset.train_indices),num_workers=2)
val_iter=DataLoader(dset,batch_size=64,shuffle=False,sampler=SubsetRandomSampler(dset.val_indices),num_workers=2)
test_iter=DataLoader(dset,batch_size=64,shuffle=False,sampler=SubsetRandomSampler(dset.test_indices),num_workers=2)

You see the parameters - like batch_size and sampler - the sampler uses the indices we computed for the training, validation and testing set - if you use a sampler shuffle has to be False. On each iteration the DataLoader object will ask the dataset for a bunch of indices (calling the getitem function we coded earlier) and then collate the data into a batch tensor. The collating can also be customized by providing collate_fn - but for now we will leave it with a default behavior. Did you notice the num_workers argument? if >0 this will enable multiprocessing - several processes will be reading examples (if supplied applying the augmentation transformation) and putting the data on queue that would be than 'consumed' by your training/evaluation process.Your 'instance' has 2 CPUs for the job so we will use that. We are beating on the same storage with all threads - so if we aren't doing much preprocessing it doesn't make sense to make this too high.

Now convince yourself that the data and labels are already tensors - which we could plug into our future model - let's iterate over first 40 batches:

num_iterations=40
trecord = loop_over_set(train_iter, num_iterations)

This calls the loop_over_set.py function defined at the top of data_iteration.py. This loop is technically all we would need to train, as we load a batch of training data from memory, we could train over it.

By the way, do you notice that roughly every 2nd iteration the time it takes to give a batch is huge? Why?

What if we want to pre-process the data? We can do that directly in the loop over the train iterator Can you modify loop_over_set function so that time is centered around 0 with a standard deviation of 1? Input data being arounds this range can help converge faster. Can you make plots using some of the plotting utitilies in data_exploration showing the difference in time? Are there any other pre-processing steps you could think of to do with our data?

As with the previous script, you can use the option -s when lauching it to skip 'Press enter to continue'.

Next let's learn how to train a model.