{"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 Neutrino Path Detector\n\nThe IceCube Neutrino Detector is a scientific instrument in the South Pole which gathers data about neutrino events. Neutrinos are very hard to capture and analyze, and the objective of this notebook is to classify the direction that a neutrino passed through the IceCube detector, given data about sensors that have detected it.","metadata":{}},{"cell_type":"markdown","source":"### Possible things that would work\n* XGBoost\n* SVD\n* PCA\n* RANSAC\n* Least Mean Squares\n* Least Median Squares","metadata":{}},{"cell_type":"markdown","source":"## Strategy\n\n1. I compute the centroid of the activated sensor readings and use it as the origin of a local coordinate system. \n2. Then, I compute the vector that best fits the sensors by selecting the first component of PCA\n    * The first principal component of the sensor readings represents the direction of maximum variance in the points, which is a great fit for the data\n3. I find the direction of the vector (whether the neutrino came from the ground or from above) using the sensor times (TODO)\n3. Finally, I use spherical coordinate conversions to compute the zenith and azimuth angles of the direction vector.    ","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import axes3d\n# %matplotlib notebook\nimport glob\nfrom IPython.display import clear_output\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.preprocessing import Normalizer\nfrom sklearn.decomposition import PCA\nimport math\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:30:35.165778Z","iopub.execute_input":"2023-02-03T01:30:35.166244Z","iopub.status.idle":"2023-02-03T01:30:36.504205Z","shell.execute_reply.started":"2023-02-03T01:30:35.166166Z","shell.execute_reply":"2023-02-03T01:30:36.502848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### EDA","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-02-03T01:30:36.506423Z","iopub.execute_input":"2023-02-03T01:30:36.507042Z","iopub.status.idle":"2023-02-03T01:30:36.554096Z","shell.execute_reply.started":"2023-02-03T01:30:36.507003Z","shell.execute_reply":"2023-02-03T01:30:36.553391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(4, 3))\nax = fig.add_subplot()\nax.scatter(sensor_geometry['x'], sensor_geometry['y'])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:30:36.555422Z","iopub.execute_input":"2023-02-03T01:30:36.555985Z","iopub.status.idle":"2023-02-03T01:30:36.750800Z","shell.execute_reply.started":"2023-02-03T01:30:36.555942Z","shell.execute_reply":"2023-02-03T01:30:36.749809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(5, 5))\nax = fig.add_subplot(projection='3d')\nax.scatter(sensor_geometry['x'], sensor_geometry['y'], sensor_geometry['z'], alpha=.25, s=2)\nplt.title(\"Neutrino sensor locations\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:30:36.753547Z","iopub.execute_input":"2023-02-03T01:30:36.754181Z","iopub.status.idle":"2023-02-03T01:30:36.996090Z","shell.execute_reply.started":"2023-02-03T01:30:36.754136Z","shell.execute_reply":"2023-02-03T01:30:36.995036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-02-03T01:30:36.997218Z","iopub.execute_input":"2023-02-03T01:30:36.997572Z","iopub.status.idle":"2023-02-03T01:31:17.357915Z","shell.execute_reply.started":"2023-02-03T01:30:36.997535Z","shell.execute_reply":"2023-02-03T01:31:17.356601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"first_batch = pd.read_parquet(\"/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_1.parquet\")\nall_non_auxiliary_events = first_batch.loc[first_batch[\"auxiliary\"] == False].reset_index()\nall_non_auxiliary_events.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:17.359143Z","iopub.execute_input":"2023-02-03T01:31:17.359440Z","iopub.status.idle":"2023-02-03T01:31:21.381667Z","shell.execute_reply.started":"2023-02-03T01:31:17.359414Z","shell.execute_reply":"2023-02-03T01:31:21.380712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This will give us all the uniqe events_ids that correspond to non-auxiliary events\nall_non_auxiliary_events['event_id'].unique()[:10]","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:21.383462Z","iopub.execute_input":"2023-02-03T01:31:21.383936Z","iopub.status.idle":"2023-02-03T01:31:21.521309Z","shell.execute_reply.started":"2023-02-03T01:31:21.383898Z","shell.execute_reply":"2023-02-03T01:31:21.519331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# compress columns into one row\nfirst_batch = first_batch.reset_index()\n# remove uneccesary headers\n# first_batch = first_batch.drop('index', axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:21.523792Z","iopub.execute_input":"2023-02-03T01:31:21.524322Z","iopub.status.idle":"2023-02-03T01:31:21.745794Z","shell.execute_reply.started":"2023-02-03T01:31:21.524288Z","shell.execute_reply":"2023-02-03T01:31:21.744161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lets visualize some events","metadata":{}},{"cell_type":"code","source":"uniqe_events = all_non_auxiliary_events['event_id'].unique()\nnum_events_to_visualize = 6\nrandom_event_nums = [uniqe_events[np.random.randint(len(uniqe_events), size=1)][0] for _ in range(num_events_to_visualize)]","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:21.747329Z","iopub.execute_input":"2023-02-03T01:31:21.747726Z","iopub.status.idle":"2023-02-03T01:31:21.945515Z","shell.execute_reply.started":"2023-02-03T01:31:21.747694Z","shell.execute_reply":"2023-02-03T01:31:21.943593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(num_events_to_visualize//2, 2, subplot_kw={'projection': '3d'}, figsize=(15, 15))\naxs = axs.reshape((1, -1))[0]\n\nfor idx, event_num in enumerate(random_event_nums):\n    event_pulses = first_batch.loc[first_batch[\"event_id\"] == event_num]\n    event_meta = train_meta.loc[train_meta[\"event_id\"] == event_num]\n    \n    # Calculate Neutrino direction line coordinates based off label azimuth and zenith\n    zenith_target = event_meta.zenith.values[0]\n    azimuth_target = event_meta.azimuth.values[0]\n    vector_target = [np.cos(azimuth_target) * np.sin(zenith_target), \n                     np.sin(azimuth_target) * np.sin(zenith_target),\n                     np.cos(zenith_target)]\n\n    vector_base = np.array([-500, 500])\n    x = vector_base * vector_target[0]\n    y = vector_base * vector_target[1]\n    z = vector_base * vector_target[2]\n    vector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n    \n    # Add sensor ID position data to the event\n    event_pulses = event_pulses.reset_index()\n    event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\n    \n    auxiliary_pulses = event_pulses.loc[event_pulses['auxiliary'] == True]\n    non_auxiliary_pulses = event_pulses.loc[event_pulses['auxiliary'] == False]\n    \n    # Plot\n    ax = axs[idx]\n\n    ax.scatter(sensor_geometry['x'], sensor_geometry['y'], sensor_geometry['z'], alpha=.25, s=1, color=\"gray\")\n    \n    ax.scatter(auxiliary_pulses['x'], auxiliary_pulses['y'], auxiliary_pulses['z'],\n               c=\"red\", label=\"Auxiliary\")\n    \n    ax.scatter(non_auxiliary_pulses['x'], non_auxiliary_pulses['y'], non_auxiliary_pulses['z'],\n               c=non_auxiliary_pulses['time'], label=\"Not Auxiliary\")\n    \n    ax.plot(vector_target_df['x'], vector_target_df['y'], vector_target_df['z'],\n            c=\"green\", alpha=1, label=\"Neutrino Path\")\n\n    ax.set_title(f\"Sample Neutrino Event {event_num}\")\n    ax.legend()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:21.950822Z","iopub.execute_input":"2023-02-03T01:31:21.952241Z","iopub.status.idle":"2023-02-03T01:31:26.598192Z","shell.execute_reply.started":"2023-02-03T01:31:21.952198Z","shell.execute_reply":"2023-02-03T01:31:26.597011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Ok lets test the algorithm on a single event to see how it works","metadata":{}},{"cell_type":"code","source":"# event_num = 1952547\n# event_num = 3129713\n# event_num=67\n# event_num=2379627\n# event_num=2647853\nevent_num=24","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:26.599943Z","iopub.execute_input":"2023-02-03T01:31:26.600506Z","iopub.status.idle":"2023-02-03T01:31:26.607102Z","shell.execute_reply.started":"2023-02-03T01:31:26.600469Z","shell.execute_reply":"2023-02-03T01:31:26.605283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_event(event_num, plot_sensors=True, plot_auxiliary=True, plot_color=False):\n    event_pulses = first_batch.loc[first_batch[\"event_id\"] == event_num]\n\n    # Calculate Neutrino direction line coordinates based off label azimuth and zenith\n    zenith_target = train_meta.loc[train_meta[\"event_id\"] == event_num].zenith.values[0]\n    azimuth_target = train_meta.loc[train_meta[\"event_id\"] == event_num].azimuth.values[0]\n    vector_target = [np.cos(azimuth_target) * np.sin(zenith_target), \n                     np.sin(azimuth_target) * np.sin(zenith_target),\n                     np.cos(zenith_target)]\n\n    vector_base = np.array([-500, 500])\n    x = vector_base * vector_target[0]\n    y = vector_base * vector_target[1]\n    z = vector_base * vector_target[2]\n    vector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n\n    # Add sensor ID position data to the event\n    event_pulses = event_pulses.reset_index()\n    event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\n\n    auxiliary_pulses = event_pulses.loc[event_pulses['auxiliary'] == True]\n    non_auxiliary_pulses = event_pulses.loc[event_pulses['auxiliary'] == False]\n    \n    # Plot\n    fig = plt.figure(figsize=(5, 5))\n    ax = fig.add_subplot(projection='3d')\n    \n    if plot_sensors == True:\n        ax.scatter(sensor_geometry['x'], sensor_geometry['y'], sensor_geometry['z'], alpha=.25, s=1, color=\"gray\")\n    \n    if plot_auxiliary == True:\n        ax.scatter(auxiliary_pulses['x'], auxiliary_pulses['y'], auxiliary_pulses['z'],\n               c=\"gray\", label=\"Auxiliary\")\n\n    if plot_color == True:\n        ax.scatter(non_auxiliary_pulses['x'], non_auxiliary_pulses['y'], non_auxiliary_pulses['z'],\n               c=non_auxiliary_pulses['time'], label=\"Not Auxiliary\")\n    else: \n        ax.scatter(non_auxiliary_pulses['x'], non_auxiliary_pulses['y'], non_auxiliary_pulses['z'],\n               c=\"blue\", label=\"Not Auxiliary\")\n\n    ax.plot(vector_target_df['x'], vector_target_df['y'], vector_target_df['z'],\n            c=\"green\", alpha=1, label=\"Neutrino Path\")\n\n    ax.set_title(f\"Sample Neutrino Event {event_num}\")\n    ax.legend()\n    \n    return ax\n    \nplot_event(event_num, plot_auxiliary=False, plot_color=True)","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:26.609257Z","iopub.execute_input":"2023-02-03T01:31:26.610049Z","iopub.status.idle":"2023-02-03T01:31:27.300297Z","shell.execute_reply.started":"2023-02-03T01:31:26.610004Z","shell.execute_reply":"2023-02-03T01:31:27.298673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses = first_batch.loc[first_batch[\"event_id\"] == event_num]\nevent_pulses = event_pulses.reset_index()\nevent_pulses[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\nnon_auxiliary_pulses = event_pulses.loc[event_pulses['auxiliary'] == False]\nnon_auxiliary_pulses.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:27.302027Z","iopub.execute_input":"2023-02-03T01:31:27.302476Z","iopub.status.idle":"2023-02-03T01:31:27.355248Z","shell.execute_reply.started":"2023-02-03T01:31:27.302438Z","shell.execute_reply":"2023-02-03T01:31:27.353987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def spherical_to_cartesian(azimuth, zenith):\n    x = np.cos(azimuth) * np.sin(zenith)\n    y = np.sin(azimuth) * np.sin(zenith)\n    z = np.cos(zenith)\n    return [x, y, z]","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:27.356505Z","iopub.execute_input":"2023-02-03T01:31:27.356802Z","iopub.status.idle":"2023-02-03T01:31:27.364237Z","shell.execute_reply.started":"2023-02-03T01:31:27.356774Z","shell.execute_reply":"2023-02-03T01:31:27.363092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 1: Compute Centroid","metadata":{}},{"cell_type":"code","source":"# The centroid is simply the mean of the points in each dimension\npoints = non_auxiliary_pulses[['x', 'y', 'z']].to_numpy()\ncentroid = np.mean(points, axis=0)\n\nax = plot_event(event_num, plot_sensors=False, plot_auxiliary=False, plot_color=True)\nax.scatter(centroid[0], centroid[1], centroid[2], s=100, c=\"red\", label=\"Centroid\", alpha=1)\nax.legend()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:27.365649Z","iopub.execute_input":"2023-02-03T01:31:27.365991Z","iopub.status.idle":"2023-02-03T01:31:27.893411Z","shell.execute_reply.started":"2023-02-03T01:31:27.365961Z","shell.execute_reply":"2023-02-03T01:31:27.891931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 2: Center non-auxiliary points around the centroid","metadata":{}},{"cell_type":"code","source":"# To center points around the centroid, subtract the centroid's coordinates from the points' corrdinates\npoints_centered = points - centroid\ncentroid_centered = [0, 0, 0]\nfig = plt.figure(figsize=(5, 5))\nax = fig.add_subplot(projection='3d')\nax.scatter(points_centered[:, 0], points_centered[:, 1], points_centered[:, 2], color=\"blue\")\nax.scatter(centroid_centered[0], centroid_centered[1], centroid_centered[2], s=100, c=\"red\", label=\"Centroid\", alpha=1)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:27.895409Z","iopub.execute_input":"2023-02-03T01:31:27.896368Z","iopub.status.idle":"2023-02-03T01:31:28.075000Z","shell.execute_reply.started":"2023-02-03T01:31:27.896329Z","shell.execute_reply":"2023-02-03T01:31:28.072961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step 3: Select first principle component as line of best fit","metadata":{}},{"cell_type":"code","source":"def cartesian_to_spherical(x, y, z):\n    r = math.sqrt(x**2 + y**2 + z**2)\n    zenith = math.acos(z / r)\n    azimuth = math.atan2(y, x)\n    return azimuth, zenith\n\npca = PCA(n_components=1).fit(points_centered)\nvector = -1 * pca.components_[0]\n\nzenith_target = train_meta.loc[train_meta[\"event_id\"] == event_num].zenith.values[0]\nazimuth_target = train_meta.loc[train_meta[\"event_id\"] == event_num].azimuth.values[0]\n\nprint(cartesian_to_spherical(*vector))\nprint(azimuth_target, zenith_target)\n\nvector_base = np.array([-500, 500])\n\nx = vector_base * vector[0]\ny = vector_base * vector[1]\nz = vector_base * vector[2]\navg_vector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n\nvector_target = spherical_to_cartesian(azimuth_target, zenith_target)\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})\n\nfig = plt.figure(figsize=(5, 5))\nax = fig.add_subplot(projection='3d')\nax.scatter(points_centered[:, 0], points_centered[:, 1], points_centered[:, 2], color=\"blue\")\nax.plot(vector_target_df['x'], vector_target_df['y'], vector_target_df['z'], color=\"blue\", label=\"True\")\nax.plot(avg_vector_target_df['x'], avg_vector_target_df['y'], avg_vector_target_df['z'], color=\"red\", label=\"Pred\")\n\nax.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:54.772311Z","iopub.execute_input":"2023-02-03T01:31:54.774570Z","iopub.status.idle":"2023-02-03T01:31:55.395444Z","shell.execute_reply.started":"2023-02-03T01:31:54.774475Z","shell.execute_reply":"2023-02-03T01:31:55.393692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Testing on more data to find the score","metadata":{}},{"cell_type":"code","source":"# Source: https://www.kaggle.com/code/sohier/mean-angular-error\ndef 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    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-02-03T01:31:59.646556Z","iopub.execute_input":"2023-02-03T01:31:59.648258Z","iopub.status.idle":"2023-02-03T01:31:59.659995Z","shell.execute_reply.started":"2023-02-03T01:31:59.648206Z","shell.execute_reply":"2023-02-03T01:31:59.658554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:32:00.041208Z","iopub.execute_input":"2023-02-03T01:32:00.041631Z","iopub.status.idle":"2023-02-03T01:32:00.056287Z","shell.execute_reply.started":"2023-02-03T01:32:00.041599Z","shell.execute_reply":"2023-02-03T01:32:00.055210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_batches = train_meta.batch_id.unique()\ntrain_meta_baseline = train_meta[train_meta.batch_id.isin(train_batches[:1])]\ntrain_meta_baseline = train_meta_baseline.iloc[:5000]\nprint(len(train_meta_baseline))\ntrain_meta_baseline.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:32:00.424983Z","iopub.execute_input":"2023-02-03T01:32:00.426443Z","iopub.status.idle":"2023-02-03T01:32:01.630536Z","shell.execute_reply.started":"2023-02-03T01:32:00.426370Z","shell.execute_reply":"2023-02-03T01:32:01.629317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cartesian_to_spherical(x, y, z):\n    r = math.sqrt(x**2 + y**2 + z**2)\n    zenith = math.acos(z / r)\n    azimuth = math.atan2(y, x)\n    return azimuth, zenith\n\ndef adjust_spherical(azimuth, zenith):\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    elif zenith < 0:\n        zenith += math.pi\n    return azimuth, zenith\n\ntrain_meta_baseline[['azimuth_pred', 'zenith_pred']] = None\nprev_batch_id = -1\n\nprogress_bar = tqdm(range(len(train_meta_baseline)))\nfor i in progress_bar:\n    # each meta row corresponds to an event\n    meta_row = train_meta_baseline.iloc[i]\n    batch_id = int(meta_row['batch_id'])\n    first_pulse_index = int(meta_row['first_pulse_index'])\n    last_pulse_index = int(meta_row['last_pulse_index'])\n    \n    if batch_id != prev_batch_id:\n        train_batch = pd.read_parquet(f\"/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_{str(meta_row.batch_id)}.parquet\")\n        prev_batch_id = batch_id\n    event_data = train_batch.iloc[first_pulse_index:last_pulse_index + 1]\n    \n    # remove auxiliary pulses if enough non-auxiliary\n    non_auxiliary_pulses = event_data.loc[event_data['auxiliary'] == False]\n    if len(non_auxiliary_pulses) > 10:\n        event_data = non_auxiliary_pulses\n\n    event_data_with_pos = pd.merge(event_data, sensor_geometry, on='sensor_id', how='left')\n    event_points = event_data_with_pos[['x', 'y', 'z']].to_numpy()\n    event_points_centered = event_points - np.mean(event_points, axis=0)\n\n    # predicting\n    pca = PCA(n_components=1).fit(event_points_centered)\n    vector = -1 * pca.components_[0]\n    \n    azimuth, zenith = adjust_spherical(*cartesian_to_spherical(*vector))\n    \n    # setting predictions\n    train_meta_baseline.loc[i, \"azimuth_pred\"] = azimuth\n    train_meta_baseline.loc[i, \"zenith_pred\"] = zenith\n\ntrain_meta_baseline.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:32:01.633063Z","iopub.execute_input":"2023-02-03T01:32:01.633675Z","iopub.status.idle":"2023-02-03T01:32:30.002407Z","shell.execute_reply.started":"2023-02-03T01:32:01.633632Z","shell.execute_reply":"2023-02-03T01:32:30.001619Z"},"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-02-03T01:32:30.003567Z","iopub.execute_input":"2023-02-03T01:32:30.004022Z","iopub.status.idle":"2023-02-03T01:32:30.015434Z","shell.execute_reply.started":"2023-02-03T01:32:30.003996Z","shell.execute_reply":"2023-02-03T01:32:30.014007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Why is the score so high (should be like 1.2x)","metadata":{"execution":{"iopub.status.busy":"2023-02-03T01:31:28.540010Z","iopub.status.idle":"2023-02-03T01:31:28.540376Z","shell.execute_reply.started":"2023-02-03T01:31:28.540218Z","shell.execute_reply":"2023-02-03T01:31:28.540234Z"},"trusted":true},"execution_count":null,"outputs":[]}]}