{"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❄️: 3D point cloud fit with RANSAC\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🧊: Neutrino🎆 fitting 3D🌪️points cloud](https://www.kaggle.com/code/jirkaborovec/icecube-neutrino-fitting-3d-points-lb-1-25)\n\n### Point-cloud 3D fitting\n\n- https://leomariga.github.io/pyRANSAC-3D/api-documentation/line\n- https://scikit-image.org/docs/stable/auto_examples/transform/plot_ransac.html\n- http://man.hubwiz.com/docset/Scikit-image.docset/Contents/Resources/Documents/auto_examples/transform/plot_ransac3D.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-02-25T00:46:44.047587Z","iopub.execute_input":"2023-02-25T00:46:44.048058Z","iopub.status.idle":"2023-02-25T00:46:44.080105Z","shell.execute_reply.started":"2023-02-25T00:46:44.047973Z","shell.execute_reply":"2023-02-25T00:46:44.079261Z"},"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-02-25T00:46:44.081602Z","iopub.execute_input":"2023-02-25T00:46:44.082047Z","iopub.status.idle":"2023-02-25T00:47:22.007184Z","shell.execute_reply.started":"2023-02-25T00:46:44.082019Z","shell.execute_reply":"2023-02-25T00:47:22.005925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(meta_train.info())","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-02-25T00:47:22.009163Z","iopub.execute_input":"2023-02-25T00:47:22.009655Z","iopub.status.idle":"2023-02-25T00:47:22.028134Z","shell.execute_reply.started":"2023-02-25T00:47:22.009624Z","shell.execute_reply":"2023-02-25T00:47:22.026785Z"},"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-02-25T00:47:22.031231Z","iopub.execute_input":"2023-02-25T00:47:22.031558Z","iopub.status.idle":"2023-02-25T00:47:38.333618Z","shell.execute_reply.started":"2023-02-25T00:47:22.031529Z","shell.execute_reply":"2023-02-25T00:47:38.331692Z"},"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-02-25T00:47:38.336787Z","iopub.execute_input":"2023-02-25T00:47:38.337318Z","iopub.status.idle":"2023-02-25T00:47:38.372623Z","shell.execute_reply.started":"2023-02-25T00:47:38.337262Z","shell.execute_reply":"2023-02-25T00:47:38.370810Z"},"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-02-25T00:47:38.374779Z","iopub.execute_input":"2023-02-25T00:47:38.375734Z","iopub.status.idle":"2023-02-25T00:47:41.209787Z","shell.execute_reply.started":"2023-02-25T00:47:38.375675Z","shell.execute_reply":"2023-02-25T00:47:41.208233Z"},"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-02-25T00:47:41.212030Z","iopub.execute_input":"2023-02-25T00:47:41.212538Z","iopub.status.idle":"2023-02-25T00:47:41.234403Z","shell.execute_reply.started":"2023-02-25T00:47:41.212497Z","shell.execute_reply":"2023-02-25T00:47:41.233465Z"},"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-02-25T00:47:41.237426Z","iopub.execute_input":"2023-02-25T00:47:41.238587Z","iopub.status.idle":"2023-02-25T00:47:45.826947Z","shell.execute_reply.started":"2023-02-25T00:47:41.238542Z","shell.execute_reply":"2023-02-25T00:47:45.825225Z"},"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,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-02-25T00:47:45.829033Z","iopub.execute_input":"2023-02-25T00:47:45.829541Z","iopub.status.idle":"2023-02-25T00:47:57.209237Z","shell.execute_reply.started":"2023-02-25T00:47:45.829495Z","shell.execute_reply":"2023-02-25T00:47:57.207325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"see the angle visualisation in https://www.kaggle.com/code/jirkaborovec/icecube-visualize-azimuth-zenith\n\n<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/a/a2/Kugelkoord-lokale-Basis-s.svg/1024px-Kugelkoord-lokale-Basis-s.svg.png\" width=\"400\" height=\"400\"/>","metadata":{}},{"cell_type":"code","source":"import math\nimport warnings\nwarnings.filterwarnings('ignore')\n\ndef cartesian_to_sphere(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_sphere(azimuth, zenith):\n    if zenith < 0:\n        zenith += math.pi\n        azimuth += math.pi\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    azimuth = azimuth % (2 * math.pi)\n    return azimuth, zenith\n\ndef sphere_to_cartesian(azimuth, zenith):\n    # see: https://stackoverflow.com/a/10868220/4521646\n    x = math.sin(zenith) * math.cos(azimuth)\n    y = math.sin(zenith) * math.sin(azimuth)\n    z = math.cos(zenith)\n    return x, y, z","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-02-25T00:47:57.213975Z","iopub.execute_input":"2023-02-25T00:47:57.215129Z","iopub.status.idle":"2023-02-25T00:47:57.226484Z","shell.execute_reply.started":"2023-02-25T00:47:57.215079Z","shell.execute_reply":"2023-02-25T00:47:57.224196Z"},"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=\"black\"), opacity=0.2\n    ), row=(i+1), col=1)\n    # sensors reading\n    evt_time_min, evt_time_max = min(evt['time']), max(evt['time'])\n    evt.loc[:, 'time'] = (evt['time'] - evt_time_min) / (evt_time_max - evt_time_min)\n    fig.add_trace(go.Scatter3d(\n        x=evt['x'], y=evt['y'], z=evt['z'], opacity=0.95,\n        mode='markers', marker=dict(\n            size=evt['charge'] * 15, color=evt['time'], colorscale='sunsetdark'\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)\n        fig.add_trace(go.Scatter3d(\n            x=[ox - x * line_len],\n            y=[oy - y * line_len],\n            z=[oz - z * line_len],\n            opacity=0.8, mode='markers', marker=dict(size=5, color=c)\n        ), row=(i+1), col=1)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-02-25T00:47:57.228961Z","iopub.execute_input":"2023-02-25T00:47:57.229524Z","iopub.status.idle":"2023-02-25T00:47:57.246978Z","shell.execute_reply.started":"2023-02-25T00:47:57.229481Z","shell.execute_reply":"2023-02-25T00:47:57.245048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample event","metadata":{}},{"cell_type":"code","source":"from plotly.subplots import make_subplots\nfrom skimage.measure import LineModelND, ransac\nfrom skspatial.objects import Line\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    print(meta)\n    \n    event.loc[:, 'charge'] = event['charge'] / max(event['charge'])\n    event = event.sort_values('time')\n    evt_ = event[~event['auxiliary']]\n    xyz = evt_[[\"x\", \"y\", \"z\"]].values\n    # robustly fit line only using inlier data with RANSAC algorithm\n    model_robust, inliers = ransac(\n        xyz, LineModelND, min_samples=2, residual_threshold=1, max_trials=1000,\n    )\n    print(f\"Estimate: {model_robust.params}\")\n    origin, direction = model_robust.params\n    line = Line(point=origin, direction=direction)\n    \n    print(f\"provided/GT of azimuth={meta['azimuth']:0.5}; zenith={meta['zenith']:0.5}\")\n    azimuth_, zenith_ = cartesian_to_sphere(*line.direction)\n    print(f\"predictions of azimuth={azimuth_:0.5} & zenith={zenith_:0.5}\")\n    azimuth_, zenith_ = adjust_sphere(azimuth_, zenith_)\n    print(f\"correction of azimuth={azimuth_:0.5} & zenith={zenith_:0.5}\")\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_ = sphere_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-02-25T00:51:35.885208Z","iopub.execute_input":"2023-02-25T00:51:35.885712Z","iopub.status.idle":"2023-02-25T00:51:36.246992Z","shell.execute_reply.started":"2023-02-25T00:51:35.885672Z","shell.execute_reply":"2023-02-25T00:51:36.245234Z"},"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-02-25T00:52:49.930022Z","iopub.execute_input":"2023-02-25T00:52:49.930497Z","iopub.status.idle":"2023-02-25T00:52:50.557915Z","shell.execute_reply.started":"2023-02-25T00:52:49.930458Z","shell.execute_reply":"2023-02-25T00:52:50.555897Z"},"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-02-25T00:54:35.802983Z","iopub.execute_input":"2023-02-25T00:54:35.803425Z","iopub.status.idle":"2023-02-25T00:54:44.973758Z","shell.execute_reply.started":"2023-02-25T00:54:35.803392Z","shell.execute_reply":"2023-02-25T00:54:44.972535Z"},"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-02-25T00:54:44.976124Z","iopub.execute_input":"2023-02-25T00:54:44.976520Z","iopub.status.idle":"2023-02-25T00:54:45.001489Z","shell.execute_reply.started":"2023-02-25T00:54:44.976484Z","shell.execute_reply":"2023-02-25T00:54:45.000223Z"},"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-02-25T00:54:45.002771Z","iopub.execute_input":"2023-02-25T00:54:45.003136Z","iopub.status.idle":"2023-02-25T00:54:45.021561Z","shell.execute_reply.started":"2023-02-25T00:54:45.003106Z","shell.execute_reply":"2023-02-25T00:54:45.019035Z"},"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) > 800:\n            dfg = dfg.sort_values('charge', ascending=False)[:800]\n        # display(dfg)\n\n        try:\n            xyz = dfg[[\"x\", \"y\", \"z\"]].values\n            model_robust, inliers = ransac(\n                xyz, LineModelND, min_samples=2, residual_threshold=1, max_trials=1000,\n            )\n            origin, direction = model_robust.params\n            line = Line(point=origin, direction=direction)\n            # TOOD: adjust by timeseriers\n            direction_ = -1 * line.direction\n            azimuth_, zenith_ = adjust_sphere(*cartesian_to_sphere(*direction_))\n        except Exception as ex:\n            azimuth_, zenith_ = 0., 0.\n            print(ex)\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-02-25T00:54:45.024584Z","iopub.execute_input":"2023-02-25T00:54:45.025299Z","iopub.status.idle":"2023-02-25T00:54:54.560077Z","shell.execute_reply.started":"2023-02-25T00:54:45.025251Z","shell.execute_reply":"2023-02-25T00:54:54.558562Z"},"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-02-25T00:54:54.561631Z","iopub.execute_input":"2023-02-25T00:54:54.562034Z","iopub.status.idle":"2023-02-25T00:54:54.970234Z","shell.execute_reply.started":"2023-02-25T00:54:54.562001Z","shell.execute_reply":"2023-02-25T00:54:54.968626Z"},"trusted":true},"execution_count":null,"outputs":[]}]}