{"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":"markdown","source":"# Summary","metadata":{}},{"cell_type":"markdown","source":"In this notebook I work with the [IceCube - Neutrinos in Deep Ice](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/data) Dataset provided by the IceCube Neutrino Observatory. \n\nI define some functions in order to create a dataframe that allows me to perform Singular Value Decomposition (SVD) on the x,y,z coordinates of the pulses related to the event ID, and sensor location from the batch parquet data file. SVD is used to find the best fit line in the 3D space. I take into consideration the auxiliary=True to perform a simple basic prediction of the neutrino trajectory (comparing it to the real trajectory located in the meta parquet data file). \n\nMost of the logic behind the functions comes from [GREYSNOW](https://www.kaggle.com/code/shlomoron/icecube-eda-pca-baseline-cv-1-28-lb-1-274) notebook. Go check him out! \n\nFor more detailed information on the neutrino observatory please refer to the following [website](https://icecube.wisc.edu/science/icecube/).\n\nThe following image contains a sketch of how the IceCube observatory is structured:","metadata":{}},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"## What does the neutrino observatory do? ","metadata":{}},{"cell_type":"markdown","source":"The IceCube Observatory can detect neutrino particles in the energy range of a few tens of GeV to a few PeV, allowing scientist to map cosmic events that emit neutrinos between the energy of those electron volt ranges. \n\nNeutrinos are very hard to detect, they have almost negligible mass, and they rarely interact with anything. They have a slight chance to collide with the ice core where the sensors (DOM's) are located. The interaction between the neutrino and the ice produces a muon which later travels along the sensor geometry structure and emits Cherenkov light which gets detected by the DOM's. We call the light readings from the DOM's pulses.\n\nTo better visualize about what we talked about let's see this simplified image about the whole mechanics:\n\n[source of image](https://www.sciencedirect.com/science/article/abs/pii/S0146641018300346)","metadata":{}},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"Now that we superficial understanding of how it all works let's jump right into predicting the neutrino trajectory for a single event. ","metadata":{}},{"cell_type":"markdown","source":"# Library import","metadata":{}},{"cell_type":"markdown","source":"As usual we import all the necessary libraries used throughout the notebook.","metadata":{}},{"cell_type":"code","source":"# Magic function that will make your plot outputs appear and be stored within the notebook\n%matplotlib inline\n\n# Function used to to render higher resolution images\n%config InlineBackend.figure_format = 'retina'\n\n# Ignore all warnings\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n# Data manipulation\nimport os\nimport pandas as pd\nimport numpy as np\n\npd.set_option('display.max_columns', None)\npd.set_option('display.expand_frame_repr', False)\n\n# Data visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly\nplotly.offline.init_notebook_mode()\nimport plotly.express as px\nimport plotly.graph_objects as go\n\n# Standardizing the style for the visualizations \nsns.set_theme()\nsns.set(font_scale=1.2)\nsns.set_palette(\"pastel\")\nplt.style.use('seaborn-whitegrid')\n\n\n# Singular Vector Decomposition\nfrom sklearn.decomposition import TruncatedSVD","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:18.742215Z","iopub.execute_input":"2023-02-05T06:13:18.742596Z","iopub.status.idle":"2023-02-05T06:13:20.762531Z","shell.execute_reply.started":"2023-02-05T06:13:18.742523Z","shell.execute_reply":"2023-02-05T06:13:20.761581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset import and description","metadata":{}},{"cell_type":"markdown","source":"In this seccion we describe each of the features available. We can see that the data parquet files are divided in to three large dataset (meta data, batch data and geometry data). \n\nFor a more detailed explanation of each feature please take a look at the [data card](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/data).\n\nLet's start by defining the principal path where the packages are located.","metadata":{}},{"cell_type":"code","source":"# Principal path \npath= \"/kaggle/input/icecube-neutrinos-in-deep-ice\"","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:20.764142Z","iopub.execute_input":"2023-02-05T06:13:20.764840Z","iopub.status.idle":"2023-02-05T06:13:20.769082Z","shell.execute_reply.started":"2023-02-05T06:13:20.764816Z","shell.execute_reply":"2023-02-05T06:13:20.767881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since we are only doing a basic prediction (not training any model) we will only use the training set for this notebook.","metadata":{}},{"cell_type":"markdown","source":"## Meta training data","metadata":{}},{"cell_type":"code","source":"# Meta train dataset path\ntrain_meta_p=\"train_meta.parquet\"","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:20.770718Z","iopub.execute_input":"2023-02-05T06:13:20.770990Z","iopub.status.idle":"2023-02-05T06:13:20.788105Z","shell.execute_reply.started":"2023-02-05T06:13:20.770967Z","shell.execute_reply":"2023-02-05T06:13:20.787051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The meta data contains the following features:","metadata":{}},{"cell_type":"markdown","source":"**[train/test]_meta.parquet**\n\n* **batch_id** (int): the ID of the batch the event was placed into.\n\n* **event_id** (int): the event ID.\n\n* **[first/last]_pulse_index** (int): index of the first/last row in the features dataframe belonging to this event.\n\n* **[azimuth/zenith]** (float32): the [azimuth/zenith] angle in radians of the neutrino. A value between 0 and 2*pi for the azimuth and 0 and pi for zenith. The target columns. Not provided for the test set. The direction vector represented by zenith and azimuth points to where the neutrino came from.","metadata":{}},{"cell_type":"markdown","source":"We proceed to read the training meta parquet data file.","metadata":{}},{"cell_type":"code","source":"# Reading the meta parquet file \ntrain_meta=pd.read_parquet(os.path.join(path,train_meta_p))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:20.789484Z","iopub.execute_input":"2023-02-05T06:13:20.789787Z","iopub.status.idle":"2023-02-05T06:13:55.935052Z","shell.execute_reply.started":"2023-02-05T06:13:20.789754Z","shell.execute_reply":"2023-02-05T06:13:55.933657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a sneak peek at the train meta data.","metadata":{}},{"cell_type":"code","source":"# Train meta dataset sneak peek\ntrain_meta.head(3).append(train_meta.tail(3))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:55.937455Z","iopub.execute_input":"2023-02-05T06:13:55.937882Z","iopub.status.idle":"2023-02-05T06:13:55.960583Z","shell.execute_reply.started":"2023-02-05T06:13:55.937858Z","shell.execute_reply":"2023-02-05T06:13:55.958969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It appears we have in total 131,954,923.00 instances. Let's confirm this by measuring the dataframe length.","metadata":{}},{"cell_type":"code","source":"len(train_meta)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:55.961943Z","iopub.execute_input":"2023-02-05T06:13:55.962331Z","iopub.status.idle":"2023-02-05T06:13:55.967642Z","shell.execute_reply.started":"2023-02-05T06:13:55.962297Z","shell.execute_reply":"2023-02-05T06:13:55.967051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Batch data","metadata":{}},{"cell_type":"code","source":"# Train batch data path \ntrain_batch_p= \"train/batch_31.parquet\"","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:55.968618Z","iopub.execute_input":"2023-02-05T06:13:55.968970Z","iopub.status.idle":"2023-02-05T06:13:55.979368Z","shell.execute_reply.started":"2023-02-05T06:13:55.968947Z","shell.execute_reply":"2023-02-05T06:13:55.978541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's take a look at what features are inside the batch_[n] parquet data file:","metadata":{}},{"cell_type":"markdown","source":"**[train/test]/batch_[n].parquet** - Each batch contains tens of thousands of events. Each event may contain thousands of pulses, each of which is the digitized output from a photomultiplier tube and occupies one row.\n\n* **event_id** (int): the event ID. Saved as the index column in parquet.\n\n* **time** (int): the time of the pulse in nanoseconds in the current event time window. \nThe absolute time of a pulse has no relevance, and only the relative time with respect to other pulses within an event is of relevance.\n\n* **sensor_id** (int): the ID of which of the 5160 IceCube photomultiplier sensors recorded this pulse.\n\n* **charge** (float32): An estimate of the amount of light in the pulse, in units of photoelectrons (p.e.). A physical photon does not exactly result in a measurement of 1 p.e. but rather can take values spread around 1 p.e. As an example, a pulse with charge 2.7 p.e. could quite likely be the result of two or three photons hitting the photomultiplier tube around the same time. This data has float16 precision but is stored as float32 due to limitations of the version of pyarrow the data was prepared with.\n\n* **auxiliary** (bool): If True, the pulse was not fully digitized, is of lower quality, and was more likely to originate from noise. If False, then this pulse was contributed to the trigger decision and the pulse was fully digitized.","metadata":{}},{"cell_type":"markdown","source":"According to the meta data we have 660 batch files to choose from. Let's choose number a random number like perhabs batch #31.","metadata":{}},{"cell_type":"markdown","source":"Let's now read the batch file.","metadata":{}},{"cell_type":"code","source":"# Reading batch file\ntrain_batch=pd.read_parquet(os.path.join(path,train_batch_p))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:13:55.980723Z","iopub.execute_input":"2023-02-05T06:13:55.980986Z","iopub.status.idle":"2023-02-05T06:14:00.445748Z","shell.execute_reply.started":"2023-02-05T06:13:55.980958Z","shell.execute_reply":"2023-02-05T06:14:00.443296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's inspect a little further.","metadata":{}},{"cell_type":"code","source":"# Showing first and last instances of batch #30\ntrain_batch.head(3).append(train_batch.tail(3))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.447710Z","iopub.execute_input":"2023-02-05T06:14:00.448257Z","iopub.status.idle":"2023-02-05T06:14:00.464217Z","shell.execute_reply.started":"2023-02-05T06:14:00.448220Z","shell.execute_reply":"2023-02-05T06:14:00.462826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's now see how many instances does this batch contains.","metadata":{}},{"cell_type":"code","source":"len(train_batch)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.465940Z","iopub.execute_input":"2023-02-05T06:14:00.466298Z","iopub.status.idle":"2023-02-05T06:14:00.481527Z","shell.execute_reply.started":"2023-02-05T06:14:00.466267Z","shell.execute_reply":"2023-02-05T06:14:00.480269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Sensor geometry / location data","metadata":{}},{"cell_type":"code","source":"# Sensor(DOM) geometry data path\ngeo_p=\"sensor_geometry.csv\"","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.482748Z","iopub.execute_input":"2023-02-05T06:14:00.484067Z","iopub.status.idle":"2023-02-05T06:14:00.492495Z","shell.execute_reply.started":"2023-02-05T06:14:00.484000Z","shell.execute_reply":"2023-02-05T06:14:00.491728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's import the sensor geometry csv file.","metadata":{}},{"cell_type":"markdown","source":"**sensor_geometry.csv** - The x, y, and z positions for each of the 5160 IceCube sensors. The row index corresponds to the sensor_idx feature of pulses. The x, y, and z coordinates are in units of meters, with the origin at the center of the IceCube detector. The coordinate system is right-handed, and the z-axis points upwards when standing at the South Pole. You can convert from these coordinates to azimuth and zenith with the following formulas (here the vector (x,y,z) is normalized):\n\nx = cos(azimuth) * sin(zenith)\ny = sin(azimuth) * sin(zenith)\nz = cos(zenith)\n","metadata":{}},{"cell_type":"markdown","source":"We start by reading the geometry data file.","metadata":{}},{"cell_type":"code","source":"# Reading geometry data\ngeo=pd.read_csv(os.path.join(path,geo_p))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.493603Z","iopub.execute_input":"2023-02-05T06:14:00.494740Z","iopub.status.idle":"2023-02-05T06:14:00.520643Z","shell.execute_reply.started":"2023-02-05T06:14:00.494707Z","shell.execute_reply":"2023-02-05T06:14:00.519234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Before doing a projection of the 3D space reflecting the DOM's position let's have a look at the first and last instances of the dataset.","metadata":{}},{"cell_type":"code","source":"# Showing first and last instances of Sensor geometry\ngeo.head(3).append(geo.tail(3))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.524523Z","iopub.execute_input":"2023-02-05T06:14:00.524899Z","iopub.status.idle":"2023-02-05T06:14:00.540339Z","shell.execute_reply.started":"2023-02-05T06:14:00.524862Z","shell.execute_reply":"2023-02-05T06:14:00.538898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I create a dataframe containing the limit of the projection space.","metadata":{}},{"cell_type":"code","source":"# Projection space limits\nindex=[\"x-min\",\"x-max\",\"y-min\",\"y-max\",\"z-min\",\"z-max\"]\n\ndata= {\n 'Limits': [geo.x.min(),geo.x.max(),\n geo.y.min(),geo.y.max(),\n geo.z.min(),geo.z.max()]\n}\n\ngeo_limits=pd.DataFrame(data=data,index=index)\n\ngeo_limits.T","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.545775Z","iopub.execute_input":"2023-02-05T06:14:00.546406Z","iopub.status.idle":"2023-02-05T06:14:00.564346Z","shell.execute_reply.started":"2023-02-05T06:14:00.546372Z","shell.execute_reply":"2023-02-05T06:14:00.563558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following plot shows a top view layout of the sensors location:","metadata":{}},{"cell_type":"code","source":"# 2D view of the Sensor distribution\n\nfig, ax =plt.subplots(figsize=(8,8))\nsns.despine()\nax.tick_params(axis='x', labelrotation=0)\nax.tick_params(axis='y', labelrotation=0)\nax.set(xlabel='x', ylabel='y')\nax.set_title('Sensor Geometry -2D', size=20)\nsns.scatterplot(x='x',y='y',hue=\"sensor_id\",data=geo,marker='s',s=100,palette='plasma',legend=False)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:00.565435Z","iopub.execute_input":"2023-02-05T06:14:00.565753Z","iopub.status.idle":"2023-02-05T06:14:01.183660Z","shell.execute_reply.started":"2023-02-05T06:14:00.565692Z","shell.execute_reply":"2023-02-05T06:14:01.180435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"3D projection of the sensors layout.","metadata":{}},{"cell_type":"code","source":"# 3D view of the Sensor distribution\n\nfig = px.scatter_3d(geo, x='x', y='y', z='z', opacity=0.6, color=\"sensor_id\",title=\"Sensor Geometry -3D\")\nfig.update_traces(marker_size=2)\nfig.update_layout(height=800, width=800)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:01.187219Z","iopub.execute_input":"2023-02-05T06:14:01.188222Z","iopub.status.idle":"2023-02-05T06:14:03.192963Z","shell.execute_reply.started":"2023-02-05T06:14:01.188101Z","shell.execute_reply":"2023-02-05T06:14:03.192195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Function definition","metadata":{}},{"cell_type":"markdown","source":"Here we define the principal functions to be used for data exploration and the basic prediction.\n\nFirst off we start by defining the metric in order to evaluate the quality of our prediction. In this case we will be using the angular error. I referenced this specific function from [SOHIER DANE](https://www.kaggle.com/code/sohier/mean-angular-error) notebook.","metadata":{}},{"cell_type":"code","source":"#Mean Angular Error - Evaluation metric\n\ndef angular_dist_score(az_true, zen_true, az_pred, zen_pred):\n '''\n calculate the MAE of the angular distance between two directions.\n The two vectors are first converted to cartesian unit vectors,\n and then their scalar product is computed, which is equal to\n the cosine of the angle between the two vectors. The inverse \n cosine (arccos) thereof is then the angle between the two input vectors\n \n Parameters:\n -----------\n \n az_true : float (or array thereof)\n true azimuth value(s) in radian\n zen_true : float (or array thereof)\n true zenith value(s) in radian\n az_pred : float (or array thereof)\n predicted azimuth value(s) in radian\n zen_pred : float (or array thereof)\n predicted zenith value(s) in radian\n \n Returns:\n --------\n \n dist : float\n mean over the angular distance(s) in radian\n '''\n \n if not (np.all(np.isfinite(az_true)) and\n np.all(np.isfinite(zen_true)) and\n np.all(np.isfinite(az_pred)) and\n np.all(np.isfinite(zen_pred))):\n raise ValueError(\"All arguments must be finite\")\n \n # pre-compute all sine and cosine values\n sa1 = np.sin(az_true)\n ca1 = np.cos(az_true)\n sz1 = np.sin(zen_true)\n cz1 = np.cos(zen_true)\n \n sa2 = np.sin(az_pred)\n ca2 = np.cos(az_pred)\n sz2 = np.sin(zen_pred)\n cz2 = np.cos(zen_pred)\n \n # scalar product of the two cartesian vectors (x = sz*ca, y = sz*sa, z = cz)\n scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n \n # scalar product of two unit vectors is always between -1 and 1, this is against nummerical instability\n # that might otherwise occure from the finite precision of the sine and cosine functions\n scalar_prod = np.clip(scalar_prod, -1, 1)\n \n # convert back to an angle (in radian)\n return np.average(np.abs(np.arccos(scalar_prod)))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.194258Z","iopub.execute_input":"2023-02-05T06:14:03.194639Z","iopub.status.idle":"2023-02-05T06:14:03.203456Z","shell.execute_reply.started":"2023-02-05T06:14:03.194615Z","shell.execute_reply":"2023-02-05T06:14:03.202212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Defining a function to get pulses coordinates depending on the event number. ","metadata":{}},{"cell_type":"code","source":"def get_pulse_coor(event_num):\n \n # get event index in the meta dataset\n event_idx=event_idx=int(\"\".join(str(i) for i in(train_meta[train_meta.event_id==event_num].index.tolist())))\n \n # slice training batch dataset to include only the instances from the selected event\n event_pulses = train_batch.iloc[train_meta.iloc[event_idx].first_pulse_index.astype(int) : \n train_meta.iloc[event_idx].last_pulse_index.astype(int)+1].copy()\n \n # reset back indexes to dataframe \n event_pulses = event_pulses.reset_index()\n \n # map sensor location (x, y, z) to corresponding sensor id \n event_pulses[['x', 'y', 'z']] = geo.loc[event_pulses.sensor_id].loc[:, ['x', 'y', 'z']].values\n\n \n #centering the data points for SVD\n event_pulses[['x', 'y', 'z']] -= event_pulses[['x', 'y', 'z']].mean()\n\n return event_pulses","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.205593Z","iopub.execute_input":"2023-02-05T06:14:03.205939Z","iopub.status.idle":"2023-02-05T06:14:03.225578Z","shell.execute_reply.started":"2023-02-05T06:14:03.205914Z","shell.execute_reply":"2023-02-05T06:14:03.224416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Function to get target vector coordinates.","metadata":{}},{"cell_type":"code","source":"def target_vector(event_num):\n \n # get event index in the meta dataset\n event_idx=event_idx=int(\"\".join(str(i) for i in(train_meta[train_meta.event_id==event_num].index.tolist())))\n \n # get neutrino true trajectory from the meta dataset\n zenith_target = train_meta.loc[event_idx, \"zenith\"]\n azimuth_target = train_meta.loc[event_idx, \"azimuth\"]\n \n # transform azimuth and zenith in to coordinates\n vector_target = [np.cos(azimuth_target) * np.sin(zenith_target), \n np.sin(azimuth_target) * np.sin(zenith_target),\n np.cos(zenith_target)]\n \n # multiply target vector by base vector (maginutude 500 and -500) for visualization purpose\n vector_base = np.array([-500, 500])\n x = vector_base*vector_target[0]\n y = vector_base*vector_target[1]\n z = vector_base*vector_target[2]\n \n # Covert the target vector to dataframe \n vector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n \n return vector_target_df","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.227119Z","iopub.execute_input":"2023-02-05T06:14:03.227582Z","iopub.status.idle":"2023-02-05T06:14:03.244643Z","shell.execute_reply.started":"2023-02-05T06:14:03.227557Z","shell.execute_reply":"2023-02-05T06:14:03.243566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Function to fit pulses coordinates to SVD and get the predicted vector dataframe. This function will be mainly used to visualize the fitted vector in the 3D scatter plot. ","metadata":{}},{"cell_type":"code","source":"def fit_vector(event_pulses):\n \n # Instantiate SVD and fit with respective \n svd = TruncatedSVD(n_components=1).fit(event_pulses[[\"x\",\"y\",\"z\"]])\n \n # Decompose SVD components \n vector= svd.components_[0]\n \n # multiply the fit vector by base vector (maginutude 500 and -500) for visualization purpose\n vector_base = np.array([-500, 500])\n x = vector_base*vector[0]\n y = vector_base*vector[1]\n z = vector_base*vector[2]\n \n # Covert the fit vector to dataframe \n vector_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n \n return vector_df","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.245788Z","iopub.execute_input":"2023-02-05T06:14:03.246677Z","iopub.status.idle":"2023-02-05T06:14:03.257985Z","shell.execute_reply.started":"2023-02-05T06:14:03.246652Z","shell.execute_reply":"2023-02-05T06:14:03.256601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The next function allows us to get the best fitted line(vector in this case) regarding the x,y,z coordinates and then convert it to Azimuth and Zenith so we can calculate the error between the target and fitted vector.","metadata":{}},{"cell_type":"code","source":"def covert_az_ze(event_pulses):\n \n # Instantiate SVD and fit with respective \n svd = TruncatedSVD(n_components=1).fit(event_pulses[[\"x\",\"y\",\"z\"]])\n \n # Decompose SVD components \n vector= svd.components_[0]\n \n # Transform to zenith/azimuth (in radians)\n zenith = np.arccos(vector[2])\n azimuth = np.arctan2(vector[1], vector[0])\n \n if azimuth<0:\n azimuth = 2*np.pi + azimuth\n return zenith,azimuth\n \n else:\n return zenith,azimuth","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.259173Z","iopub.execute_input":"2023-02-05T06:14:03.260082Z","iopub.status.idle":"2023-02-05T06:14:03.276555Z","shell.execute_reply.started":"2023-02-05T06:14:03.260047Z","shell.execute_reply":"2023-02-05T06:14:03.275368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Function used to plot neutrinos pulses and trajectory in the 3D space.","metadata":{}},{"cell_type":"code","source":"def plot_event(event_pulses,vector,target_vector,auxiliary):\n \n # 3D scaterplot of auxiliary pulses \n fig_auxiliary = px.scatter_3d(event_pulses.loc[event_pulses.auxiliary],\n x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'],size='charge')\n \n # 3D scaterplot of non-auxiliary pulses \n fig_non_auxiliary = px.scatter_3d(event_pulses.loc[~event_pulses.auxiliary],\n x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'],size='charge')\n \n # Target neutrino trajectory\n fig_line_target = px.line_3d(target_vector, x=\"x\", y=\"y\", z=\"z\")\n \n # Predicted neutrino trajectory (best fit auxiliary = True)\n fig_line_svd= px.line_3d(vector, x=\"x\", y=\"y\", z=\"z\",color_discrete_sequence=['magenta'])\n \n # Plot auxiliary pulses\n if auxiliary == True:\n \n fig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line_svd.data + fig_line_target.data)\n \n \n return fig \n \n # Remove auxiliary pulses \n else:\n \n fig = go.Figure(data = fig_non_auxiliary.data + fig_line_svd.data + fig_line_target.data)\n \n return fig \n ","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.277740Z","iopub.execute_input":"2023-02-05T06:14:03.278979Z","iopub.status.idle":"2023-02-05T06:14:03.290318Z","shell.execute_reply.started":"2023-02-05T06:14:03.278928Z","shell.execute_reply":"2023-02-05T06:14:03.289067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploratory Analysis","metadata":{}},{"cell_type":"markdown","source":"Before doing a prediction let's explore the data a bit more. But before let's make sure what azimuth and zenith represent. The following image does a good job explaining both of them:","metadata":{}},{"cell_type":"markdown","source":"
","metadata":{}},{"cell_type":"markdown","source":"From the meta data file let's see the target azimuth and zenith distribution. We want to know what’s the most prominent direction of the neutrino in the dataset. ","metadata":{}},{"cell_type":"code","source":"# Plot Azimuth and Zenith \n\nfig, ax =plt.subplots(2,2,gridspec_kw={'height_ratios': [3, 1]},figsize=(16,8))\n\nsns.despine(right=True, top=True, left=True)\n\nax[0,0].tick_params(axis='x', labelrotation=0)\nax[0,0].tick_params(axis='y', labelrotation=0)\nax[0,0].set_title('Azimuth - Distribution', size=20)\nax[0,0].set(xlabel='Azimuth', ylabel='Count')\n\nax[0,1].tick_params(axis='x', labelrotation=0)\nax[0,1].tick_params(axis='y', labelrotation=0)\nax[0,1].set_title('Zenith - Distribution', size=20)\nax[0,1].set(xlabel='Zenith', ylabel='Count')\n\nsns.histplot(data=train_meta, x='azimuth',y=None,bins=50,ax=ax[0,0])\nsns.histplot(data=train_meta, x='zenith',y=None,bins=50,ax=ax[0,1])\n\nsns.boxenplot(data=train_meta, x='azimuth',y=None,ax=ax[1,0])\nsns.boxenplot(data=train_meta, x='zenith',y=None,ax=ax[1,1])\n","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:14:03.292301Z","iopub.execute_input":"2023-02-05T06:14:03.293829Z","iopub.status.idle":"2023-02-05T06:17:51.889719Z","shell.execute_reply.started":"2023-02-05T06:14:03.293796Z","shell.execute_reply":"2023-02-05T06:17:51.886800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The azimuth ranges from 0° to 360° and it seems consistent throughout the whole meta data. This is normally what we expected.\n\nNow the Zenith which ranges from 0° to 180° tells us that we are most likely to find a neutrino coming at an angle between 60° and 120° (90° being the most common). It appears that horizontal(trajectories with zenith around 0 or 3 radians) are quite rare.\n\nLastly for a better visualization let's compare both targets.\n\n","metadata":{}},{"cell_type":"code","source":"# Overlaping the two graphs Blue is Azimuth Red is Zenith\n\nfig, ax =plt.subplots(figsize=(8,8))\n\nsns.despine(right=True, top=True, left=True)\n\nax.set_title('Azimuth vs Zenith', size=20)\nax.set(xlabel='Azimuth / Zenith', ylabel='Count')\n\nsns.histplot(data=train_meta,x='azimuth',color='r',bins=50,ax=ax)\nsns.histplot(data=train_meta,x='zenith',color='b',bins=50,ax=ax)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:17:51.892810Z","iopub.execute_input":"2023-02-05T06:17:51.894540Z","iopub.status.idle":"2023-02-05T06:19:45.235451Z","shell.execute_reply.started":"2023-02-05T06:17:51.894492Z","shell.execute_reply":"2023-02-05T06:19:45.234359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Basic prediction","metadata":{}},{"cell_type":"markdown","source":"We are ready to make the basic prediction using SVD to get the best fitted line regarding the event pulses coordinates x, y and z.\n\nLet's explore a random event located in the training batch #31.\n\nThe comments in each code cell will be explaining the process for the prediction step by step.","metadata":{}},{"cell_type":"code","source":"# Calculating target vector for event #100867570\ntv =target_vector(100867570)\ntv","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:45.236939Z","iopub.execute_input":"2023-02-05T06:19:45.237292Z","iopub.status.idle":"2023-02-05T06:19:45.620329Z","shell.execute_reply.started":"2023-02-05T06:19:45.237257Z","shell.execute_reply":"2023-02-05T06:19:45.619105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Obtaining pulse coordinates \nfit_coor=get_pulse_coor(100867570)\nfit_coor","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:45.621925Z","iopub.execute_input":"2023-02-05T06:19:45.622259Z","iopub.status.idle":"2023-02-05T06:19:45.819318Z","shell.execute_reply.started":"2023-02-05T06:19:45.622227Z","shell.execute_reply":"2023-02-05T06:19:45.817811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating best fit line\nv_aux=fit_vector(fit_coor)\nv_aux","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:45.821261Z","iopub.execute_input":"2023-02-05T06:19:45.821602Z","iopub.status.idle":"2023-02-05T06:19:45.857940Z","shell.execute_reply.started":"2023-02-05T06:19:45.821579Z","shell.execute_reply":"2023-02-05T06:19:45.857221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting event # 100867570 auxiliary = True \n#Magenta line is predicted trajectory / Blue line is true trajectory\n\nplot_event(fit_coor,v_aux,tv,True)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:45.859101Z","iopub.execute_input":"2023-02-05T06:19:45.859926Z","iopub.status.idle":"2023-02-05T06:19:46.093114Z","shell.execute_reply.started":"2023-02-05T06:19:45.859892Z","shell.execute_reply":"2023-02-05T06:19:46.092071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pulse coordinates auxiliary \nno_aux=fit_coor.loc[~fit_coor.auxiliary]\nno_aux","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.094897Z","iopub.execute_input":"2023-02-05T06:19:46.095864Z","iopub.status.idle":"2023-02-05T06:19:46.118718Z","shell.execute_reply.started":"2023-02-05T06:19:46.095811Z","shell.execute_reply":"2023-02-05T06:19:46.117655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating best fit line\nv_non_aux=fit_vector(no_aux)\nv_non_aux","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.120474Z","iopub.execute_input":"2023-02-05T06:19:46.121696Z","iopub.status.idle":"2023-02-05T06:19:46.145229Z","shell.execute_reply.started":"2023-02-05T06:19:46.121659Z","shell.execute_reply":"2023-02-05T06:19:46.144091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting event # 100867570 auxiliary = False\n#Magenta line is predicted trajectory / Blue line is true trajectory\n\nplot_event(fit_coor,v_non_aux,tv,False)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.146653Z","iopub.execute_input":"2023-02-05T06:19:46.147969Z","iopub.status.idle":"2023-02-05T06:19:46.354600Z","shell.execute_reply.started":"2023-02-05T06:19:46.147936Z","shell.execute_reply":"2023-02-05T06:19:46.353226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Calculating error","metadata":{}},{"cell_type":"code","source":"def get_event_idx(event_num):\n event_idx=event_idx=int(\"\".join(str(i) for i in(train_meta[train_meta.event_id==event_num].index.tolist())))\n return event_idx","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.355838Z","iopub.execute_input":"2023-02-05T06:19:46.356099Z","iopub.status.idle":"2023-02-05T06:19:46.361693Z","shell.execute_reply.started":"2023-02-05T06:19:46.356068Z","shell.execute_reply":"2023-02-05T06:19:46.360829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_event_idx(100867570)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.362828Z","iopub.execute_input":"2023-02-05T06:19:46.363102Z","iopub.status.idle":"2023-02-05T06:19:46.549192Z","shell.execute_reply.started":"2023-02-05T06:19:46.363047Z","shell.execute_reply":"2023-02-05T06:19:46.548053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zt = train_meta.iloc[6199999].zenith\nat = train_meta.iloc[6199999].azimuth","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.550494Z","iopub.execute_input":"2023-02-05T06:19:46.550772Z","iopub.status.idle":"2023-02-05T06:19:46.556850Z","shell.execute_reply.started":"2023-02-05T06:19:46.550746Z","shell.execute_reply":"2023-02-05T06:19:46.555690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Error with auxiliary pulses","metadata":{}},{"cell_type":"code","source":"z,a= covert_az_ze(fit_coor)\n\nprint('Zenith Prediction=',round(z,4))\nprint(\"Target Zenith=\",round(zt,4), \"\\n\")\n\nprint('Azimuth Prediction=',round(a,4))\nprint(\"Target Azimuth=\",round(at,4),\"\\n\")\n\nerror=angular_dist_score(at, zt, a, z)\nprint('Error w/ auxiliary pulses=',round(error,4))\n","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.558203Z","iopub.execute_input":"2023-02-05T06:19:46.558438Z","iopub.status.idle":"2023-02-05T06:19:46.574275Z","shell.execute_reply.started":"2023-02-05T06:19:46.558416Z","shell.execute_reply":"2023-02-05T06:19:46.573429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Error without auxiliary pulses","metadata":{}},{"cell_type":"code","source":"z,a= covert_az_ze(no_aux)\n\nprint('Zenith Prediction=',round(z,4))\nprint(\"Target Zenith=\",round(zt,4), \"\\n\")\n\nprint('Azimuth Prediction=',round(a,4))\nprint(\"Target Azimuth=\",round(at,4),\"\\n\")\n\nerror=angular_dist_score(at, zt, a, z)\nprint('Error w/o auxiliary pulses=',round(error,4))","metadata":{"execution":{"iopub.status.busy":"2023-02-05T06:19:46.575843Z","iopub.execute_input":"2023-02-05T06:19:46.576147Z","iopub.status.idle":"2023-02-05T06:19:46.594686Z","shell.execute_reply.started":"2023-02-05T06:19:46.576117Z","shell.execute_reply":"2023-02-05T06:19:46.593792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Afterword","metadata":{}},{"cell_type":"markdown","source":"Even if we drop the auxiliary pulses for this particular event the best fitted 3D line had quite noticeable discrepancy from the target trajectory. We could try changing the type of model to a more complex one plus input some constraints on the data we feed to the model.\n\nObviously to evaluate the true potential of this method we will have to train it on all the batches and calculate the respective error with the meta data.\n\nFeel free to fork this notebook and try with other events and see how SVD performs. \n","metadata":{}}]}