{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport plotly.express as px","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-03-01T19:55:15.960924Z","iopub.execute_input":"2023-03-01T19:55:15.961695Z","iopub.status.idle":"2023-03-01T19:55:19.543490Z","shell.execute_reply.started":"2023-03-01T19:55:15.961651Z","shell.execute_reply":"2023-03-01T19:55:19.541756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Paper Overview: Graph Neural Networks in IceCube\n\n","metadata":{}},{"cell_type":"markdown","source":"This notebook is based off [Graph Neural Networks for low-energy event\nclassification & reconstruction in IceCube](https://iopscience.iop.org/article/10.1088/1748-0221/17/11/P11003/pdf) paper by R. Abbasi et al 2022 JINST 17 P11003. It is one of the most recent works on neutrino detection predictions, which can be relevant for the current [IceCube - Neutrinos in Deep Ice](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice) competition.\n\nThe paper is some 30+ pages long and some of its parts are superfluous for the competiton in question. As a result, this notebook is a more focused summary of the paper and includes only the parts relevant for those who would like to compete.\n\nThe purpose of this notebook is to:\n\n* create an overview of the paper in question;\n* identify any useful insights;\n* understand the imperatives behind the use of GNNs;\n* understand how GNNs are implemented in the paper;\n* identify possible use of GNNs for the IceCube competition.","metadata":{}},{"cell_type":"markdown","source":"### Preparation of code for illustrations\n\nFirst of all, let's quickly read a few files from the competition that will be used for illustration purposes later on.","metadata":{}},{"cell_type":"code","source":"train_parquet = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_1.parquet')\nsensor_geometry = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv')\ntest_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ntrain_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train_meta.parquet')\nsubmission_sample = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/sample_submission.parquet')","metadata":{"execution":{"iopub.status.busy":"2023-03-01T19:55:19.545618Z","iopub.execute_input":"2023-03-01T19:55:19.546024Z","iopub.status.idle":"2023-03-01T19:56:07.952873Z","shell.execute_reply.started":"2023-03-01T19:55:19.545960Z","shell.execute_reply":"2023-03-01T19:56:07.951475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### IceCube detector engineering\n\nThe first part of the paper is a brief overview of the detector itself. While there is a bit of information on the topic in the competion description, paper adds quite some intriguing information. Please, have a closer look to this diagram.","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"The diagram [and the text afterwards] contains quite a lot of potentially usefull information. Main takeaways are:\n\n* Ice around detectors aren't homogeneous - their absorption properties differs, depending on the depth. \n* [As a result] Different detectors have different performance/precision\n* *DeepCore enhanced quantum efficiency detectors* (marked green on the picture) are **the most precise**, according to the authors, they're at least 1.35 times more precise than standard detectors \n* Other *detectors beneath the Dust layer* [layer of ice with dust impurities] perform well, but worse than DeepCore ones\n* *Detectors above the Dust layer* perform worse than those beneath it\n* *Detectors inside the Dust layer* are **the worst**, precision wise\n\nDatasets for the competiton contain positional data for the detectors. Let's check how they look like and if there are any discrepancies between the diagram and actual available data.","metadata":{}},{"cell_type":"code","source":"scatter3d = px.scatter_3d(sensor_geometry, x='x', y='y', z='z', opacity=0.7)\nscatter3d.update_traces(marker = dict(size = 2, symbol = \"diamond-open\"))\nscatter3d.update_coloraxes(showscale = False)\nscatter3d.update_layout(template = \"plotly\", font = dict(family = \"Arial\", size = 12, color = \"#9e97ff\"))\nscatter3d.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-01T19:56:07.954485Z","iopub.execute_input":"2023-03-01T19:56:07.954852Z","iopub.status.idle":"2023-03-01T19:56:10.364726Z","shell.execute_reply.started":"2023-03-01T19:56:07.954815Z","shell.execute_reply":"2023-03-01T19:56:10.363607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see on the scatter plot, detectors in the dataset are situated in the same way as those on the diagram - we can spot the DeepCore detectors, dust layer (empty space in DeepCore), etc. \n\n**Main conclusion regarding engineering**\n\nDetectors perform differently and it can be potentially useful to take this into account for better predictions.\n\nMore on detectors, their use and possible application of this information to the competition can be found in [this notebook](https://www.kaggle.com/code/antonsevostianov/icecube-sensor-efficiency-feature-engineering/notebook).","metadata":{}},{"cell_type":"markdown","source":"### The challenge of IceCube event classification and reconstruction and data input\n\nThe second chapter of the paper talks about problems that a researcher can bump into while trying to predict and classify IceCube events.\n\nBut the first question is what are we trying to predict and classify?\n\nPrediction is the main goal of this competition and the main target we want to predict is the direction of neutrino event. \n\nClassification is beyond the scope of this competition, but according to the paper, being able to classify and distinguish between neutrino and muon events has some scientific value (events detected are sometimes neutrino events and sometimes muon event, so classification is trying to guess which kind of event ws it).\n\nThe paper mentions several parameters of interest:\n\n* depositied energy\n* direction of the neutrino candidate\n* the interaction vertex \n\nLet's check training dataset for the Kaggle competition.","metadata":{}},{"cell_type":"code","source":"train_parquet.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-01T19:56:10.366800Z","iopub.execute_input":"2023-03-01T19:56:10.367230Z","iopub.status.idle":"2023-03-01T19:56:10.412228Z","shell.execute_reply.started":"2023-03-01T19:56:10.367191Z","shell.execute_reply":"2023-03-01T19:56:10.411074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see from our training dataset, only charge data is available in raw data. It is probably possible to reconstruct the direction of the event, based on sensor_id and time. But we don't have information regarding interaction vertex, at least in direct form.\n\nIt's interesting to compare it to the input data used by the authors in the paper. Below is the input data used by the authors.","metadata":{}},{"cell_type":"code","source":"features = ['Dxyz', 't', 'q', 'QE']\ndescription = ['Position of DOMs in IceCube coordinates', 'Pulse time relative to trigger time', 'Charge of a pulse', 'Quantum efficiency of PMT']\nunit = ['m', 'ns', 'P.E', '-']\ndataframe = {'Feature':features,'Description':description, 'Unit': unit}\nmy_df = pd.DataFrame(dataframe)\ndisplay(my_df)","metadata":{"execution":{"iopub.status.busy":"2023-03-01T20:13:02.552500Z","iopub.execute_input":"2023-03-01T20:13:02.553269Z","iopub.status.idle":"2023-03-01T20:13:02.575835Z","shell.execute_reply.started":"2023-03-01T20:13:02.553207Z","shell.execute_reply":"2023-03-01T20:13:02.574474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"DOM - is the abbreviation for detectors, PMT - tube that is a part of the detector. As we can see, the only data readily unavailable to us in Kaggle dataset, that the authors used in the paper is the quantum efficiency of the detectors (PMT).\n\n**Main conclusion regarding data input**\n\nOur kaggle datasets lack some of the data used by the researchers, but some of that data can be derived and added based on what is already available.","metadata":{}},{"cell_type":"markdown","source":"### Non-GNN Reconstruction Methods\n\nThe third chapter talks about methods that were historically used for neutrino events reconstruction and classification.\n\n*Maximum likelihood estimation* is one such technique. As authors point out this method, even in the most sophisticated form, can give **only a rough approximation**. However, there is an upside as well - those techniques are fast, computationally wise.\n\nAuthors mention *RETRO algorithm*, used by IceCube team, which is quite sophisticated, but **requires up to 40 seconds of CPU time** for proper reconstruction of just one event. Considering the number of the events, it's impossible to use this algorithm at the observatory (due to amount of computational infrustructure needed) or analyse real time data. \n\n*Convolutional Neural Networks* (CNNs) is another method used for the task. Paper points out that due to inherent architecture characteristics of CNNs, this approach lead to serious degradation of information for low-energy events. In other words, CNNs **don't work particularly well with low-energy events**.\n\n","metadata":{}},{"cell_type":"markdown","source":"### Graph Neural Networks\n\nNext, paper talks about GNNs and why they were chosen for the research. Let's see what GNN is and how it works \"under the hood\".","metadata":{}},{"cell_type":"markdown","source":"#### What is GNN\n\nGenerally, GNNs are NNs that work on **graph representations of data**. A graph consists of nodes\ninterconnected by edges. Nodes are associated with data, and the edges specify the relationship\namong the nodes. By adopting graphs as the input data structure, the idea of convolution generalizes\nfrom the application of filters on the rigid structure of grids to abstract mathematical operators\nthat utilize the interconnection of nodes in its computation.\n\nAuthors chose to represent **each event by a single graph**. Each observed pulse is represented by a\nnode in the graph and contains the per-DOM information shown in table 1. Each node in the graph is\nconnected to its **8 nearest neighbors** based on the Euclidean distance, and for this reason we consider the interconnectivity of the nodes in the graphs to be spatial.\n\nIn a nutshell, GNN creates a graph that clusters similar events (based on the Euclidean distance) and uses it to guess a target. Consider the picture below.\n","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"GNN clusters events based on certain features (kind of sport activity in this example, charge/time/direction/detector in our neutrino detection case) and uses this graphs to guess the target.","metadata":{}},{"cell_type":"markdown","source":"#### Event preprocessing and model architecture","metadata":{}},{"cell_type":"markdown","source":"Before using feeding data to GNN it is important to consider several facts:\n\n* While neural networks can in theory process data in any range of real numbers, **the complexity of the loss landscape** that a model navigates during training is highly dependent on the relative scale of the input data. \n* Distributions not centered around zero can lead to a **slower convergence time**. \n\nEach input variable π‘₯ is therefore transformed into π‘₯Β― using: \n\n\n\nwhere *𝑃𝑖thβ€žπ‘₯”* is the 𝑖-th percentile of the distribution of input feature π‘₯. This transformation brings\nthe input variables into roughly similar orders of magnitude, gives a median of zero and makes them\nunitless.\n","metadata":{}},{"cell_type":"markdown","source":"Some facts about the GNN used for IceCube:\n\n* It's called DYNEDGE. \n* It uses graph learning to extract features from pulses.\n* Implemented using GraphNeT.\n* Framework is built using PyTorch Geometric.\n* Uses a convolutional operator EdgeConv. \n\n\n\nThis picture shows the way model works. What happens on this diagram? \n\n1. First, the **following 5 global statistics are calculated** from the input graph: node homophily ratio of 𝐷xyz and 𝑑, and number of pulses in the graph. The homophily ratio is the ratio of connected node pairs that share the same node feature, and thus a number between 0 and 1. The homophily ratio of 𝐷xyz indicates what fraction of connected pulses originate from the same PMT.\n\n2. Second, the **input graph is propagated through 4 different EdgeConv blocks** such that the output from one flows into the next.\n\n3. Third, the **input graph and each state graph are concatenated together** into a [𝑛, 1030]-dimensional array that is passed through a fully connected Multilayer Perceptron (MLP) block containing two layers.\n\n4. This **array is aggregated node-wise into summary statistics** in four parallel ways; mean, min, max, and sum.\n\n5. **Aggregations are concatenated together** to minimize information loss, which produces an array with dimension [1, 4 Β· 256] that is concatenated together with the initially calculated 5 global statistics, producing a [1, 1029]-dimensional input array.\n\n6. This array is subsequently moved to **2-layer MLP block that makes the final prediction** by mapping the [1, 1029]-dimensional array to a [1, 𝑛outputs]-dimensional output.\n\nA DYNEDGE network is trained for each of the reconstruction variables: \n* deposited energy 𝐸\n* zenith\n* azimuth angles (πœƒ, πœ™)\n* the interaction vertex 𝑉xyz\n* 𝜈/πœ‡ classification\n* T /C classification.\n\nThis totals 6 independently trained models.","metadata":{}},{"cell_type":"markdown","source":"#### Main takeaways from GNN theory\n\n* GNNs are most probably better than CNNs for the task.\n* There are pretrained models and frameworks that can be used for the competition.\n* Hyperparameters of the model in question is optimised for low-energy events, considering that the competition isn't restricted to only such events, it might be a good idea to tweak the parameters","metadata":{}},{"cell_type":"markdown","source":"### Conclusions\n\nIn this notebook we have skimmed through a paper regarding IceCube neutrino detection that can be relevant to the Kaggle competiton, learned about IceCube engineering, data used by the scientists and GNNs.\n\nWe have found out that: \n\n* Detectors perform differently and it can be potentially useful to take this into account for better predictions.\n* Our Kaggle datasets lack some of the data used by the researchers, but some of that data can be derived and added based on what is already available.\n* GNNs are most probably better than CNNs for the task.\n* There are pretrained models and frameworks that can be used for the competition.\n* Hyperparameters of the model in question is optimised for low-energy events, considering that the competition isn't restricted to only such events, it might be a good idea to tweak the parameters","metadata":{}}]}