{"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":"### Intro","metadata":{}},{"cell_type":"markdown","source":"In this notebook we try to predict a neutrino particle’s direction by using PCA. \\\nPrincipal Component Analysis (PCA) projects the data from a higher dimensional space into a lower dimensional space by maximizing the variance of each component **via lenear combinations of the original space**, getting rid of linear correlations and defining a new orthogonal basis known as principal components. So we will use the first principal component to reconstruct the three-dimentional neutrino's track by finding the best lenear combination of the original space of features (the x, y, z coordinates of triggered sensors).","metadata":{}},{"cell_type":"code","source":"import os\nimport 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\n\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:13.643433Z","iopub.execute_input":"2023-03-04T08:26:13.644010Z","iopub.status.idle":"2023-03-04T08:26:15.446551Z","shell.execute_reply.started":"2023-03-04T08:26:13.643898Z","shell.execute_reply":"2023-03-04T08:26:15.445223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH_DATASETS = '/kaggle/input/icecube-neutrinos-in-deep-ice/'","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:15.448815Z","iopub.execute_input":"2023-03-04T08:26:15.449677Z","iopub.status.idle":"2023-03-04T08:26:15.455668Z","shell.execute_reply.started":"2023-03-04T08:26:15.449626Z","shell.execute_reply":"2023-03-04T08:26:15.454305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-04T08:26:15.457559Z","iopub.execute_input":"2023-03-04T08:26:15.457977Z","iopub.status.idle":"2023-03-04T08:26:15.469946Z","shell.execute_reply.started":"2023-03-04T08:26:15.457939Z","shell.execute_reply":"2023-03-04T08:26:15.468670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Train data loading","metadata":{}},{"cell_type":"code","source":"train_meta = pd.read_parquet(os.path.join(PATH_DATASETS,'train_meta.parquet'))\ntrain_meta_batch_125 = train_meta[train_meta.batch_id == 125]\ntrain_meta_batch_125.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:15.472792Z","iopub.execute_input":"2023-03-04T08:26:15.473273Z","iopub.status.idle":"2023-03-04T08:26:54.595294Z","shell.execute_reply.started":"2023-03-04T08:26:15.473236Z","shell.execute_reply":"2023-03-04T08:26:54.594107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry = pd.read_csv(os.path.join(PATH_DATASETS, 'sensor_geometry.csv'))\nsensor_geometry.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:54.597280Z","iopub.execute_input":"2023-03-04T08:26:54.598185Z","iopub.status.idle":"2023-03-04T08:26:54.624888Z","shell.execute_reply.started":"2023-03-04T08:26:54.598133Z","shell.execute_reply":"2023-03-04T08:26:54.623596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_batch_125 = pd.read_parquet( \\\n    os.path.join(PATH_DATASETS, 'train/batch_125.parquet')).reset_index()\ntrain_batch_125.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:54.626191Z","iopub.execute_input":"2023-03-04T08:26:54.626521Z","iopub.status.idle":"2023-03-04T08:26:58.393062Z","shell.execute_reply.started":"2023-03-04T08:26:54.626492Z","shell.execute_reply":"2023-03-04T08:26:58.391826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"auxiliary_ratio_df = train_batch_125.pivot_table(index='event_id', values='auxiliary', aggfunc='mean') \\\n                                    .rename(columns={'auxiliary': 'auxiliary_ratio'})\ntrain_meta_batch_125 = train_meta_batch_125.merge(auxiliary_ratio_df, on='event_id')","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:26:58.395173Z","iopub.execute_input":"2023-03-04T08:26:58.396096Z","iopub.status.idle":"2023-03-04T08:27:00.676275Z","shell.execute_reply.started":"2023-03-04T08:26:58.396040Z","shell.execute_reply":"2023-03-04T08:27:00.675176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_batch_125['auxiliary_ratio'].describe()","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.677545Z","iopub.execute_input":"2023-03-04T08:27:00.677859Z","iopub.status.idle":"2023-03-04T08:27:00.705109Z","shell.execute_reply.started":"2023-03-04T08:27:00.677832Z","shell.execute_reply":"2023-03-04T08:27:00.703813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PCA fitting for \"noisy\" event","metadata":{}},{"cell_type":"code","source":"train_meta_batch_125[train_meta_batch_125.auxiliary_ratio > 0.95].sample(1, random_state=1)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.706803Z","iopub.execute_input":"2023-03-04T08:27:00.707558Z","iopub.status.idle":"2023-03-04T08:27:00.731981Z","shell.execute_reply.started":"2023-03-04T08:27:00.707507Z","shell.execute_reply":"2023-03-04T08:27:00.730712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_id_1 = 404418940","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.737232Z","iopub.execute_input":"2023-03-04T08:27:00.738001Z","iopub.status.idle":"2023-03-04T08:27:00.743262Z","shell.execute_reply.started":"2023-03-04T08:27:00.737951Z","shell.execute_reply":"2023-03-04T08:27:00.742145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses = train_batch_125[train_batch_125.event_id == event_id_1].copy()\nevent_pulses = event_pulses.reset_index()\nprint(event_pulses.shape)\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.744940Z","iopub.execute_input":"2023-03-04T08:27:00.745331Z","iopub.status.idle":"2023-03-04T08:27:00.874766Z","shell.execute_reply.started":"2023-03-04T08:27:00.745289Z","shell.execute_reply":"2023-03-04T08:27:00.873539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Adding zero-centered coordinates of each sensor of the event","metadata":{}},{"cell_type":"code","source":"event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\n\nevent_pulses.x = event_pulses.x - event_pulses.x.mean()\nevent_pulses.y = event_pulses.y - event_pulses.y.mean()\nevent_pulses.z = event_pulses.z - event_pulses.z.mean()\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.876070Z","iopub.execute_input":"2023-03-04T08:27:00.876407Z","iopub.status.idle":"2023-03-04T08:27:00.901138Z","shell.execute_reply.started":"2023-03-04T08:27:00.876377Z","shell.execute_reply":"2023-03-04T08:27:00.899782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Calculating the true vector of the event","metadata":{}},{"cell_type":"code","source":"zenith_true = train_meta_batch_125.loc[train_meta_batch_125.event_id == event_id_1, 'zenith'].values\nazimuth_true = train_meta_batch_125.loc[train_meta_batch_125.event_id == event_id_1, 'azimuth'].values\n\nvector = [np.cos(azimuth_true) * np.sin(zenith_true),\n          np.sin(azimuth_true) * np.sin(zenith_true),\n          np.cos(azimuth_true)]\n\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})\nvector_df","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.902630Z","iopub.execute_input":"2023-03-04T08:27:00.903726Z","iopub.status.idle":"2023-03-04T08:27:00.924314Z","shell.execute_reply.started":"2023-03-04T08:27:00.903685Z","shell.execute_reply":"2023-03-04T08:27:00.922981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Event visualizition","metadata":{}},{"cell_type":"markdown","source":"Now let's see if the first principal component fits the neutrino's diraction.","metadata":{}},{"cell_type":"code","source":"#Creating a PCA object and fitting it to pulses with highest charge\npca = PCA(n_components=1).fit(event_pulses.loc[event_pulses.charge > event_pulses.charge.mean()][['x', 'y', 'z']])\n\n# Calculating the first principal component vector\nvector_pca = pca.components_[0]\nvector_base = np.array([-500, 500])\nx_pca = vector_base*vector_pca[0]\ny_pca = vector_base*vector_pca[1]\nz_pca = vector_base*vector_pca[2]\nvector_pca_df = pd.DataFrame({'x_pca': x_pca, 'y_pca': y_pca, 'z_pca': z_pca})\n\n# Plotting auxiliary pulses\nfig_auxiliary = px.scatter_3d(event_pulses.loc[event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\n# Plotting non-auxiliary pulses\nfig_non_auxiliary = px.scatter_3d(event_pulses.loc[~event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\n\n# Plotting the true vector\nfig_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\", color_discrete_sequence=['black'])\n\n# Plotting the first principal component vector\nfig_pca_line = px.line_3d(vector_pca_df, x=\"x_pca\", y=\"y_pca\", z=\"z_pca\", color_discrete_sequence=['green'])\n\nfig_auxiliary.update_traces(marker_size=event_pulses.charge*10)\nfig_non_auxiliary.update_traces(marker_size=event_pulses.charge*10)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data + fig_pca_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:00.926245Z","iopub.execute_input":"2023-03-04T08:27:00.926848Z","iopub.status.idle":"2023-03-04T08:27:02.403138Z","shell.execute_reply.started":"2023-03-04T08:27:00.926800Z","shell.execute_reply":"2023-03-04T08:27:02.402160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check out the amount of variance explained by first principal component.","metadata":{}},{"cell_type":"code","source":"for i, component in enumerate(pca.components_):\n    print(\"{} component: {}% of initial variance\".format(i + 1, \n          round(100 * pca.explained_variance_ratio_[i], 2)))","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.404864Z","iopub.execute_input":"2023-03-04T08:27:02.405372Z","iopub.status.idle":"2023-03-04T08:27:02.413073Z","shell.execute_reply.started":"2023-03-04T08:27:02.405326Z","shell.execute_reply":"2023-03-04T08:27:02.411869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's calculate the error.","metadata":{}},{"cell_type":"code","source":"zenith = np.arccos(vector_pca[2])\nazimuth = np.arctan2(vector_pca[1], vector_pca[0])\nif azimuth < 0:\n    azimuth = 2*np.pi + azimuth\n\nevent_error = angular_dist_score(azimuth_true, zenith_true, azimuth, zenith)\nprint(f'The value of error: {event_error:.2f}')","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.414755Z","iopub.execute_input":"2023-03-04T08:27:02.415275Z","iopub.status.idle":"2023-03-04T08:27:02.427898Z","shell.execute_reply.started":"2023-03-04T08:27:02.415231Z","shell.execute_reply":"2023-03-04T08:27:02.426668Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PCA fitting for non-auxiliary pulses","metadata":{}},{"cell_type":"markdown","source":"Let's try to use PCA fitting again, but, this time, let's observe the event with larger number of non-auxiliary pulses.","metadata":{}},{"cell_type":"code","source":"train_meta_batch_125[train_meta_batch_125.event_id == 403533241]","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.429479Z","iopub.execute_input":"2023-03-04T08:27:02.429897Z","iopub.status.idle":"2023-03-04T08:27:02.448560Z","shell.execute_reply.started":"2023-03-04T08:27:02.429852Z","shell.execute_reply":"2023-03-04T08:27:02.447300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_id_2 = 403533241","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.449861Z","iopub.execute_input":"2023-03-04T08:27:02.450230Z","iopub.status.idle":"2023-03-04T08:27:02.455419Z","shell.execute_reply.started":"2023-03-04T08:27:02.450173Z","shell.execute_reply":"2023-03-04T08:27:02.454224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses = train_batch_125[train_batch_125.event_id == event_id_2].copy()\nevent_pulses = event_pulses.reset_index()\nprint(event_pulses.shape)\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.456738Z","iopub.execute_input":"2023-03-04T08:27:02.457663Z","iopub.status.idle":"2023-03-04T08:27:02.515070Z","shell.execute_reply.started":"2023-03-04T08:27:02.457627Z","shell.execute_reply":"2023-03-04T08:27:02.513944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Adding zero-centered coordinates of each sensor of the event","metadata":{"execution":{"iopub.status.busy":"2023-02-25T14:14:48.907766Z","iopub.execute_input":"2023-02-25T14:14:48.908334Z","iopub.status.idle":"2023-02-25T14:14:48.914747Z","shell.execute_reply.started":"2023-02-25T14:14:48.908286Z","shell.execute_reply":"2023-02-25T14:14:48.913296Z"}}},{"cell_type":"code","source":"event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\nevent_pulses.x = event_pulses.x - event_pulses.x.mean()\nevent_pulses.y = event_pulses.y - event_pulses.y.mean()\nevent_pulses.z = event_pulses.z - event_pulses.z.mean()\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.516564Z","iopub.execute_input":"2023-03-04T08:27:02.516942Z","iopub.status.idle":"2023-03-04T08:27:02.544159Z","shell.execute_reply.started":"2023-03-04T08:27:02.516904Z","shell.execute_reply":"2023-03-04T08:27:02.542845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Calculating the true vector of the event","metadata":{}},{"cell_type":"code","source":"zenith_true = train_meta_batch_125.loc[train_meta_batch_125.event_id == event_id_2, 'zenith'].values\nazimuth_true = train_meta_batch_125.loc[train_meta_batch_125.event_id == event_id_2, 'azimuth'].values\n\nvector = [np.cos(azimuth_true) * np.sin(zenith_true),\n          np.sin(azimuth_true) * np.sin(zenith_true),\n          np.cos(azimuth_true)]\n\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})\nvector_df","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.546139Z","iopub.execute_input":"2023-03-04T08:27:02.546603Z","iopub.status.idle":"2023-03-04T08:27:02.566284Z","shell.execute_reply.started":"2023-03-04T08:27:02.546557Z","shell.execute_reply":"2023-03-04T08:27:02.565080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Event visualizition\n","metadata":{}},{"cell_type":"code","source":"#Creating a PCA object and fitting it to only non-auxiliary pulses\npca = PCA(n_components=1).fit(event_pulses.loc[~event_pulses.auxiliary][['x', 'y', 'z']])\n\n# Calculating the first principal component vector\nvector_pca = pca.components_[0]\nvector_base = np.array([-500, 500])\nx_pca = vector_base*vector_pca[0]\ny_pca = vector_base*vector_pca[1]\nz_pca = vector_base*vector_pca[2]\nvector_pca_df = pd.DataFrame({'x_pca': x_pca, 'y_pca': y_pca, 'z_pca': z_pca})\n\n# Plotting auxiliary pulses\nfig_auxiliary = px.scatter_3d(event_pulses.loc[event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\n# Plotting non-auxiliary pulses\nfig_non_auxiliary = px.scatter_3d(event_pulses.loc[~event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\n\n# Plotting the true vector\nfig_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\", color_discrete_sequence=['black'])\n\n# Plotting the first principal component vector\nfig_pca_line = px.line_3d(vector_pca_df, x=\"x_pca\", y=\"y_pca\", z=\"z_pca\", color_discrete_sequence=['green'])\n\nfig_auxiliary.update_traces(marker_size=event_pulses.charge*10)\nfig_non_auxiliary.update_traces(marker_size=event_pulses.charge*10)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data + fig_pca_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.567798Z","iopub.execute_input":"2023-03-04T08:27:02.568161Z","iopub.status.idle":"2023-03-04T08:27:02.794904Z","shell.execute_reply.started":"2023-03-04T08:27:02.568129Z","shell.execute_reply":"2023-03-04T08:27:02.794003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check out the amount of variance explained by first principal component.","metadata":{}},{"cell_type":"code","source":"for i, component in enumerate(pca.components_):\n    print(\"{} component: {}% of initial variance\".format(i + 1, \n          round(100 * pca.explained_variance_ratio_[i], 2)))","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.796028Z","iopub.execute_input":"2023-03-04T08:27:02.796800Z","iopub.status.idle":"2023-03-04T08:27:02.802973Z","shell.execute_reply.started":"2023-03-04T08:27:02.796763Z","shell.execute_reply":"2023-03-04T08:27:02.801862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's calculate the error.","metadata":{}},{"cell_type":"code","source":"zenith = np.arccos(vector_pca[2])\nazimuth = np.arctan2(vector_pca[1], vector_pca[0])\nif azimuth < 0:\n    azimuth = 2*np.pi + azimuth\n\nevent_error = angular_dist_score(azimuth_true, zenith_true, azimuth, zenith)\nprint(f'The value of error: {event_error:.2f}')","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.804522Z","iopub.execute_input":"2023-03-04T08:27:02.805284Z","iopub.status.idle":"2023-03-04T08:27:02.816658Z","shell.execute_reply.started":"2023-03-04T08:27:02.805235Z","shell.execute_reply":"2023-03-04T08:27:02.815528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PCA fitting with ranked pulses","metadata":{}},{"cell_type":"code","source":"train_meta_sample = train_meta_batch_125[:10000].copy()\ntrain_meta_sample.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.817765Z","iopub.execute_input":"2023-03-04T08:27:02.818521Z","iopub.status.idle":"2023-03-04T08:27:02.841006Z","shell.execute_reply.started":"2023-03-04T08:27:02.818484Z","shell.execute_reply":"2023-03-04T08:27:02.839628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sample = train_batch_125.loc[:train_meta_sample.iloc[-1].last_pulse_index.astype(int)].reset_index()\ntrain_sample[['x', 'y', 'z']] = sensor_geometry.loc[train_sample.sensor_id].reset_index()[['x', 'y', 'z']]\n\ntrain_sample['x'] = train_sample['x'] - train_sample.groupby('event_id').x.transform('mean')\ntrain_sample['y'] = train_sample['y'] - train_sample.groupby('event_id').y.transform('mean')\ntrain_sample['z'] = train_sample['z'] - train_sample.groupby('event_id').z.transform('mean')\n\nprint(len(train_sample))\ntrain_sample.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:02.842202Z","iopub.execute_input":"2023-03-04T08:27:02.842549Z","iopub.status.idle":"2023-03-04T08:27:03.326745Z","shell.execute_reply.started":"2023-03-04T08:27:02.842518Z","shell.execute_reply":"2023-03-04T08:27:03.325553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's rank our pulses, using this ranking system (<a href=\"https://www.kaggle.com/code/seungmoklee/3-lstms-with-data-picking-and-shifting#Data-I/O-(for-CPU)-&-Normalization\">Notebook link</a>):\n- Rank 3 will get non-auxiliary pulses within the valid time window;\n- Rank 2 - non-auxiliary pulses out of the valid time window;\n- Rank 1 - auxiliary pulses within the valid time window;\n- Rank 0 - auxiliary pulses out of the valid time window.\n\nThe information about the valid time window: \"as neutrinos travel with the speed of light, it must take less than 6200 ns to transverse the detector. So the pulses further than 6200 ns may not be in our interest\".","metadata":{}},{"cell_type":"markdown","source":"Below you can find the code that calculates the value of this threshold transit time for neutrinos - 6200 ns.","metadata":{}},{"cell_type":"code","source":"#the speed of light\nc_const = 0.299792458\n\nx_min, x_max = sensor_geometry.x.agg([min, max]).values\ny_min, y_max = sensor_geometry.y.agg([min, max]).values\nz_min, z_max = sensor_geometry.z.agg([min, max]).values\n\ndetector_length = np.sqrt((x_max - x_min)**2 + (y_max - y_min)**2 + (z_max - z_min)**2)\ntime_valid_length = detector_length / c_const\ntime_valid_length","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:03.328276Z","iopub.execute_input":"2023-03-04T08:27:03.329215Z","iopub.status.idle":"2023-03-04T08:27:03.342468Z","shell.execute_reply.started":"2023-03-04T08:27:03.329167Z","shell.execute_reply":"2023-03-04T08:27:03.341073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_ranked_pulses(event_pulses):\n    event_pulses[\"time_from_start\"] = event_pulses.time.values - event_pulses.time.min()\n    time_peak_charge = event_pulses.loc[event_pulses['charge'].argmax(), 'time_from_start']\n    time_valid_min = time_peak_charge - time_valid_length\n    time_valid_max = time_peak_charge + time_valid_length\n    time_valid_window = ((event_pulses['time_from_start'] > time_valid_min) * \n                  (event_pulses['time_from_start'] < time_valid_max))\n    event_pulses['pulse_rank'] = 2 * (1 - event_pulses['auxiliary'].values) + (time_valid_window)\n    ranked_pulses = event_pulses.drop('time_from_start', axis=1)\n    return ranked_pulses","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:03.347305Z","iopub.execute_input":"2023-03-04T08:27:03.347916Z","iopub.status.idle":"2023-03-04T08:27:03.355417Z","shell.execute_reply.started":"2023-03-04T08:27:03.347857Z","shell.execute_reply":"2023-03-04T08:27:03.354060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_angles(event_pulses):\n    event_pulses = (get_ranked_pulses(event_pulses)).query('pulse_rank == 3')\n    try:\n        pca = PCA(n_components=1).fit(event_pulses[['x', 'y', 'z']])\n        vector = pca.components_[0]\n\n        # zero-centering of the coordinates\n        event_pulses.x = event_pulses.x - event_pulses.x.mean()\n        event_pulses.y = event_pulses.y - event_pulses.y.mean()\n        event_pulses.z = event_pulses.z - event_pulses.z.mean()\n\n        # calculating a distance of a point from the first principal component, that goes through the origin\n        event_pulses['distance'] = (event_pulses.x * vector[0] + event_pulses.y * vector[1] +\n                                              event_pulses.z * vector[2]) / np.linalg.norm(vector)\n\n        # Flip the vector direction if it points away from the neutrino origin\n        if (event_pulses.loc[event_pulses.distance > 0].time.mean() >\n            event_pulses.loc[event_pulses.distance < 0].time.mean()):\n            vector = -1 * vector\n            \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\n    except:\n        zenith, azimuth = 0., 0.\n        \n    return zenith, azimuth","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:03.356934Z","iopub.execute_input":"2023-03-04T08:27:03.357259Z","iopub.status.idle":"2023-03-04T08:27:03.371715Z","shell.execute_reply.started":"2023-03-04T08:27:03.357221Z","shell.execute_reply":"2023-03-04T08:27:03.370023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ntrain_meta_sample[['azimuth_pred', 'zenith_pred']] = None\nfor i in range(len(train_meta_sample)):\n    row = train_meta_sample.iloc[i]\n    event_pulses = train_sample.iloc[ \\\n        row.first_pulse_index.astype(int):row.last_pulse_index.astype(int)+1].reset_index().copy()\n    zenith, azimuth = get_angles(event_pulses)\n    train_meta_sample.loc[i, 'zenith_pred'] = zenith\n    train_meta_sample.loc[i, 'azimuth_pred'] = azimuth","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:27:03.373664Z","iopub.execute_input":"2023-03-04T08:27:03.374058Z","iopub.status.idle":"2023-03-04T08:29:40.592436Z","shell.execute_reply.started":"2023-03-04T08:27:03.374020Z","shell.execute_reply":"2023-03-04T08:29:40.590108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_error = angular_dist_score(train_meta_sample.azimuth.to_list(), train_meta_sample.zenith.to_list(),\n                   train_meta_sample.azimuth_pred.astype(float).to_list(),\n                   train_meta_sample.zenith_pred.astype(float).to_list())\n\nprint(f'The value of error: {sample_error:.3f}')","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:29:40.595701Z","iopub.execute_input":"2023-03-04T08:29:40.596526Z","iopub.status.idle":"2023-03-04T08:29:40.623438Z","shell.execute_reply.started":"2023-03-04T08:29:40.596460Z","shell.execute_reply":"2023-03-04T08:29:40.621931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test data loading","metadata":{}},{"cell_type":"code","source":"test_meta = pd.read_parquet(os.path.join(PATH_DATASETS,'test_meta.parquet'))\ntest_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:29:40.625427Z","iopub.execute_input":"2023-03-04T08:29:40.625792Z","iopub.status.idle":"2023-03-04T08:29:40.656245Z","shell.execute_reply.started":"2023-03-04T08:29:40.625761Z","shell.execute_reply":"2023-03-04T08:29:40.654823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def load_test_batch(batch_id):\n    test_batch = pd.read_parquet(os.path.join(PATH_DATASETS, f'test/batch_{batch_id}.parquet')) \\\n                                .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-03-04T08:29:40.658245Z","iopub.execute_input":"2023-03-04T08:29:40.658994Z","iopub.status.idle":"2023-03-04T08:29:40.669762Z","shell.execute_reply.started":"2023-03-04T08:29:40.658944Z","shell.execute_reply":"2023-03-04T08:29:40.668827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for batch in test_meta.batch_id.unique():\n    test_batch = load_test_batch(batch)\n    test_meta_batch = test_meta.loc[test_meta.batch_id == batch]\n\n    pulses_info = []\n    for i in range(len(test_meta_batch)):\n        row = test_meta_batch.iloc[i]\n        pulses_info.append(test_batch.iloc[ \\\n            row.first_pulse_index.astype(int):row.last_pulse_index.astype(int)+1].reset_index().copy())\n        \n    results = []\n    for result in map(get_angles, pulses_info):\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    \ntest_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-04T08:29:40.672065Z","iopub.execute_input":"2023-03-04T08:29:40.672586Z","iopub.status.idle":"2023-03-04T08:29:40.756032Z","shell.execute_reply.started":"2023-03-04T08:29:40.672539Z","shell.execute_reply":"2023-03-04T08:29:40.754710Z"},"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-03-04T08:29:40.758034Z","iopub.execute_input":"2023-03-04T08:29:40.758516Z","iopub.status.idle":"2023-03-04T08:29:40.779694Z","shell.execute_reply.started":"2023-03-04T08:29:40.758468Z","shell.execute_reply":"2023-03-04T08:29:40.778822Z"},"trusted":true},"execution_count":null,"outputs":[]}]}