{"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":"# IceCube❄️: EDA & 3D🔎 interactive viewer\n\nThe goal of this competition is to predict a neutrino particle’s direction. You will develop a model based on data from the \"IceCube\" detector, which observes the cosmos from deep within the South Pole ice.\n\nThe IceCube Neutrino Observatory is the first detector of its kind, encompassing a cubic kilometer of ice and designed to search for the nearly massless neutrinos. An international group of scientists is responsible for the scientific research that makes up the IceCube Collaboration.\n\nBy making the process faster and more precise, you'll help improve the reconstruction of neutrinos. As a result, we could gain a clearer image of our universe.\n\n![illustration](https://storage.googleapis.com/kaggle-media/competitions/IceCube/icecube_detector.jpg)","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nPATH_DATASET = \"/kaggle/input/icecube-neutrinos-in-deep-ice\"","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-28T14:33:36.822768Z","iopub.execute_input":"2023-01-28T14:33:36.823877Z","iopub.status.idle":"2023-01-28T14:33:36.854833Z","shell.execute_reply.started":"2023-01-28T14:33:36.823775Z","shell.execute_reply":"2023-01-28T14:33:36.853525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse meta data\n\n**[train/test]_meta.parquet**\n\n- **batch_id** (int): the ID of the batch the event was placed into.\n- **event_id** (int): the event ID.\n- **[first/last]_pulse_index** (int): index of the first/last row in the features dataframe belonging to this event.\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":"code","source":"meta_test = pd.read_parquet(os.path.join(PATH_DATASET, \"test_meta.parquet\"))\nprint(f\"length: {len(meta_test)}\")\nmeta_test.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:33:36.857230Z","iopub.execute_input":"2023-01-28T14:33:36.858310Z","iopub.status.idle":"2023-01-28T14:33:37.008345Z","shell.execute_reply.started":"2023-01-28T14:33:36.858262Z","shell.execute_reply":"2023-01-28T14:33:37.007168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_train = pd.read_parquet(os.path.join(PATH_DATASET, \"train_meta.parquet\"))\nprint(f\"length: {len(meta_train)}\")\ndisplay(meta_train.head())\ndisplay(meta_train.info())","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:33:37.010117Z","iopub.execute_input":"2023-01-28T14:33:37.010851Z","iopub.status.idle":"2023-01-28T14:34:17.062827Z","shell.execute_reply.started":"2023-01-28T14:33:37.010806Z","shell.execute_reply":"2023-01-28T14:34:17.061834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(6, 2))\nmeta_train.groupby(\"batch_id\").size().plot.hist(bins=50)\nplt.xlabel(\"nb. events in batch\"), plt.ylabel(\"nb. cases\")\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:34:17.064882Z","iopub.execute_input":"2023-01-28T14:34:17.065320Z","iopub.status.idle":"2023-01-28T14:34:20.289738Z","shell.execute_reply.started":"2023-01-28T14:34:17.065292Z","shell.execute_reply":"2023-01-28T14:34:20.288653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"seems that event id is unique per whole dataset as grouping is 1","metadata":{}},{"cell_type":"code","source":"meta_train.groupby(\"event_id\").size().max()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:34:20.291193Z","iopub.execute_input":"2023-01-28T14:34:20.291538Z","iopub.status.idle":"2023-01-28T14:35:34.167630Z","shell.execute_reply.started":"2023-01-28T14:34:20.291506Z","shell.execute_reply":"2023-01-28T14:35:34.166421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Hitogram of first and last pulses","metadata":{}},{"cell_type":"code","source":"meta_train[[\"first_pulse_index\", \"last_pulse_index\"]].plot.hist(bins=100, alpha=0.5, figsize=(8, 3))\nplt.xlabel('time')","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:35:34.169218Z","iopub.execute_input":"2023-01-28T14:35:34.170298Z","iopub.status.idle":"2023-01-28T14:35:52.682865Z","shell.execute_reply.started":"2023-01-28T14:35:34.170262Z","shell.execute_reply":"2023-01-28T14:35:52.681677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Delay between first and last pulse","metadata":{}},{"cell_type":"code","source":"meta_train[\"delay_pulse_index\"] = meta_train[\"last_pulse_index\"] - meta_train[\"first_pulse_index\"]\n_ = plt.hist(meta_train[\"delay_pulse_index\"], bins=100, log=True)\nplt.ylabel('count cases'), plt.xlabel('duration')\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:35:52.684167Z","iopub.execute_input":"2023-01-28T14:35:52.684486Z","iopub.status.idle":"2023-01-28T14:35:56.552038Z","shell.execute_reply.started":"2023-01-28T14:35:52.684440Z","shell.execute_reply":"2023-01-28T14:35:56.550872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### azimuth and zenith\n\nThe solar azimuth and solar zenith express the position of the sun. The solar azimuth is the angle of the direction of the sun measured clockwise north from the horizon. The solar zenith is the angle measured from the local zenith and the line of sight of the sun.\n\n![illustration](https://www.researchgate.net/publication/354755033/figure/fig2/AS:1070891857362948@1632331923300/Schematic-depicting-the-solar-zenith-angle-solar-altitude-angle-and-solar-azimuth-angle_W640.jpg)","metadata":{}},{"cell_type":"code","source":"meta_train[[\"azimuth\", \"zenith\"]].plot.hist(bins=50, alpha=0.5, figsize=(8, 3))\nplt.xlabel('azimuth/zenith')\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:35:56.553549Z","iopub.execute_input":"2023-01-28T14:35:56.554001Z","iopub.status.idle":"2023-01-28T14:36:42.210530Z","shell.execute_reply.started":"2023-01-28T14:35:56.553972Z","shell.execute_reply":"2023-01-28T14:36:42.209129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist2d(meta_train[\"azimuth\"], meta_train[\"zenith\"], bins=(50, 50), cmap=plt.cm.jet)\nplt.xlabel('azimuth'), plt.ylabel('zenith')\nplt.colorbar()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:36:42.212941Z","iopub.execute_input":"2023-01-28T14:36:42.213282Z","iopub.status.idle":"2023-01-28T14:36:58.837443Z","shell.execute_reply.started":"2023-01-28T14:36:42.213250Z","shell.execute_reply":"2023-01-28T14:36:58.836419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse geometry\n\nThe 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.","metadata":{}},{"cell_type":"code","source":"geometry = pd.read_csv(os.path.join(PATH_DATASET, \"sensor_geometry.csv\"))\nprint(f\"length: {len(geometry)}\")\ngeometry.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:36:58.841214Z","iopub.execute_input":"2023-01-28T14:36:58.841573Z","iopub.status.idle":"2023-01-28T14:36:58.867154Z","shell.execute_reply.started":"2023-01-28T14:36:58.841542Z","shell.execute_reply":"2023-01-28T14:36:58.866165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\n\nfig = px.scatter_3d(geometry, x='x', y='y', z='z', opacity=0.6, color=\"sensor_id\")\nfig.update_traces(marker_size=2)\nfig.update_layout(height=600, width=600)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:36:58.868282Z","iopub.execute_input":"2023-01-28T14:36:58.868610Z","iopub.status.idle":"2023-01-28T14:37:01.166413Z","shell.execute_reply.started":"2023-01-28T14:36:58.868583Z","shell.execute_reply":"2023-01-28T14:37:01.165304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse training/test data\n\n**[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- **time** (int): the time of the pulse in nanoseconds in the current event time window. The 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- **sensor_id** (int): the ID of which of the 5160 IceCube photomultiplier sensors recorded this pulse.\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- **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":"code","source":"test = pd.read_parquet(os.path.join(PATH_DATASET, \"test/batch_661.parquet\"))\nprint(f\"length: {len(test)}\")\nprint(f\"events: {len(test.index.unique())}\")\ntest.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:01.168142Z","iopub.execute_input":"2023-01-28T14:37:01.169165Z","iopub.status.idle":"2023-01-28T14:37:01.200993Z","shell.execute_reply.started":"2023-01-28T14:37:01.169120Z","shell.execute_reply":"2023-01-28T14:37:01.199843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_parquet(os.path.join(PATH_DATASET, \"train/batch_15.parquet\"))\nprint(f\"length: {len(train)}\")\nprint(f\"events: {len(train.index.unique())}\")\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:01.202286Z","iopub.execute_input":"2023-01-28T14:37:01.202612Z","iopub.status.idle":"2023-01-28T14:37:05.399740Z","shell.execute_reply.started":"2023-01-28T14:37:01.202583Z","shell.execute_reply":"2023-01-28T14:37:05.398722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"seems that not all event have same number of samples","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10, 4))\n_ = plt.hist(train.groupby(level=0).size(), bins=50, log=True)\nplt.xlabel(\"nb. ALL measument per event\"), plt.ylabel(\"nb. of event\")\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:05.401115Z","iopub.execute_input":"2023-01-28T14:37:05.401472Z","iopub.status.idle":"2023-01-28T14:37:06.830387Z","shell.execute_reply.started":"2023-01-28T14:37:05.401423Z","shell.execute_reply":"2023-01-28T14:37:06.829369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 4))\n_ = plt.hist(train[~train['auxiliary']].groupby(level=0).size(), bins=50, log=True)\nplt.xlabel(\"nb. (aux==False) measument per event\"), plt.ylabel(\"nb. of event\")\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:06.832196Z","iopub.execute_input":"2023-01-28T14:37:06.832618Z","iopub.status.idle":"2023-01-28T14:37:08.895232Z","shell.execute_reply.started":"2023-01-28T14:37:06.832578Z","shell.execute_reply":"2023-01-28T14:37:08.894154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Chages per auxiliary","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\n_ = plt.hist(train[~train['auxiliary']]['charge'], bins=50, log=True, label=\"False\")\n_ = plt.hist(train[train['auxiliary']]['charge'], bins=50, log=True, label=\"True\")\nplt.ylabel('count cases'), plt.xlabel('charge')\nplt.grid(), plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:08.896825Z","iopub.execute_input":"2023-01-28T14:37:08.897145Z","iopub.status.idle":"2023-01-28T14:37:11.878706Z","shell.execute_reply.started":"2023-01-28T14:37:08.897116Z","shell.execute_reply":"2023-01-28T14:37:11.877705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample event\n\ninspiration from: https://github.com/plotly/plotly.py/blob/master/doc/python/subplots.md","metadata":{}},{"cell_type":"code","source":"! pip download -q --prefer-binary scikit-spatial --dest frozen-packages/\n! ls -lh frozen-packages/","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-28T14:37:11.879865Z","iopub.execute_input":"2023-01-28T14:37:11.880170Z","iopub.status.idle":"2023-01-28T14:37:20.475277Z","shell.execute_reply.started":"2023-01-28T14:37:11.880142Z","shell.execute_reply":"2023-01-28T14:37:20.473641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_train[meta_train['event_id'].isin([46528394, 2135637939, 2084362251])]","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:20.479779Z","iopub.execute_input":"2023-01-28T14:37:20.480488Z","iopub.status.idle":"2023-01-28T14:37:21.545047Z","shell.execute_reply.started":"2023-01-28T14:37:20.480423Z","shell.execute_reply":"2023-01-28T14:37:21.544009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\n\ndef charge_center(event):\n    evt_sub = event[~event['auxiliary']][['x', 'y', 'z', 'charge']]\n    evt_sub.loc[:, 'coef'] = evt_sub['charge'] / evt_sub['charge'].sum()\n    # display(evt_sub)\n    for c in ['x', 'y', 'z']:\n        evt_sub.loc[:, c] *= evt_sub['coef']\n    cx, cy, cz = evt_sub[['x', 'y', 'z']].sum().values\n    return cx, cy, cz\n\ndef direction(meta):\n    azimuth = meta[\"azimuth\"]\n    zenith = meta[\"zenith\"]\n    dx = math.sin(zenith) * math.cos(azimuth)\n    dy = math.sin(zenith) * math.sin(azimuth)\n    dz = math.cos(zenith)\n    return dx, dy, dz","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:21.546340Z","iopub.execute_input":"2023-01-28T14:37:21.546743Z","iopub.status.idle":"2023-01-28T14:37:21.555072Z","shell.execute_reply.started":"2023-01-28T14:37:21.546714Z","shell.execute_reply":"2023-01-28T14:37:21.554100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.graph_objects as go\n\ndef draw_subplot(fig, i, evt, sensors, cx, cy, cz, dx, dy, dz, scale=400):\n    # sensors as background\n    fig.add_trace(\n        go.Scatter3d(\n            x=sensors['x'], y=sensors['y'], z=sensors['z'], \n            mode='markers', marker=dict(size=1, color=\"black\"), opacity=0.2\n        ), row=(i+1), col=1)\n    # sensors reading\n    fig.add_trace(\n        go.Scatter3d(\n            x=evt['x'], y=evt['y'], z=evt['z'], opacity=0.8,\n            mode='markers', marker=dict(\n                size=evt['charge'] * 15,\n                color=evt['time'],\n                colorscale='sunsetdark',\n            )\n        ), row=(i+1), col=1)\n    # direction from metad data\n    fig.add_trace(\n        go.Scatter3d(\n            x=[cx - dx * scale, cx + dx * scale],\n            y=[cy - dy * scale, cy + dy * scale],\n            z=[cz - dz * scale, cz + dz * scale],\n            opacity=0.8, mode='lines', line=dict(color='red', width=3)\n        ), row=(i+1), col=1)","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:40:32.301616Z","iopub.execute_input":"2023-01-28T14:40:32.302012Z","iopub.status.idle":"2023-01-28T14:40:32.313024Z","shell.execute_reply.started":"2023-01-28T14:40:32.301980Z","shell.execute_reply":"2023-01-28T14:40:32.311891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from plotly.subplots import make_subplots\n\n\ndef show_event(event_id=46528394, data=train, sensors=geometry, metadata=meta_train):\n    meta = dict(metadata[metadata[\"event_id\"]==event_id].iloc[0])\n    display(meta)\n    event = data[data.index == event_id]\n    event = event.merge(sensors, on=\"sensor_id\")\n    event.loc[:, 'charge'] /= event['charge'].max()\n    dx, dy, dz = direction(meta)\n    cx, cy, cz = charge_center(event)\n    #print(cx, cy, cz)\n    \n    auxiliaries = [False, True]\n    fig = make_subplots(\n        rows=2, specs=[[{'type': 'scene'}], [{'type': 'scene'}]],\n        subplot_titles=[f\"auxiliary={aux}\" for aux in auxiliaries],\n        vertical_spacing=0.05,\n    )\n    for i, aux in enumerate(auxiliaries):\n        evt_ = event[event['auxiliary'] == aux]\n        draw_subplot(fig, i, evt_, sensors, cx, cy, cz, dx, dy, dz)\n    fig.update_layout(\n        height=800, width=600, showlegend=False,\n        title_text=f\"Event #{event_id}\"\n        f\" / azimuth={meta['azimuth']:0.3}; zenith={meta['zenith']:0.3}\\n\"\n        f\"-> x={dx:0.2}; y={dy:0.2}; z={dz:0.2}\",\n    )\n    # fig = px.scatter_3d(event, x='x', y='y', z='z', size=\"charge\", color=\"auxiliary\",\n    #     opacity=0.8, title=f\"Event #{event_id}\\n -> x={x_}; y={y_}; z={z_}\")\n    return fig\n\nshow_event().show()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:40:50.469994Z","iopub.execute_input":"2023-01-28T14:40:50.470435Z","iopub.status.idle":"2023-01-28T14:40:50.810767Z","shell.execute_reply.started":"2023-01-28T14:40:50.470398Z","shell.execute_reply":"2023-01-28T14:40:50.809599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**This one is an interactive chart when you can slide over all event in the training batch**\n\nbut to get in you need to ioen the notebook in edior mode (so to copy as your own)","metadata":{}},{"cell_type":"code","source":"from ipywidgets import interact, IntSlider\n\ndef interactive_show(events):\n    interact(\n        lambda i: show_event(events[i]).show(),\n        i=IntSlider(min=0, max=len(events), step=1, value=len(events) // 2),\n    )\n\nevents = train.index.unique().tolist()\ninteractive_show(events)","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:40:58.066318Z","iopub.execute_input":"2023-01-28T14:40:58.067373Z","iopub.status.idle":"2023-01-28T14:40:58.724582Z","shell.execute_reply.started":"2023-01-28T14:40:58.067333Z","shell.execute_reply":"2023-01-28T14:40:58.723398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse submission\n\nAn example submission with the correct columns and properly ordered event IDs. The sample submission is provided in the parquet format so it can be read quickly but your final submission must be a csv.","metadata":{}},{"cell_type":"code","source":"ssub = pd.read_parquet(os.path.join(PATH_DATASET, \"sample_submission.parquet\"))\nprint(f\"length: {len(ssub)}\")\nssub.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:22.203062Z","iopub.status.idle":"2023-01-28T14:37:22.203605Z","shell.execute_reply.started":"2023-01-28T14:37:22.203313Z","shell.execute_reply":"2023-01-28T14:37:22.203337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ssub['zenith'] = meta_train['zenith'].mean()\nssub['azimuth'] = 3.21\nssub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:22.204716Z","iopub.status.idle":"2023-01-28T14:37:22.205313Z","shell.execute_reply.started":"2023-01-28T14:37:22.204980Z","shell.execute_reply":"2023-01-28T14:37:22.205007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-01-28T14:37:22.207241Z","iopub.status.idle":"2023-01-28T14:37:22.207824Z","shell.execute_reply.started":"2023-01-28T14:37:22.207547Z","shell.execute_reply":"2023-01-28T14:37:22.207573Z"},"trusted":true},"execution_count":null,"outputs":[]}]}