{"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":"This notebook builds upon the excellent work of https://www.kaggle.com/code/gyozzza/icecube-3d-linear-regression","metadata":{}},{"cell_type":"markdown","source":"## Method\nFitting the relationship between $x$ and $z$ and the relationship between $y$ and $z$ with straight lines.\n\n$$\nx = m_x z + c_x \\\\\ny = m_y z + c_y\n$$\n\nCalculating azimuthal angle $\\phi$ and zenith angle $\\theta$ using the slopes $m_x$, $m_y$.\n\n$$\n\\phi = \\arctan \\frac{m_y}{m_x} \\\\\n\\theta = \\arccos \\frac{1}{\\sqrt{1+m_x^2+m_y^2}}\n$$","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom tqdm.auto import tqdm\nfrom numba import njit\nimport gc\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:41:33.673045Z","iopub.execute_input":"2023-01-22T13:41:33.673479Z","iopub.status.idle":"2023-01-22T13:41:33.693425Z","shell.execute_reply.started":"2023-01-22T13:41:33.673445Z","shell.execute_reply":"2023-01-22T13:41:33.692110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load a single batch","metadata":{}},{"cell_type":"code","source":"%%time\nsensor_geometry = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv', index_col='sensor_id')\nmeta_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train_meta.parquet')\nbatch_id = meta_df.batch_id.unique()[0]\n\nbatch_df = meta_df[meta_df.batch_id == batch_id].set_index('event_id', drop=True)\nbatch_features = pd.read_parquet(f'/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_{batch_id}.parquet')\n\nevent_id = batch_df.index[6]\n\nevent_features = batch_features.iloc[batch_df.loc[event_id, 'first_pulse_index']:batch_df.loc[event_id, 'last_pulse_index']+1]\nazimuth = batch_df.loc[event_id, 'azimuth']\nzenith = batch_df.loc[event_id, 'zenith']\n\nposition = sensor_geometry.loc[event_features.sensor_id].values\ntime = event_features.time.values\ncharge = event_features.charge.values\nauxiliary = event_features.auxiliary.values","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:38:26.279585Z","iopub.execute_input":"2023-01-22T13:38:26.280262Z","iopub.status.idle":"2023-01-22T13:39:11.701046Z","shell.execute_reply.started":"2023-01-22T13:38:26.280224Z","shell.execute_reply":"2023-01-22T13:39:11.699903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Numba compiled angle function","metadata":{}},{"cell_type":"code","source":"@njit(cache=True)\ndef compute_angle_numba(x, y, z):\n    covx = np.cov(z, x)\n    mx = covx[0,1]/covx[0,0]\n    covy = np.cov(z, y)\n    my = covy[0,1]/covy[0,0]\n    zr = 1\n    xr = mx*zr\n    yr = my*zr\n    r = np.sqrt(zr**2+xr**2+yr**2)\n    azimuth = np.arctan2(yr, xr)\n    if azimuth < 0:\n        azimuth = 2*np.pi + azimuth\n    zenith = np.arccos(zr/r)\n    return azimuth, zenith","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:39:11.703272Z","iopub.execute_input":"2023-01-22T13:39:11.703697Z","iopub.status.idle":"2023-01-22T13:39:12.008463Z","shell.execute_reply.started":"2023-01-22T13:39:11.703654Z","shell.execute_reply":"2023-01-22T13:39:12.007397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test on a single batch","metadata":{}},{"cell_type":"code","source":"mask = ~auxiliary\nx = position[:, 0]\ny = position[:, 1]\nz = position[:, 2]\n\npred_azimuth, pred_zenith = compute_angle_numba(x[mask], y[mask], z[mask])\nazimuth_error = azimuth - pred_azimuth\nzenith_error = zenith - pred_zenith\n\nprint(f\"azimuth: {azimuth:.4f} pred_azimuth: {pred_azimuth:.4f}\")\nprint(f\"zenith: {zenith:.4f} pred_zenith: {pred_zenith:.4f}\")\nprint(f'Azimuth Error: {azimuth_error:.4f}\\nZenith Error: {zenith_error:.4f}')","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:39:27.894506Z","iopub.execute_input":"2023-01-22T13:39:27.894905Z","iopub.status.idle":"2023-01-22T13:39:32.656880Z","shell.execute_reply.started":"2023-01-22T13:39:27.894874Z","shell.execute_reply":"2023-01-22T13:39:32.656044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualise Single Batch","metadata":{}},{"cell_type":"code","source":"ax = plt.figure(figsize=(10, 10)).add_subplot(projection='3d')\n\nax.set_xlabel('x')\nax.set_ylabel('y')\nax.set_zlabel('z')\nax.set_box_aspect([1,1,1])\nax.view_init(azim=-35, elev=25)\nax.scatter(sensor_geometry.x, sensor_geometry.y, sensor_geometry.z, s=0.3, color='black', alpha=0.2)\n\nzmax = sensor_geometry.z.max()\nzmin = sensor_geometry.z.min()\n\ncovx = np.cov(z[mask], x[mask])\nmx = covx[0,1]/covx[0,0]\ncx = x[mask].mean() - mx*z[mask].mean()\ncovy = np.cov(z[mask], y[mask])\nmy = covy[0,1]/covy[0,0]\ncy = y[mask].mean() - my*z[mask].mean()\n\nxzmax = mx*zmax + cx\nyzmax = my*zmax + cy\nxzmin = mx*zmin + cx\nyzmin = my*zmin + cy\n\nzrange = zmax - zmin\n\npr = (zrange)/np.cos(pred_zenith)\npred_x = pr*np.sin(pred_zenith)*np.cos(pred_azimuth) + xzmin\npred_y = pr*np.sin(pred_zenith)*np.sin(pred_azimuth) + yzmin\n\ntr = (zrange)/np.cos(zenith)\ntrue_x = tr*np.sin(zenith)*np.cos(azimuth) + xzmin\ntrue_y = tr*np.sin(zenith)*np.sin(azimuth) + yzmin\n\nax.scatter(x[~mask], y[~mask], z[~mask], c='black', s=100.0, alpha=0.1, label='Auxiliary')\nax.scatter(x[mask], y[mask], z[mask], c='blue', s=100.0, alpha=0.7, label='Digitised Hits')\n\nax.plot([xzmin, true_x], [yzmin, true_y], [zmin, zmax], c='red', linewidth=3.0, label=\"True Neutrino Path\")\nax.plot([xzmin, pred_x], [yzmin, pred_y], [zmin, zmax], c='orange', linewidth=3.0, label=\"Predicted Neutrino Path\")\n\nax.legend()\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:41:07.899626Z","iopub.execute_input":"2023-01-22T13:41:07.900060Z","iopub.status.idle":"2023-01-22T13:41:08.444001Z","shell.execute_reply.started":"2023-01-22T13:41:07.900022Z","shell.execute_reply":"2023-01-22T13:41:08.442756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del meta_df, batch_df, batch_features, event_features\n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:41:35.610026Z","iopub.execute_input":"2023-01-22T13:41:35.610437Z","iopub.status.idle":"2023-01-22T13:41:36.202331Z","shell.execute_reply.started":"2023-01-22T13:41:35.610405Z","shell.execute_reply":"2023-01-22T13:41:36.200775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Run Inference","metadata":{}},{"cell_type":"code","source":"%%time\nsub_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/sample_submission.parquet')\nsub_df = sub_df.set_index('event_id', drop=True)\ndisplay(sub_df)\n\nmeta_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ndisplay(meta_df)","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:41:59.503667Z","iopub.execute_input":"2023-01-22T13:41:59.504136Z","iopub.status.idle":"2023-01-22T13:41:59.542681Z","shell.execute_reply.started":"2023-01-22T13:41:59.504098Z","shell.execute_reply":"2023-01-22T13:41:59.541763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfor batch_id in tqdm(meta_df.batch_id.unique()):\n    batch_df = meta_df[meta_df.batch_id == batch_id].set_index('event_id', drop=True)\n    batch_features = pd.read_parquet(f'/kaggle/input/icecube-neutrinos-in-deep-ice/test/batch_{batch_id}.parquet')\n    for event_id in batch_df.index:\n        event_features = batch_features.iloc[batch_df.loc[event_id, 'first_pulse_index']:batch_df.loc[event_id, 'last_pulse_index']+1]\n        position = sensor_geometry.loc[event_features.sensor_id].values\n        time = event_features.time.values\n        charge = event_features.charge.values\n        auxiliary = event_features.auxiliary.values\n        \n        mask = ~auxiliary\n        if np.unique(position[mask], axis=0).shape[0] <= 1:\n            sub_df.loc[event_id, 'azimuth'] = np.nan\n            sub_df.loc[event_id, 'zenith'] = np.nan\n            continue\n        x = position[:, 0]\n        y = position[:, 1]\n        z = position[:, 2]\n        pred_azimuth, pred_zenith = compute_angle_numba(x[mask], y[mask], z[mask])\n        sub_df.loc[event_id, 'azimuth'] = pred_azimuth\n        sub_df.loc[event_id, 'zenith'] = pred_zenith\n        \n    del batch_df\n    _ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:42:22.844088Z","iopub.execute_input":"2023-01-22T13:42:22.844459Z","iopub.status.idle":"2023-01-22T13:42:23.072959Z","shell.execute_reply.started":"2023-01-22T13:42:22.844430Z","shell.execute_reply":"2023-01-22T13:42:23.071865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_df = sub_df.fillna(sub_df.mean())\nsub_df","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:42:37.248399Z","iopub.execute_input":"2023-01-22T13:42:37.248787Z","iopub.status.idle":"2023-01-22T13:42:37.261813Z","shell.execute_reply.started":"2023-01-22T13:42:37.248756Z","shell.execute_reply":"2023-01-22T13:42:37.260907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-22T13:42:40.243879Z","iopub.execute_input":"2023-01-22T13:42:40.244275Z","iopub.status.idle":"2023-01-22T13:42:40.251327Z","shell.execute_reply.started":"2023-01-22T13:42:40.244244Z","shell.execute_reply":"2023-01-22T13:42:40.250355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}