{"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":"code","source":"import numpy as np\nimport pandas as pd\nimport plotly.express as px\nfrom sklearn.decomposition import PCA\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nimport multiprocessing","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-22T11:24:02.441208Z","iopub.execute_input":"2023-01-22T11:24:02.441551Z","iopub.status.idle":"2023-01-22T11:24:02.447659Z","shell.execute_reply.started":"2023-01-22T11:24:02.441522Z","shell.execute_reply":"2023-01-22T11:24:02.445629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Update (version 8+)\nIt's better to always apply the PCA only on the non-auxiliary, regardless of their number. This change alone raises the score to VAL 1.23 LB 1.1218. Also, check [this notebook](https://www.kaggle.com/code/jirkaborovec/icecube-neutrino-fitting-3d-points-cloud) that uses skspatial: it gives very similar fitted lines but is a bit faster than PCA.","metadata":{}},{"cell_type":"markdown","source":"# TL;DR\nI perform a PCA on the pulses coordinates and take the first principal component as the neutrino's direction. The PCA is performed only on the non-auxiliary ('true') pulses if there are enough (>15); otherwise, I perform it on all the data. Then, I calculate the average time of the 'true' pulses on the negative and positive sides of the vector and flip its direction if it points to the wrong one. (lower times should be on the positive side since it is supposed to point to the source of the neutrino). This analysis is done only for the 'true' pulses since the auxiliary pulses times are too noisy for this analysis without further refinement. Validation score 1.282, LB score 1.274\n","metadata":{}},{"cell_type":"markdown","source":"# Train meta","metadata":{}},{"cell_type":"code","source":"train_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train_meta.parquet')\ntrain_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:23:27.372241Z","iopub.execute_input":"2023-01-22T11:23:27.372596Z","iopub.status.idle":"2023-01-22T11:24:02.439548Z","shell.execute_reply.started":"2023-01-22T11:23:27.372564Z","shell.execute_reply":"2023-01-22T11:24:02.438047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Total number of events: ' + str(len(train_meta)))","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:09.490607Z","iopub.execute_input":"2023-01-22T11:24:09.490965Z","iopub.status.idle":"2023-01-22T11:24:09.497175Z","shell.execute_reply.started":"2023-01-22T11:24:09.490935Z","shell.execute_reply":"2023-01-22T11:24:09.495699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train batch 1","metadata":{}},{"cell_type":"code","source":"train_batch_1 = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_1.parquet')\ntrain_batch_1.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:11.891862Z","iopub.execute_input":"2023-01-22T11:24:11.892580Z","iopub.status.idle":"2023-01-22T11:24:14.829571Z","shell.execute_reply.started":"2023-01-22T11:24:11.892550Z","shell.execute_reply":"2023-01-22T11:24:14.828531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Number of pulses: ' + str(len(train_batch_1)))\nprint('Number of events: ' + str(len(train_batch_1.index.unique())))","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:17.165906Z","iopub.execute_input":"2023-01-22T11:24:17.166277Z","iopub.status.idle":"2023-01-22T11:24:17.687353Z","shell.execute_reply.started":"2023-01-22T11:24:17.166246Z","shell.execute_reply":"2023-01-22T11:24:17.685986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sensor geometry","metadata":{}},{"cell_type":"code","source":"sensor_geometry = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv')\nsensor_geometry.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:23.044719Z","iopub.execute_input":"2023-01-22T11:24:23.045036Z","iopub.status.idle":"2023-01-22T11:24:23.069019Z","shell.execute_reply.started":"2023-01-22T11:24:23.045012Z","shell.execute_reply":"2023-01-22T11:24:23.067995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(sensor_geometry.x, sensor_geometry.y)","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:25.110209Z","iopub.execute_input":"2023-01-22T11:24:25.110696Z","iopub.status.idle":"2023-01-22T11:24:25.322193Z","shell.execute_reply.started":"2023-01-22T11:24:25.110652Z","shell.execute_reply":"2023-01-22T11:24:25.321319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_3d(sensor_geometry, x='x', y='y', z='z', opacity=0.5)\nfig.update_traces(marker_size=2)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:27.160973Z","iopub.execute_input":"2023-01-22T11:24:27.161325Z","iopub.status.idle":"2023-01-22T11:24:28.460528Z","shell.execute_reply.started":"2023-01-22T11:24:27.161297Z","shell.execute_reply":"2023-01-22T11:24:28.459472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test meta","metadata":{}},{"cell_type":"code","source":"test_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ntest_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ntest_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:24:33.212707Z","iopub.execute_input":"2023-01-22T11:24:33.213243Z","iopub.status.idle":"2023-01-22T11:24:33.242822Z","shell.execute_reply.started":"2023-01-22T11:24:33.213186Z","shell.execute_reply":"2023-01-22T11:24:33.241154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(test_meta))","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:25:02.582023Z","iopub.execute_input":"2023-01-22T11:25:02.583757Z","iopub.status.idle":"2023-01-22T11:25:02.590596Z","shell.execute_reply.started":"2023-01-22T11:25:02.583664Z","shell.execute_reply":"2023-01-22T11:25:02.589585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Understanding the data\n### I will explore event 67 (number three in train_meta) as a case event","metadata":{}},{"cell_type":"code","source":"case_event_idx = 3\ncase_event_pulses = train_batch_1.iloc[train_meta.iloc[case_event_idx].first_pulse_index.astype(int) : \n                                         train_meta.iloc[case_event_idx].last_pulse_index.astype(int)+1].copy()\nprint(len(case_event_pulses))\ncase_event_pulses.tail()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:25:50.516730Z","iopub.execute_input":"2023-01-22T11:25:50.517206Z","iopub.status.idle":"2023-01-22T11:25:50.537919Z","shell.execute_reply.started":"2023-01-22T11:25:50.517163Z","shell.execute_reply":"2023-01-22T11:25:50.535929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Finding the coordinates of each pulse from the sensor_geometry data and mean normalizing","metadata":{}},{"cell_type":"code","source":"case_event_pulses = case_event_pulses.reset_index()\ncase_event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[case_event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\n\ncase_event_pulses.x = case_event_pulses.x - case_event_pulses.x.mean()\ncase_event_pulses.y = case_event_pulses.y - case_event_pulses.y.mean()\ncase_event_pulses.z = case_event_pulses.z - case_event_pulses.z.mean()\ncase_event_pulses.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:26:59.108769Z","iopub.execute_input":"2023-01-22T11:26:59.109141Z","iopub.status.idle":"2023-01-22T11:26:59.132021Z","shell.execute_reply.started":"2023-01-22T11:26:59.109110Z","shell.execute_reply":"2023-01-22T11:26:59.130867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Calculating the vector of the case event from the target angles and preparing for visualisation","metadata":{}},{"cell_type":"code","source":"zenith_target = train_meta.iloc[case_event_idx].zenith\nazimuth_target = train_meta.iloc[case_event_idx].azimuth\n\nvector_target = [np.cos(azimuth_target) * np.sin(zenith_target), np.sin(azimuth_target) * np.sin(zenith_target),\n                 np.cos(zenith_target)]\n\nvector_base = np.array([-500, 500])\nx = vector_base*vector_target[0]\ny = vector_base*vector_target[1]\nz = vector_base*vector_target[2]\nvector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\nvector_target_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:27:09.781574Z","iopub.execute_input":"2023-01-22T11:27:09.782097Z","iopub.status.idle":"2023-01-22T11:27:09.802557Z","shell.execute_reply.started":"2023-01-22T11:27:09.782056Z","shell.execute_reply":"2023-01-22T11:27:09.799810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Many thanks to JIRKA BOROVEC for demonstrating plotly [here](https://www.kaggle.com/code/jirkaborovec/icecube-eda-3d-interactive-viewer). I was quite frustrated by the lack of pyplot interactivity in kaggle.\n","metadata":{}},{"cell_type":"code","source":"fig_auxiliary = px.scatter_3d(case_event_pulses.loc[case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\nfig_non_auxiliary = px.scatter_3d(case_event_pulses.loc[~case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\nfig_line = px.line_3d(vector_target_df, x=\"x\", y=\"y\", z=\"z\")\n\n\nfig_auxiliary.update_traces(marker_size=2)\nfig_non_auxiliary.update_traces(marker_size=2)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:27:18.373270Z","iopub.execute_input":"2023-01-22T11:27:18.373783Z","iopub.status.idle":"2023-01-22T11:27:18.548500Z","shell.execute_reply.started":"2023-01-22T11:27:18.373748Z","shell.execute_reply":"2023-01-22T11:27:18.546697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The figure demonstrates the alignment of the angles with the 'true' pulses. Unfortunately, many events (probably most) do not have a clear pulse line like this one (hence I chose it for demonstration). Try to change case_event_idx to other numbers, then rerun the last four cells and see for yourself :)\nAn important thing to note is that most non-auxiliary pulses appear as a short vertical line of close dots, not because they come straight from the top but due to the structure of the detector. Watch [this](https://www.youtube.com/watch?v=iv-Rz3-s4BM) video for a better explanation, and in particular, 3:12-3:25.\n","metadata":{}},{"cell_type":"markdown","source":"### Trying a naive fitting: first principal component","metadata":{}},{"cell_type":"code","source":"pca = PCA(n_components=1).fit(case_event_pulses[['x', 'y', 'z']])\nvector = pca.components_[0]\nvector_base = np.array([-500, 500])\nx = vector_base*vector[0]\ny = vector_base*vector[1]\nz = vector_base*vector[2]\nvector_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n\nfig_auxiliary = px.scatter_3d(case_event_pulses.loc[case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\nfig_non_auxiliary = px.scatter_3d(case_event_pulses.loc[~case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\nfig_line = px.line_3d(vector_target_df, x=\"x\", y=\"y\", z=\"z\")\nfig_fit_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\", color_discrete_sequence=['magenta'])\n\nfig_auxiliary.update_traces(marker_size=2)\nfig_non_auxiliary.update_traces(marker_size=2)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data + fig_fit_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:30:40.450327Z","iopub.execute_input":"2023-01-22T11:30:40.450754Z","iopub.status.idle":"2023-01-22T11:30:40.665354Z","shell.execute_reply.started":"2023-01-22T11:30:40.450724Z","shell.execute_reply":"2023-01-22T11:30:40.663601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Let's calculate the error using [this scheme](https://www.kaggle.com/code/sohier/mean-angular-error)","metadata":{}},{"cell_type":"code","source":"def angular_dist_score(az_true, zen_true, az_pred, zen_pred):\n    '''\n    calculate the MAE of the angular distance between two directions.\n    The two vectors are first converted to cartesian unit vectors,\n    and then their scalar product is computed, which is equal to\n    the cosine of the angle between the two vectors. The inverse \n    cosine (arccos) thereof is then the angle between the two input vectors\n    \n    Parameters:\n    -----------\n    \n    az_true : float (or array thereof)\n        true azimuth value(s) in radian\n    zen_true : float (or array thereof)\n        true zenith value(s) in radian\n    az_pred : float (or array thereof)\n        predicted azimuth value(s) in radian\n    zen_pred : float (or array thereof)\n        predicted zenith value(s) in radian\n    \n    Returns:\n    --------\n    \n    dist : float\n        mean over the angular distance(s) in radian\n    '''\n    \n    if not (np.all(np.isfinite(az_true)) and\n            np.all(np.isfinite(zen_true)) and\n            np.all(np.isfinite(az_pred)) and\n            np.all(np.isfinite(zen_pred))):\n        raise ValueError(\"All arguments must be finite\")\n    \n    # pre-compute all sine and cosine values\n    sa1 = np.sin(az_true)\n    ca1 = np.cos(az_true)\n    sz1 = np.sin(zen_true)\n    cz1 = np.cos(zen_true)\n    \n    sa2 = np.sin(az_pred)\n    ca2 = np.cos(az_pred)\n    sz2 = np.sin(zen_pred)\n    cz2 = np.cos(zen_pred)\n    \n    # scalar product of the two cartesian vectors (x = sz*ca, y = sz*sa, z = cz)\n    scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n    \n    # scalar product of two unit vectors is always between -1 and 1, this is against nummerical instability\n    # that might otherwise occure from the finite precision of the sine and cosine functions\n    scalar_prod =  np.clip(scalar_prod, -1, 1)\n    \n    # convert back to an angle (in radian)\n    return np.average(np.abs(np.arccos(scalar_prod)))","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:31:26.606880Z","iopub.execute_input":"2023-01-22T11:31:26.607372Z","iopub.status.idle":"2023-01-22T11:31:26.619918Z","shell.execute_reply.started":"2023-01-22T11:31:26.607331Z","shell.execute_reply":"2023-01-22T11:31:26.617855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zenith = np.arccos(vector[2])\nazimuth = np.arctan2(vector[1], vector[0])\nif azimuth<0:\n    azimuth = 2*np.pi + azimuth\n\nprint('Predictions:')\nprint(f'zenith: {zenith}')\nprint(f'azimuth: {azimuth}')\nprint('Target:')\nprint(f'zenith: {zenith_target}')\nprint(f'azimuth: {azimuth_target}')\nprint('Error:')\nprint(angular_dist_score(azimuth_target, zenith_target, azimuth, zenith))","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:31:31.039150Z","iopub.execute_input":"2023-01-22T11:31:31.039608Z","iopub.status.idle":"2023-01-22T11:31:31.049078Z","shell.execute_reply.started":"2023-01-22T11:31:31.039572Z","shell.execute_reply":"2023-01-22T11:31:31.047807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Baseline validation","metadata":{}},{"cell_type":"markdown","source":"The baseline strategy will be fitting a vector through only the 'true' pulses if there are enough (>15); otherwise, fitting a vector through all the data. Additionally, I will calculate the average time of the 'true' pulses on the negative and positive sides of the vector and flip its direction if it points to the wrong one. (lower times should be on the positive side since it is supposed to point to the source of the neutrino). This analysis is done only for the 'true' pulses since the auxiliary pulses times are too noisy for this analysis without further refinement.","metadata":{}},{"cell_type":"code","source":"train_batches = train_meta.batch_id.unique()\ntrain_meta_baseline = train_meta.loc[train_meta.batch_id == train_batches[0]].iloc[:10000].copy()\ntrain_meta_baseline.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:33:26.258468Z","iopub.execute_input":"2023-01-22T11:33:26.258984Z","iopub.status.idle":"2023-01-22T11:33:27.361432Z","shell.execute_reply.started":"2023-01-22T11:33:26.258940Z","shell.execute_reply":"2023-01-22T11:33:27.360122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_baseline = pd.read_parquet(f'/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_{str(train_batches[0])}.parquet')\ntrain_baseline = train_baseline.iloc[:train_meta_baseline.iloc[-1].last_pulse_index.astype(int)]\ntrain_baseline = train_baseline.reset_index()\ntrain_baseline[['x', 'y', 'z']] = sensor_geometry.loc[train_baseline.sensor_id].reset_index()[['x', 'y', 'z']]\n\ntrain_baseline['x'] = train_baseline['x'] - train_baseline.groupby('event_id').x.transform('mean')\ntrain_baseline['y'] = train_baseline['y'] - train_baseline.groupby('event_id').y.transform('mean')\ntrain_baseline['z'] = train_baseline['z'] - train_baseline.groupby('event_id').z.transform('mean')\n\nprint(len(train_baseline))\ntrain_baseline.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:33:33.547851Z","iopub.execute_input":"2023-01-22T11:33:33.548249Z","iopub.status.idle":"2023-01-22T11:33:35.820209Z","shell.execute_reply.started":"2023-01-22T11:33:33.548217Z","shell.execute_reply":"2023-01-22T11:33:35.818544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:33:38.999431Z","iopub.execute_input":"2023-01-22T11:33:39.000147Z","iopub.status.idle":"2023-01-22T11:33:39.171062Z","shell.execute_reply.started":"2023-01-22T11:33:39.000100Z","shell.execute_reply":"2023-01-22T11:33:39.168971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_angles(event_df):\n\n    try:\n        pca = PCA(n_components=2).fit(event_df.loc[~event_df.auxiliary][['x', 'y', 'z']])\n        vector = pca.components_[0]\n\n        # In order to perform time direction analysis on the non-auxiliary ('true') pulses, we neeed to mean-normalize them separately\n        event_df_non_auxiliary = event_df.loc[~event_df.auxiliary].copy()\n        event_df_non_auxiliary.x = event_df_non_auxiliary.x - event_df_non_auxiliary.x.mean()\n        event_df_non_auxiliary.y = event_df_non_auxiliary.y - event_df_non_auxiliary.y.mean()\n        event_df_non_auxiliary.z = event_df_non_auxiliary.z - event_df_non_auxiliary.z.mean()\n\n        # This is the formula for calculating a distance of a point from a plane. The plane direction is the first principal component, and it goes through the origin.\n        event_df_non_auxiliary['distance'] = (event_df_non_auxiliary.x*vector[0]+\n                                              event_df_non_auxiliary.y*vector[1]+\n                                              event_df_non_auxiliary.z*vector[2])/np.linalg.norm(vector)\n\n        # Flip the vector direction if it points away from the neutrino origin\n        if (event_df_non_auxiliary.loc[(~event_df_non_auxiliary.auxiliary)&(event_df_non_auxiliary.distance>0)].time.mean() >\n            event_df_non_auxiliary.loc[(~event_df_non_auxiliary.auxiliary)&(event_df_non_auxiliary.distance<0)].time.mean()):\n            vector = -1*vector\n\n        # This is important since minor numerical deviations can give a vector component like 1.00001, and then arccos fail.\n        vector = np.clip(vector, -1, 1)\n\n        zenith = np.arccos(vector[2])\n        azimuth = np.arctan2(vector[1], vector[0])\n        if azimuth<0:\n            azimuth = 2*np.pi + azimuth\n    except:\n        zenith, azimuth = 0., 0.\n        \n    return zenith, azimuth","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:34:10.761912Z","iopub.execute_input":"2023-01-22T11:34:10.762254Z","iopub.status.idle":"2023-01-22T11:34:10.772724Z","shell.execute_reply.started":"2023-01-22T11:34:10.762229Z","shell.execute_reply":"2023-01-22T11:34:10.771698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ntrain_meta_baseline[['azimuth_pred', 'zenith_pred']] = None\nfor i in range(len(train_meta_baseline)):\n    row = train_meta_baseline.iloc[i]\n    event_df = train_baseline.iloc[row.first_pulse_index.astype(int):row.last_pulse_index.astype(int)+1].copy()\n    zenith, azimuth = get_angles(event_df)\n    train_meta_baseline.loc[i, 'zenith_pred'] = zenith\n    train_meta_baseline.loc[i, 'azimuth_pred'] = azimuth\n    \ntrain_meta_baseline.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:34:57.874432Z","iopub.execute_input":"2023-01-22T11:34:57.874849Z","iopub.status.idle":"2023-01-22T11:36:27.044782Z","shell.execute_reply.started":"2023-01-22T11:34:57.874818Z","shell.execute_reply":"2023-01-22T11:36:27.042746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"angular_dist_score(train_meta_baseline.azimuth.to_numpy(), train_meta_baseline.zenith.to_numpy(),\n                   train_meta_baseline.azimuth_pred.astype(float).to_numpy(),\n                   train_meta_baseline.astype(float).zenith_pred.to_numpy())","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:45:54.451694Z","iopub.execute_input":"2023-01-22T11:45:54.453889Z","iopub.status.idle":"2023-01-22T11:45:54.472201Z","shell.execute_reply.started":"2023-01-22T11:45:54.453826Z","shell.execute_reply":"2023-01-22T11:45:54.469300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_meta_baseline,train_baseline, train_batch_1, train_meta\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:45:57.342853Z","iopub.execute_input":"2023-01-22T11:45:57.343312Z","iopub.status.idle":"2023-01-22T11:45:57.671161Z","shell.execute_reply.started":"2023-01-22T11:45:57.343277Z","shell.execute_reply":"2023-01-22T11:45:57.669134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Making a submission","metadata":{}},{"cell_type":"code","source":"sample_submission = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/sample_submission.parquet')\nsample_submission.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:47:06.374469Z","iopub.execute_input":"2023-01-22T11:47:06.374875Z","iopub.status.idle":"2023-01-22T11:47:06.406636Z","shell.execute_reply.started":"2023-01-22T11:47:06.374848Z","shell.execute_reply":"2023-01-22T11:47:06.405022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ntest_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:47:08.326886Z","iopub.execute_input":"2023-01-22T11:47:08.327754Z","iopub.status.idle":"2023-01-22T11:47:08.349390Z","shell.execute_reply.started":"2023-01-22T11:47:08.327708Z","shell.execute_reply":"2023-01-22T11:47:08.348374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_test_batch(test_batch_id):\n    test_batch = pd.read_parquet(f'/kaggle/input/icecube-neutrinos-in-deep-ice/test/batch_{str(int(test_batch_id))}.parquet')\n    test_batch = test_batch.reset_index()\n    test_batch[['x', 'y', 'z']] = sensor_geometry.loc[test_batch.sensor_id].reset_index()[['x', 'y', 'z']]\n    test_batch['x'] = test_batch['x'] - test_batch.groupby('event_id').x.transform('mean')\n    test_batch['y'] = test_batch['y'] - test_batch.groupby('event_id').y.transform('mean')\n    test_batch['z'] = test_batch['z'] - test_batch.groupby('event_id').z.transform('mean')\n    return test_batch","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:47:12.038981Z","iopub.execute_input":"2023-01-22T11:47:12.039415Z","iopub.status.idle":"2023-01-22T11:47:12.049009Z","shell.execute_reply.started":"2023-01-22T11:47:12.039364Z","shell.execute_reply":"2023-01-22T11:47:12.047184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batches = test_meta.batch_id.unique()\nfor batch in batches:\n    test_batch = load_test_batch(batch)\n    test_meta_batch = test_meta.loc[test_meta.batch_id == batch]\n    \n    # multiprocessing ftw\n    items = []\n    for i in range(len(test_meta_batch)):\n        row = test_meta_batch.iloc[i]\n        items.append(test_batch.iloc[row.first_pulse_index.astype(int):row.last_pulse_index.astype(int)+1].copy())\n    results = []\n    with multiprocessing.Pool() as pool:\n        for result in pool.map(get_angles, items):\n            results.append(result)\n            \n    test_meta.loc[test_meta_batch.index, 'zenith'] = [results[i][0] for i in range(len(results))]\n    test_meta.loc[test_meta_batch.index, 'azimuth'] = [results[i][1] for i in range(len(results))]\n    \n    del items, test_batch, test_meta_batch,results\n    gc.collect()\ntest_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:51:36.294809Z","iopub.execute_input":"2023-01-22T11:51:36.295326Z","iopub.status.idle":"2023-01-22T11:51:37.036326Z","shell.execute_reply.started":"2023-01-22T11:51:36.295290Z","shell.execute_reply":"2023-01-22T11:51:37.034032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = test_meta[['event_id', 'azimuth', 'zenith']]\nsubmission.to_csv('submission.csv', index=False)\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-22T11:51:39.645187Z","iopub.execute_input":"2023-01-22T11:51:39.645807Z","iopub.status.idle":"2023-01-22T11:51:39.669422Z","shell.execute_reply.started":"2023-01-22T11:51:39.645763Z","shell.execute_reply":"2023-01-22T11:51:39.666962Z"},"trusted":true},"execution_count":null,"outputs":[]}]}