{"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❄️: Fitting line to 3D point cloud\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### See related notebooks:\n\n- [IceCube🧊: Neutrino🎆EDA & 3D🔭interactive viewer](https://www.kaggle.com/code/jirkaborovec/icecube-neutrino-eda-3d-interactive-viewer)\n- [IceCube❄️: baseline with XGBoost regression](https://www.kaggle.com/code/jirkaborovec/icecube-neutrino-baseline-xgboost-regression)\n\n### Point-cloud 3D fitting\n\n- https://scikit-spatial.readthedocs.io/en/stable/gallery/fitting/plot_line_3d.html","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-27T08:44:09.268352Z","iopub.execute_input":"2023-01-27T08:44:09.268854Z","iopub.status.idle":"2023-01-27T08:44:09.276916Z","shell.execute_reply.started":"2023-01-27T08:44:09.268819Z","shell.execute_reply":"2023-01-27T08:44:09.275905Z"},"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_train = pd.read_parquet(os.path.join(PATH_DATASET, \"train_meta.parquet\"))\ndisplay(meta_train.head())","metadata":{"execution":{"iopub.status.busy":"2023-01-27T08:45:04.950564Z","iopub.execute_input":"2023-01-27T08:45:04.951060Z","iopub.status.idle":"2023-01-27T08:45:48.405709Z","shell.execute_reply.started":"2023-01-27T08:45:04.951025Z","shell.execute_reply":"2023-01-27T08:45:48.402687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(meta_train.info())","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-27T08:46:44.270058Z","iopub.execute_input":"2023-01-27T08:46:44.270637Z","iopub.status.idle":"2023-01-27T08:46:44.295437Z","shell.execute_reply.started":"2023-01-27T08:46:44.270586Z","shell.execute_reply":"2023-01-27T08:46:44.294356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### azimuth and zenith","metadata":{}},{"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-27T08:46:56.286417Z","iopub.execute_input":"2023-01-27T08:46:56.286939Z","iopub.status.idle":"2023-01-27T08:47:14.508806Z","shell.execute_reply.started":"2023-01-27T08:46:56.286901Z","shell.execute_reply":"2023-01-27T08:47:14.507874Z"},"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-27T08:47:52.104749Z","iopub.execute_input":"2023-01-27T08:47:52.107511Z","iopub.status.idle":"2023-01-27T08:47:52.155453Z","shell.execute_reply.started":"2023-01-27T08:47:52.107425Z","shell.execute_reply":"2023-01-27T08:47:52.153921Z"},"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-27T08:55:00.357122Z","iopub.execute_input":"2023-01-27T08:55:00.357630Z","iopub.status.idle":"2023-01-27T08:55:00.440029Z","shell.execute_reply.started":"2023-01-27T08:55:00.357579Z","shell.execute_reply":"2023-01-27T08:55:00.438694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"geometry.set_index(\"sensor_id\", inplace=True)\ngeometry = geometry.apply(np.float32)\ngeometry.info()","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-27T09:05:24.823811Z","iopub.execute_input":"2023-01-27T09:05:24.826136Z","iopub.status.idle":"2023-01-27T09:05:24.855769Z","shell.execute_reply.started":"2023-01-27T09:05:24.826061Z","shell.execute_reply":"2023-01-27T09:05:24.853880Z"},"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- **sensor_id** (int): the ID of which of the 5160 IceCube photomultiplier sensors recorded this pulse.\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":"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-27T09:05:31.543681Z","iopub.execute_input":"2023-01-27T09:05:31.544516Z","iopub.status.idle":"2023-01-27T09:05:37.396987Z","shell.execute_reply.started":"2023-01-27T09:05:31.544475Z","shell.execute_reply":"2023-01-27T09:05:37.395671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fitting line to points","metadata":{}},{"cell_type":"code","source":"! pip install -q scikit-spatial --no-index -f /kaggle/input/icecube-neutrino-eda-3d-interactive-viewer/frozen-packages/","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-27T09:07:18.109162Z","iopub.execute_input":"2023-01-27T09:07:18.109690Z","iopub.status.idle":"2023-01-27T09:07:30.528346Z","shell.execute_reply.started":"2023-01-27T09:07:18.109647Z","shell.execute_reply":"2023-01-27T09:07:30.526786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\n\ndef cartesian_to_polar(x, y, z):\n    # https://en.wikipedia.org/wiki/Spherical_coordinate_system\n    x2y2 = x**2 + y**2\n    r = math.sqrt(x2y2 + z**2)\n    azimuth = math.acos(x / math.sqrt(x2y2)) * np.sign(y)\n    zenith = math.acos(z / r)\n    return azimuth, zenith\n\ndef adjust_polar(azimuth, zenith):\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    elif zenith < 0:\n        zenith += math.pi\n    return azimuth, zenith\n\ndef polar_to_cartesian(azimuth, zenith):\n    # see: https://stackoverflow.com/a/10868220/4521646\n    x = math.cos(azimuth) * math.sin(zenith)\n    y = math.sin(azimuth) * math.sin(zenith)\n    z = math.cos(zenith)\n    return x, y, z","metadata":{"execution":{"iopub.status.busy":"2023-01-27T09:08:36.448727Z","iopub.execute_input":"2023-01-27T09:08:36.449257Z","iopub.status.idle":"2023-01-27T09:08:36.463669Z","shell.execute_reply.started":"2023-01-27T09:08:36.449217Z","shell.execute_reply":"2023-01-27T09:08:36.462248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.graph_objects as go\n\ndef draw_subplot(\n    fig, i, evt, sensors, meta, pos_offset, pos_direction, gt_direction, line_len=250\n):\n    # sensors as background\n    fig.add_trace(go.Scatter3d(\n        x=sensors['x'], y=sensors['y'], z=sensors['z'], \n        mode='markers', marker=dict(size=1, color=0), opacity=0.2\n    ), row=(i+1), col=1)\n    # sensors reading\n    fig.add_trace(go.Scatter3d(\n        x=evt['x'], y=evt['y'], z=evt['z'], opacity=0.8,\n        mode='markers', marker=dict(\n            size=evt['charge'] * 10, color=evt['time'], colorscale='temps'\n        )\n    ), row=(i+1), col=1)\n    # direction from metad data\n    ox, oy, oz = pos_offset\n    \n    for c, (x, y, z) in [('red', pos_direction), ('green', gt_direction)]:\n        fig.add_trace(go.Scatter3d(\n            x=[ox - x * line_len, ox + x * line_len],\n            y=[oy - y * line_len, oy + y * line_len],\n            z=[oz - z * line_len, oz + z * line_len],\n            opacity=0.8, mode='lines', line=dict(color=c, width=3)\n        ), row=(i+1), col=1)","metadata":{"execution":{"iopub.status.busy":"2023-01-27T09:10:26.911940Z","iopub.execute_input":"2023-01-27T09:10:26.912489Z","iopub.status.idle":"2023-01-27T09:10:26.929341Z","shell.execute_reply.started":"2023-01-27T09:10:26.912453Z","shell.execute_reply":"2023-01-27T09:10:26.927833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample event","metadata":{}},{"cell_type":"code","source":"from plotly.subplots import make_subplots\nfrom skspatial.objects import Line, Points\n\n\ndef show_event(event_id=46528394, data=train, sensors=geometry, metadata=meta_train):\n    event = data[data.index == event_id]\n    event = event.merge(sensors, left_on=\"sensor_id\", right_index=True)\n    meta = dict(metadata[metadata[\"event_id\"]==event_id].iloc[0])\n    display(meta)\n    \n    evt_ = event[~event['auxiliary']]\n    points = Points(evt_[[\"x\", \"y\", \"z\"]].values)\n    line = Line.best_fit(points)\n    print(line)\n    azimuth_, zenith_ = adjust_polar(*cartesian_to_polar(*line.direction))\n    print(f\"predictions of azimuth={azimuth_} & zenith={zenith_}\")\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    x_, y_, z_ = polar_to_cartesian(meta['azimuth'], meta['zenith'])\n    for i, aux in enumerate(auxiliaries):\n        draw_subplot(\n            fig, i, event[event['auxiliary'] == aux], sensors, meta,\n            line.point, line.direction, gt_direction=(x_, y_, z_),\n        )\n    desc_gt = f\"azimuth={meta['azimuth']:0.3}; zenith={meta['zenith']:0.3}\"\n    desc_gt += f\"\\n -> x={x_:0.2}; y={y_:0.2}; z={z_:0.2}\"\n    fig.update_layout(\n        height=800, width=600, showlegend=False,\n        title_text=f\"Event #{event_id} / {desc_gt}\",\n    )\n    return fig\n\nshow_event().show()","metadata":{"execution":{"iopub.status.busy":"2023-01-27T09:11:30.201053Z","iopub.execute_input":"2023-01-27T09:11:30.201560Z","iopub.status.idle":"2023-01-27T09:11:31.516918Z","shell.execute_reply.started":"2023-01-27T09:11:30.201521Z","shell.execute_reply":"2023-01-27T09:11:31.515718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Interactive - browse events\n\n**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-27T09:12:11.824149Z","iopub.execute_input":"2023-01-27T09:12:11.824793Z","iopub.status.idle":"2023-01-27T09:12:12.595350Z","shell.execute_reply.started":"2023-01-27T09:12:11.824744Z","shell.execute_reply":"2023-01-27T09:12:12.594433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc, time\n\ndel train, meta_train, events\ngc.collect(), time.sleep(9)","metadata":{"execution":{"iopub.status.busy":"2023-01-27T09:20:31.637902Z","iopub.execute_input":"2023-01-27T09:20:31.638433Z","iopub.status.idle":"2023-01-27T09:20:40.816836Z","shell.execute_reply.started":"2023-01-27T09:20:31.638397Z","shell.execute_reply":"2023-01-27T09:20:40.815455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare 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\"))\nssub.set_index(\"event_id\", inplace=True)\ndisplay(ssub.head())","metadata":{"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-01-27T09:20:53.798717Z","iopub.execute_input":"2023-01-27T09:20:53.799152Z","iopub.status.idle":"2023-01-27T09:20:53.828699Z","shell.execute_reply.started":"2023-01-27T09:20:53.799112Z","shell.execute_reply":"2023-01-27T09:20:53.827210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ssub = ssub.apply(np.float32)\ndisplay(ssub.info())","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-27T09:21:06.056458Z","iopub.execute_input":"2023-01-27T09:21:06.057734Z","iopub.status.idle":"2023-01-27T09:21:06.076011Z","shell.execute_reply.started":"2023-01-27T09:21:06.057674Z","shell.execute_reply":"2023-01-27T09:21:06.074836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fitting to all data","metadata":{}},{"cell_type":"code","source":"ls = glob.glob(os.path.join(PATH_DATASET, \"test\", \"*.parquet\"))\n\nfor batch_file in ls:\n    print(f\"processing: {batch_file}\")\n    df = pd.read_parquet(batch_file)\n    # del df['time'], df['charge']\n    # display(df.head())\n    gc.collect()\n    \n    for eid, dfg in df.groupby(\"event_id\"):\n        dfg = dfg[~dfg['auxiliary']]\n        dfg = dfg.merge(geometry, left_on=\"sensor_id\", right_index=True)\n        # TODO: for memory reason subsample point cloud\n        if len(dfg) > 1200:\n            dfg = dfg.sort_values('charge', ascending=False)[:1200]\n        # display(dfg)\n\n        try:\n            points = Points(dfg[[\"x\", \"y\", \"z\"]].values)\n            line = Line.best_fit(points)\n            direction_ = -1 * line.direction\n            azimuth_, zenith_ = adjust_polar(*cartesian_to_polar(*direction_))\n        except Exception:\n            azimuth_, zenith_ = 0., 0.\n        else:\n            ssub.at[eid, \"azimuth\"] = azimuth_\n            ssub.at[eid, \"zenith\"] = zenith_\n        if len(df) < 1e5:\n            print(f\"Estimation {line} with azimuth={azimuth_} & zenith={zenith_}\")\n\n    del df, dfg\n    gc.collect(), time.sleep(9)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-27T09:21:15.190442Z","iopub.execute_input":"2023-01-27T09:21:15.190992Z","iopub.status.idle":"2023-01-27T09:21:24.617334Z","shell.execute_reply.started":"2023-01-27T09:21:15.190949Z","shell.execute_reply":"2023-01-27T09:21:24.615740Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Finalize submission","metadata":{}},{"cell_type":"code","source":"ssub.fillna(0).to_csv('submission.csv', index=True)\n\n!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-01-27T09:21:49.807554Z","iopub.execute_input":"2023-01-27T09:21:49.808135Z","iopub.status.idle":"2023-01-27T09:21:51.125204Z","shell.execute_reply.started":"2023-01-27T09:21:49.808092Z","shell.execute_reply":"2023-01-27T09:21:51.123403Z"},"trusted":true},"execution_count":null,"outputs":[]}]}