{"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":"# 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$$\n\nSome neutrinos come from the sky above Antarctica and go downward. Others come from underground through the Earth and go upward.\nThis method cannot determine which direction the neutrino went on the regression line. For this reason, this notebook assumes that all neutrinos went downward.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom numba import njit\n\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:04:27.594320Z","iopub.execute_input":"2023-01-22T18:04:27.594775Z","iopub.status.idle":"2023-01-22T18:04:28.268107Z","shell.execute_reply.started":"2023-01-22T18:04:27.594686Z","shell.execute_reply":"2023-01-22T18:04:28.267231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@njit(\"Tuple((f8, f8))(f8[:], f8[:], f8[:])\", cache=True)\ndef compute_angle(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 # downward-going\n    #zr = -1 # upward-going\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-22T18:04:28.269886Z","iopub.execute_input":"2023-01-22T18:04:28.270486Z","iopub.status.idle":"2023-01-22T18:04:33.084467Z","shell.execute_reply.started":"2023-01-22T18:04:28.270454Z","shell.execute_reply":"2023-01-22T18:04:33.083317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv', index_col='sensor_id')","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:04:33.085853Z","iopub.execute_input":"2023-01-22T18:04:33.086224Z","iopub.status.idle":"2023-01-22T18:04:33.126392Z","shell.execute_reply.started":"2023-01-22T18:04:33.086180Z","shell.execute_reply":"2023-01-22T18:04:33.125265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_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-22T18:04:33.129446Z","iopub.execute_input":"2023-01-22T18:04:33.129924Z","iopub.status.idle":"2023-01-22T18:05:42.838123Z","shell.execute_reply.started":"2023-01-22T18:04:33.129878Z","shell.execute_reply":"2023-01-22T18:05:42.836245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mask = ~auxiliary\nx = position[:, 0]\ny = position[:, 1]\nz = position[:, 2]\n\npred_azimuth, pred_zenith = compute_angle(x[mask], y[mask], z[mask])\n\nprint(f\"azimuth:{azimuth} pred_azimuth:{pred_azimuth}\")\nprint(f\"zenith:{zenith} pred_zenith:{pred_zenith}\")","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:42.841065Z","iopub.execute_input":"2023-01-22T18:05:42.842350Z","iopub.status.idle":"2023-01-22T18:05:42.859603Z","shell.execute_reply.started":"2023-01-22T18:05:42.842291Z","shell.execute_reply":"2023-01-22T18:05:42.858079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ax = plt.figure(figsize=(16, 12)).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=-30, elev=30)\nax.scatter(sensor_geometry.x, sensor_geometry.y, sensor_geometry.z, s=0.3, color='black', alpha=0.3)\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\npx1 = xzmin\npy1 = yzmin\npz1 = zmin\npz2 = zmax\npr = (pz2-pz1)/np.cos(pred_zenith)\npx2 = pr*np.sin(pred_zenith)*np.cos(pred_azimuth) + px1\npy2 = pr*np.sin(pred_zenith)*np.sin(pred_azimuth) + py1\n\ntx1 = xzmin\nty1 = yzmin\ntz1 = zmin\ntz2 = zmax\ntr = (tz2-tz1)/np.cos(zenith)\ntx2 = tr*np.sin(zenith)*np.cos(azimuth) + tx1\nty2 = tr*np.sin(zenith)*np.sin(azimuth) + ty1\n\nax.scatter(x[~mask], y[~mask], z[~mask], c='black', s=100.0, alpha=0.1)\nax.scatter(x[mask], y[mask], z[mask], c='blue', s=100.0, alpha=0.7)\n\nax.plot([tx1, tx2], [ty1, ty2], [tz1, tz2], c='red', linewidth=3.0, label=\"true\")\nax.plot([px1, px2], [py1, py2], [pz1, pz2], c='orange', linewidth=3.0, label=\"pred\")\n\nax.legend()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:42.861593Z","iopub.execute_input":"2023-01-22T18:05:42.862086Z","iopub.status.idle":"2023-01-22T18:05:43.783112Z","shell.execute_reply.started":"2023-01-22T18:05:42.862036Z","shell.execute_reply":"2023-01-22T18:05:43.782191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del meta_df, batch_df, batch_features, event_features","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:43.784427Z","iopub.execute_input":"2023-01-22T18:05:43.785070Z","iopub.status.idle":"2023-01-22T18:05:44.042824Z","shell.execute_reply.started":"2023-01-22T18:05:43.785034Z","shell.execute_reply":"2023-01-22T18:05:44.041709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference","metadata":{}},{"cell_type":"code","source":"sub_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/sample_submission.parquet')\nsub_df = sub_df.set_index('event_id', drop=True)\nsub_df","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:44.044543Z","iopub.execute_input":"2023-01-22T18:05:44.045080Z","iopub.status.idle":"2023-01-22T18:05:44.069760Z","shell.execute_reply.started":"2023-01-22T18:05:44.044940Z","shell.execute_reply":"2023-01-22T18:05:44.068905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\nmeta_df","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:44.072394Z","iopub.execute_input":"2023-01-22T18:05:44.073215Z","iopub.status.idle":"2023-01-22T18:05:44.091097Z","shell.execute_reply.started":"2023-01-22T18:05:44.073182Z","shell.execute_reply":"2023-01-22T18:05:44.090198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for batch_id in 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(x[mask], y[mask], z[mask])\n        sub_df.loc[event_id, 'azimuth'] = pred_azimuth\n        sub_df.loc[event_id, 'zenith'] = pred_zenith","metadata":{"execution":{"iopub.status.busy":"2023-01-22T18:05:44.092515Z","iopub.execute_input":"2023-01-22T18:05:44.093130Z","iopub.status.idle":"2023-01-22T18:05:44.122758Z","shell.execute_reply.started":"2023-01-22T18:05:44.093096Z","shell.execute_reply":"2023-01-22T18:05:44.121463Z"},"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-22T18:05:44.124176Z","iopub.execute_input":"2023-01-22T18:05:44.125034Z","iopub.status.idle":"2023-01-22T18:05:44.137878Z","shell.execute_reply.started":"2023-01-22T18:05:44.124988Z","shell.execute_reply":"2023-01-22T18:05:44.136992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-21T09:16:00.162770Z","iopub.execute_input":"2023-01-21T09:16:00.163371Z","iopub.status.idle":"2023-01-21T09:16:00.177622Z","shell.execute_reply.started":"2023-01-21T09:16:00.163315Z","shell.execute_reply":"2023-01-21T09:16:00.175063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}