{"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":"**Introduction**","metadata":{}},{"cell_type":"markdown","source":"The basic idea of this notebook is to focus on a comparison between Standard and Weighted PCA algorithm perfomance for predicting the neutrino particle’s diraction and give some demonstration of how these two algorithms work.\n\n**Standard PCA** is a well-known technique initially designed to reduce the dimensionality of an original dataset by finding the orthogonal projection of the data into a lower-dimensional linear space such that the variance of the projected data is maximised. So the new orthogonal basis known as principal components aim to explain at best the whole variance of a given data set. Nevetheless, the main limitation of Standard PCA comes from the fact that this algorithm cannot handle datasets with noisy entries. It means that the classical PCA implementations make no difference between variance coming from a genuine underlying signal and variance coming from noise.\n\nOn the other hand, the basic idea of **the Weighted PCA** algorithm is to focus on the maximization of the weighted variance explained by each principal component through the diagonalization of the associated weighted covariance matrix. The resulting principal components will then be those that are the most significant in identifying pattern within the dataset even if their linear combination is not necessarily the best at explaining the total dataset variance.\n\nLet's take a look at these methods in action and determine how these two algorithms differ in their ability to predict which direction neutrinos came from, and whether we can choose only one of them with better predictive power.","metadata":{}},{"cell_type":"code","source":"!mkdir -p /tmp/pip/cache/\n!cp /kaggle/input/wpca-01-py3-whl/wpca-0.1-py3-none-any.whl /tmp/pip/cache/\n!pip install --no-index --find-links /tmp/pip/cache/ wpca","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:12:31.518158Z","iopub.execute_input":"2023-03-14T13:12:31.518609Z","iopub.status.idle":"2023-03-14T13:12:47.082275Z","shell.execute_reply.started":"2023-03-14T13:12:31.518567Z","shell.execute_reply":"2023-03-14T13:12:47.080664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport plotly.express as px\nfrom plotly.subplots import make_subplots\nfrom sklearn.decomposition import PCA\nfrom wpca import WPCA\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nimport gc\n\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:12:47.086034Z","iopub.execute_input":"2023-03-14T13:12:47.086566Z","iopub.status.idle":"2023-03-14T13:12:49.828144Z","shell.execute_reply.started":"2023-03-14T13:12:47.086509Z","shell.execute_reply":"2023-03-14T13:12:49.826522Z"},"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-14T13:12:49.830066Z","iopub.execute_input":"2023-03-14T13:12:49.830461Z","iopub.status.idle":"2023-03-14T13:12:49.836085Z","shell.execute_reply.started":"2023-03-14T13:12:49.830421Z","shell.execute_reply":"2023-03-14T13:12:49.834516Z"},"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-14T13:12:49.839481Z","iopub.execute_input":"2023-03-14T13:12:49.840032Z","iopub.status.idle":"2023-03-14T13:12:49.851516Z","shell.execute_reply.started":"2023-03-14T13:12:49.839992Z","shell.execute_reply":"2023-03-14T13:12:49.850360Z"},"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-14T13:12:49.853282Z","iopub.execute_input":"2023-03-14T13:12:49.854155Z","iopub.status.idle":"2023-03-14T13:13:31.208146Z","shell.execute_reply.started":"2023-03-14T13:12:49.854108Z","shell.execute_reply":"2023-03-14T13:13:31.206518Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_meta","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:31.209635Z","iopub.execute_input":"2023-03-14T13:13:31.210948Z","iopub.status.idle":"2023-03-14T13:13:31.222630Z","shell.execute_reply.started":"2023-03-14T13:13:31.210891Z","shell.execute_reply":"2023-03-14T13:13:31.221279Z"},"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-14T13:13:31.224289Z","iopub.execute_input":"2023-03-14T13:13:31.225189Z","iopub.status.idle":"2023-03-14T13:13:31.255134Z","shell.execute_reply.started":"2023-03-14T13:13:31.225142Z","shell.execute_reply":"2023-03-14T13:13:31.253953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The code below will add to sensor_geometry dataframe the column with value 1 if sensor is located on one of the DeepCore strings, otherwise 0.","metadata":{}},{"cell_type":"code","source":"sensor_geometry['line_id'] = sensor_geometry.sensor_id // 60 + 1\nsensor_geometry['core'] = (sensor_geometry.line_id > 78).astype(np.int8)\nsensor_geometry.drop('line_id', axis=1, inplace=True)\nsensor_geometry.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:31.256524Z","iopub.execute_input":"2023-03-14T13:13:31.256910Z","iopub.status.idle":"2023-03-14T13:13:31.278648Z","shell.execute_reply.started":"2023-03-14T13:13:31.256873Z","shell.execute_reply":"2023-03-14T13:13:31.277158Z"},"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-14T13:13:31.280568Z","iopub.execute_input":"2023-03-14T13:13:31.281558Z","iopub.status.idle":"2023-03-14T13:13:40.906412Z","shell.execute_reply.started":"2023-03-14T13:13:31.281503Z","shell.execute_reply":"2023-03-14T13:13:40.904726Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The Standard PCA fitting","metadata":{}},{"cell_type":"code","source":"event_id_1 = 403533054","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:40.912952Z","iopub.execute_input":"2023-03-14T13:13:40.913378Z","iopub.status.idle":"2023-03-14T13:13:40.922097Z","shell.execute_reply.started":"2023-03-14T13:13:40.913341Z","shell.execute_reply":"2023-03-14T13:13:40.920065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses = train_batch_125[train_batch_125.event_id == event_id_1].reset_index(drop=True).copy()\nprint(event_pulses.shape)\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:40.923650Z","iopub.execute_input":"2023-03-14T13:13:40.924044Z","iopub.status.idle":"2023-03-14T13:13:42.093836Z","shell.execute_reply.started":"2023-03-14T13:13:40.923999Z","shell.execute_reply":"2023-03-14T13:13:42.086475Z"},"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[ \\\n                                            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-14T13:13:42.100772Z","iopub.execute_input":"2023-03-14T13:13:42.102487Z","iopub.status.idle":"2023-03-14T13:13:42.210128Z","shell.execute_reply.started":"2023-03-14T13:13:42.102389Z","shell.execute_reply":"2023-03-14T13:13:42.204459Z"},"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 \\\n                        .loc[train_meta_batch_125.event_id == event_id_1, 'zenith'].values\nazimuth_true = train_meta_batch_125 \\\n                        .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-14T13:13:42.216257Z","iopub.execute_input":"2023-03-14T13:13:42.220946Z","iopub.status.idle":"2023-03-14T13:13:42.311275Z","shell.execute_reply.started":"2023-03-14T13:13:42.220508Z","shell.execute_reply":"2023-03-14T13:13:42.306112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The Standard PCA fitting for all pulses of the event","metadata":{}},{"cell_type":"code","source":"#Creating a PCA object and fitting it to all pulses\npca = PCA(n_components=1).fit(event_pulses[['x', 'y', 'z']])\n\n# Calculating the first principal component vector\nvector_pca_all = pca.components_[0]\n\n#Calculating a distance of a point from the first principal component,that goes through the origin\nevent_pulses['distance'] = (event_pulses.x * vector_pca_all[0] +\n                            event_pulses.y * vector_pca_all[1] +\n                            event_pulses.z * vector_pca_all[2]) / np.linalg.norm(vector_pca_all)\n\n#Flip the vector direction if it points away from the neutrino origin\nif (event_pulses.loc[event_pulses.distance > 0].time.mean() >\n    event_pulses.loc[event_pulses.distance < 0].time.mean()):\n    vector_pca_all = -1 * np.array(vector_pca_all)\n\nvector_pca_all = np.clip(vector_pca_all, -1, 1)\n\n# Calculating the angles\nzen_pca_all = np.arccos(vector_pca_all[2])\naz_pca_all = np.arctan2(vector_pca_all[1], vector_pca_all[0])\nif az_pca_all < 0:\n    az_pca_all = 2*np.pi + az_pca_all","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.316381Z","iopub.execute_input":"2023-03-14T13:13:42.317509Z","iopub.status.idle":"2023-03-14T13:13:42.406565Z","shell.execute_reply.started":"2023-03-14T13:13:42.317310Z","shell.execute_reply":"2023-03-14T13:13:42.402098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating the value of error\npca_all_error = angular_dist_score(azimuth_true, zenith_true, az_pca_all, zen_pca_all)\nprint(f'The value of error: {pca_all_error}')","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.412491Z","iopub.execute_input":"2023-03-14T13:13:42.415676Z","iopub.status.idle":"2023-03-14T13:13:42.443270Z","shell.execute_reply.started":"2023-03-14T13:13:42.415580Z","shell.execute_reply":"2023-03-14T13:13:42.437388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating predicted vector of the event\nx_pca_all = vector_base * vector_pca_all[0]\ny_pca_all = vector_base * vector_pca_all[1]\nz_pca_all = vector_base * vector_pca_all[2]\nvector_pca_all_df = pd.DataFrame({'x_pca': x_pca_all, 'y_pca': y_pca_all, 'z_pca': z_pca_all})\nvector_pca_all_df","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.450442Z","iopub.execute_input":"2023-03-14T13:13:42.452213Z","iopub.status.idle":"2023-03-14T13:13:42.500904Z","shell.execute_reply.started":"2023-03-14T13:13:42.452009Z","shell.execute_reply":"2023-03-14T13:13:42.496272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The Standard PCA fitting for non-auxiliary pulses of the event","metadata":{}},{"cell_type":"code","source":"#Creating a PCA object and fitting it to non-aux 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_non_aux = pca.components_[0]\n\n#Calculating a distance of a point from the first principal component,that goes through the origin\nevent_pulses['distance'] = (event_pulses.x * vector_pca_non_aux[0] +\n                            event_pulses.y * vector_pca_non_aux[1] +\n                            event_pulses.z * vector_pca_non_aux[2]) / np.linalg.norm(vector_pca_non_aux)\n\n#Flip the vector direction if it points away from the neutrino origin\nif (event_pulses.loc[event_pulses.distance > 0].time.mean() >\n    event_pulses.loc[event_pulses.distance < 0].time.mean()):\n    vector_pca_non_aux = -1 * np.array(vector_pca_non_aux)\n\nvector_pca_non_aux = np.clip(vector_pca_non_aux, -1, 1)\n\n# Calculating the angles\nzen_pca_non = np.arccos(vector_pca_non_aux[2])\naz_pca_non = np.arctan2(vector_pca_non_aux[1], vector_pca_non_aux[0])\nif az_pca_non < 0:\n    az_pca_non = 2*np.pi + az_pca_non","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.505843Z","iopub.execute_input":"2023-03-14T13:13:42.509315Z","iopub.status.idle":"2023-03-14T13:13:42.572424Z","shell.execute_reply.started":"2023-03-14T13:13:42.509261Z","shell.execute_reply":"2023-03-14T13:13:42.563814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating the value of error\npca_non_error = angular_dist_score(azimuth_true, zenith_true, az_pca_non, zen_pca_non)\nprint(f'The value of error: {pca_non_error}')","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.581572Z","iopub.execute_input":"2023-03-14T13:13:42.585896Z","iopub.status.idle":"2023-03-14T13:13:42.618226Z","shell.execute_reply.started":"2023-03-14T13:13:42.585666Z","shell.execute_reply":"2023-03-14T13:13:42.608729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating predicted vector of the event\nx_pca_non = vector_base * vector_pca_non_aux[0]\ny_pca_non = vector_base * vector_pca_non_aux[1]\nz_pca_non = vector_base * vector_pca_non_aux[2]\nvector_pca_non_aux_df = pd.DataFrame({'x_pca': x_pca_non, 'y_pca': y_pca_non, 'z_pca': z_pca_non})\nvector_pca_non_aux_df","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.624078Z","iopub.execute_input":"2023-03-14T13:13:42.625397Z","iopub.status.idle":"2023-03-14T13:13:42.682665Z","shell.execute_reply.started":"2023-03-14T13:13:42.625224Z","shell.execute_reply":"2023-03-14T13:13:42.679683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The standard PCA fitting visualization","metadata":{}},{"cell_type":"code","source":"aux_pulses = event_pulses[event_pulses.auxiliary]\nnon_aux_pulses = event_pulses[~event_pulses.auxiliary]\naux_df = aux_pulses[['x','y','z','charge','time']]\nnon_aux_df = non_aux_pulses[['x','y','z','charge','time']]\n\nfig = make_subplots(rows=1, cols=2,\n        specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],\n        subplot_titles=(f'All pulses fitting (error:{pca_all_error:.2f})',\n                        f'Non-auxiliary pulses fitting (error: {pca_non_error:.2f})')\n)\n\n# Plotting auxiliary pulses\nfig.add_trace(\n        go.Scatter3d(\n            x=aux_df.x, y=aux_df.y, z=aux_df.z,\n            opacity=0.75, mode='markers',\n            marker_size=aux_df.charge*10,\n            text=aux_df.charge,\n            marker=dict(color='red'),\n            name='Aux pulses'\n        ),\n        row=1, col=1\n)\n# Plotting non-auxiliary pulses\nfig.add_trace(\n        go.Scatter3d(\n            x=non_aux_df.x, y=non_aux_df.y, z=non_aux_df.z,\n            opacity=0.75, mode='markers',\n            marker_size=non_aux_df.charge*10,\n            text=non_aux_df.charge,\n            marker=dict(color='blue'),\n            name='Non-aux pulses'\n        ),\n        row='all', col='all'\n)\n \n# Plotting the true vector\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_df.x, y=vector_df.y, z=vector_df.z,\n            opacity=0.8, mode='lines', \n            line=dict(color='gray', width=4),\n            name='True vector'\n        ),\n        \n    row='all', col='all'\n)\n\n# Plotting the PC1 vector (fitted to all pulses)\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_pca_all_df.x_pca, y=vector_pca_all_df.y_pca, z=vector_pca_all_df.z_pca,\n            opacity=0.8, mode='lines', \n            line=dict(color='red', width=3),\n            name='Standard PCA (all pulses)'\n        ),      \n    row=1, col=1\n)\n# Plotting the PC1 vector (fitted to non-aux pulses)\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_pca_non_aux_df.x_pca, y=vector_pca_non_aux_df.y_pca, z=vector_pca_non_aux_df.z_pca,\n            opacity=0.8, mode='lines', \n            line=dict(color='purple', width=3),\n            name='Standard PCA (non-aux pulses)'\n        ),      \n    row=1, col=2\n)\n\nfig.update_layout(title=dict(text=\"Standard PCA\",\n                             font_family='Arial', pad_l=300,\n                             font_color='black', font_size=25), legend=dict(font_color='black', font_size=10))\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:42.689173Z","iopub.execute_input":"2023-03-14T13:13:42.692030Z","iopub.status.idle":"2023-03-14T13:13:44.403821Z","shell.execute_reply.started":"2023-03-14T13:13:42.691963Z","shell.execute_reply":"2023-03-14T13:13:44.399303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The Weighted PCA fitting","metadata":{}},{"cell_type":"markdown","source":"#### Assigning ranks to all pulses in the event","metadata":{}},{"cell_type":"markdown","source":"Let's assign ranks to pulses, according to the following 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":{"_kg_hide-input":true}},{"cell_type":"markdown","source":"Calculating the value of this threshold transit time for neutrinos - 6200 ns.","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"c_const = 0.299792458  # speed of light [m/ns]\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-14T13:13:44.414708Z","iopub.execute_input":"2023-03-14T13:13:44.418398Z","iopub.status.idle":"2023-03-14T13:13:44.494044Z","shell.execute_reply.started":"2023-03-14T13:13:44.417890Z","shell.execute_reply":"2023-03-14T13:13:44.485971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Making a function that assigns suitable ranks to all pulses in the event.\n","metadata":{}},{"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 = ((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)\n    event_pulses = event_pulses.drop('time_from_start', axis=1)\n    return event_pulses","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.507263Z","iopub.execute_input":"2023-03-14T13:13:44.511314Z","iopub.status.idle":"2023-03-14T13:13:44.548671Z","shell.execute_reply.started":"2023-03-14T13:13:44.510863Z","shell.execute_reply":"2023-03-14T13:13:44.543448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Applying the function to our event.","metadata":{}},{"cell_type":"code","source":"event_pulses = get_ranked_pulses(event_pulses)\nevent_pulses.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.561357Z","iopub.execute_input":"2023-03-14T13:13:44.564611Z","iopub.status.idle":"2023-03-14T13:13:44.672353Z","shell.execute_reply.started":"2023-03-14T13:13:44.564256Z","shell.execute_reply":"2023-03-14T13:13:44.667191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In order to incorporate weights into PCA it should be noted that the weight matrix has to be the same size as the corresponding data set. It means that the WPCA includes weights associated with each variable within each observation (in our case the variables are x, y and z coordinates of each active sensor). So that, only if each cell of our data matrix has an associated weight $w_{i}$, then we can use wpca.WPCA as intended.\n\n","metadata":{}},{"cell_type":"code","source":"#Making the weight matrix\nweights = np.stack([event_pulses.pulse_rank, event_pulses.pulse_rank, event_pulses.pulse_rank], axis=1)\n\n#Creating a WPCA object and fitting it to all pulses\nwpca = WPCA(n_components=1).fit(event_pulses[['x', 'y', 'z']], **{'weights': weights})\n\n# Calculating the first principal component vector\nvector_wpca_all = wpca.components_[0]\n\n#Calculating a distance of a point from the first principal component,that goes through the origin\nevent_pulses['distance'] = (event_pulses.x * vector_wpca_all[0] +\n                            event_pulses.y * vector_wpca_all[1] +\n                            event_pulses.z * vector_wpca_all[2]) / np.linalg.norm(vector_wpca_all)\n\n#Flip the vector direction if it points away from the neutrino origin\nif (event_pulses.loc[event_pulses.distance > 0].time.mean() >\n    event_pulses.loc[event_pulses.distance < 0].time.mean()):\n    vector_wpca_all = -1 * np.array(vector_wpca_all)\n\nvector_wpca_all = np.clip(vector_wpca_all, -1, 1)\n\n# Calculating the angles\nzen_wpca_all = np.arccos(vector_wpca_all[2])\naz_wpca_all = np.arctan2(vector_wpca_all[1], vector_wpca_all[0])\nif az_wpca_all < 0:\n    az_wpca_all = 2 * np.pi + az_wpca_all","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.677479Z","iopub.execute_input":"2023-03-14T13:13:44.678632Z","iopub.status.idle":"2023-03-14T13:13:44.770060Z","shell.execute_reply.started":"2023-03-14T13:13:44.678433Z","shell.execute_reply":"2023-03-14T13:13:44.763994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating the value of error\nwpca_all_error = angular_dist_score(azimuth_true, zenith_true, az_wpca_all, zen_wpca_all)\nprint(f'The value of error: {wpca_all_error}')","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.779809Z","iopub.execute_input":"2023-03-14T13:13:44.787244Z","iopub.status.idle":"2023-03-14T13:13:44.840629Z","shell.execute_reply.started":"2023-03-14T13:13:44.786686Z","shell.execute_reply":"2023-03-14T13:13:44.832443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating predicted vector of the event\nx_wpca_all = vector_base * vector_wpca_all[0]\ny_wpca_all = vector_base * vector_wpca_all[1]\nz_wpca_all = vector_base * vector_wpca_all[2]\nvector_wpca_all_df = pd.DataFrame({'x_pca': x_wpca_all, 'y_pca': y_wpca_all, 'z_pca': z_wpca_all})\nvector_wpca_all_df","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.851678Z","iopub.execute_input":"2023-03-14T13:13:44.855283Z","iopub.status.idle":"2023-03-14T13:13:44.911703Z","shell.execute_reply.started":"2023-03-14T13:13:44.855022Z","shell.execute_reply":"2023-03-14T13:13:44.906686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses = event_pulses.loc[~event_pulses.auxiliary]\n\n#Making the weight matrix\nweights = np.stack([event_pulses.pulse_rank, event_pulses.pulse_rank, event_pulses.pulse_rank], axis=1)\n\n#Creating a WPCA object and fitting it to non-aux pulses\nwpca = WPCA(n_components=1).fit(event_pulses[['x', 'y', 'z']], **{'weights': weights})\n\n# Calculating the first principal component vector\nvector_wpca_non = wpca.components_[0]\n\n#Calculating a distance of a point from the first principal component,that goes through the origin\nevent_pulses['distance'] = (event_pulses.x * vector_wpca_non[0] +\n                            event_pulses.y * vector_wpca_non[1] +\n                            event_pulses.z * vector_wpca_non[2]) / np.linalg.norm(vector_wpca_non)\n\n#Flip the vector direction if it points away from the neutrino origin\nif (event_pulses.loc[event_pulses.distance > 0].time.mean() >\n    event_pulses.loc[event_pulses.distance < 0].time.mean()):\n    vector_wpca_non = -1 * np.array(vector_wpca_non)\n\nvector_wpca_non = np.clip(vector_wpca_non, -1, 1)\n\n# Calculating the angles\nzen_wpca_non = np.arccos(vector_wpca_non[2])\naz_wpca_non = np.arctan2(vector_wpca_non[1], vector_wpca_non[0])\nif az_wpca_non < 0:\n    az_wpca_non = 2 * np.pi + az_wpca_non","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:44.917348Z","iopub.execute_input":"2023-03-14T13:13:44.918590Z","iopub.status.idle":"2023-03-14T13:13:45.005998Z","shell.execute_reply.started":"2023-03-14T13:13:44.918403Z","shell.execute_reply":"2023-03-14T13:13:45.000374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's calculate the error.","metadata":{}},{"cell_type":"code","source":"# Calculating the value of error\nwpca_non_error = angular_dist_score(azimuth_true, zenith_true, az_wpca_non, zen_wpca_non)\nprint(f'The value of error: {wpca_non_error}')","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:45.043884Z","iopub.execute_input":"2023-03-14T13:13:45.046205Z","iopub.status.idle":"2023-03-14T13:13:45.075132Z","shell.execute_reply.started":"2023-03-14T13:13:45.045901Z","shell.execute_reply":"2023-03-14T13:13:45.070226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check out the amount of variance explained by the weighted PC1.","metadata":{}},{"cell_type":"code","source":"for i, component in enumerate(wpca.components_):\n    print(\"{} component: {}% of initial variance\".format(i + 1, \n          round(100 * wpca.explained_variance_ratio_[i], 2)))","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:45.082382Z","iopub.execute_input":"2023-03-14T13:13:45.084549Z","iopub.status.idle":"2023-03-14T13:13:45.120501Z","shell.execute_reply.started":"2023-03-14T13:13:45.084309Z","shell.execute_reply":"2023-03-14T13:13:45.110141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating predicted vector of the event\nx_wpca_non = vector_base * vector_wpca_non[0]\ny_wpca_non = vector_base * vector_wpca_non[1]\nz_wpca_non = vector_base * vector_wpca_non[2]\nvector_wpca_non_aux_df = pd.DataFrame({'x_pca': x_wpca_non, 'y_pca': y_wpca_non, 'z_pca': z_wpca_non})\nvector_wpca_non_aux_df","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:45.129480Z","iopub.execute_input":"2023-03-14T13:13:45.131184Z","iopub.status.idle":"2023-03-14T13:13:45.191003Z","shell.execute_reply.started":"2023-03-14T13:13:45.131139Z","shell.execute_reply":"2023-03-14T13:13:45.184040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The Weighted PCA fitting visualization","metadata":{}},{"cell_type":"code","source":"fig = make_subplots(rows=1, cols=2,\n        specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],\n        subplot_titles=(f'All pulses fitting (error: {wpca_all_error:.2f})',\n                        f'Non-auxiliary pulses fitting (error: {wpca_non_error:.2f})')\n)\n\n# Plotting auxiliary pulses\nfig.add_trace(\n        go.Scatter3d(\n            x=aux_df.x, y=aux_df.y, z=aux_df.z,\n            opacity=0.75, mode='markers',\n            marker_size=aux_df.charge*10,\n            text=aux_df.charge,\n            marker=dict(color='red'),\n            name='Aux pulses'\n        ),\n        row=1, col=1\n)\n# Plotting non-auxiliary pulses\nfig.add_trace(\n        go.Scatter3d(\n            x=non_aux_df.x, y=non_aux_df.y, z=non_aux_df.z,\n            opacity=0.75, mode='markers',\n            marker_size=non_aux_df.charge*10,\n            text=non_aux_df.charge,\n            marker=dict(color='blue'),\n            name='Non-aux pulses'\n        ),\n        row='all', col='all'\n)\n\n# Plotting the true vector\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_df.x, y=vector_df.y, z=vector_df.z,\n            opacity=0.8, mode='lines', \n            line=dict(color='gray', width=4),\n            name='True vector'\n        ),  \n        row='all', col='all'\n)\n\n# Plotting the weighted PC1 vector (fitted to all pulses)\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_wpca_all_df.x_pca,\n            y=vector_wpca_all_df.y_pca,\n            z=vector_wpca_all_df.z_pca,\n            opacity=0.8, mode='lines',\n            line=dict(color='red', width=3),\n            name='Weighted PCA (all pulses)'\n        ),\n        row=1, col=1\n)\n\n# Plotting the weighted PC1 vector (fitted to non-aux pulses)\nfig.add_trace(\n        go.Scatter3d(\n            x=vector_wpca_non_aux_df.x_pca,\n            y=vector_wpca_non_aux_df.y_pca,\n            z=vector_wpca_non_aux_df.z_pca,\n            opacity=0.8, mode='lines',\n            line=dict(color='purple', width=3),\n            name='Weighted PCA (non-aux pulses)'\n        ),\n        row=1, col=2\n)\n\nfig.update_layout(\n    title=dict(text=\"Weighted PCA\", font_family='Arial',\n                pad_l=300,font_color='black', font_size=25),\n    legend=dict(font_color='black', font_size=10))\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-14T13:13:45.196599Z","iopub.execute_input":"2023-03-14T13:13:45.199363Z","iopub.status.idle":"2023-03-14T13:13:45.463872Z","shell.execute_reply.started":"2023-03-14T13:13:45.199300Z","shell.execute_reply":"2023-03-14T13:13:45.457246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, as you can see, for this event the Weighted PCA algorithm performs substantially better, than the Standard one, either we take into account all pulses in this event (aux and non-aux), or only non-auxiliary pulses. Let's take a look at the performance of these algorithms after applying them to a sample of dataset.","metadata":{}},{"cell_type":"markdown","source":"### Fitting to a sample of the training data","metadata":{}},{"cell_type":"markdown","source":"Making a function that calculates the target angles of the event.","metadata":{}},{"cell_type":"code","source":"def get_angles(event_pulses, ThisPCA=WPCA, drop_aux=True):\n    # Computing the weighted PCA\n    if ThisPCA == WPCA:\n        event_pulses = get_ranked_pulses(event_pulses)\n        \n        #if it is necessary to throw out auxiliary pulses\n        if drop_aux:\n            event_pulses = event_pulses.loc[~event_pulses.auxiliary]\n            weights = np.stack([event_pulses.pulse_rank, event_pulses.pulse_rank, \\\n                                event_pulses.pulse_rank], axis=1)\n        else:\n            weights = np.stack([event_pulses.pulse_rank, event_pulses.pulse_rank, \\\n                                event_pulses.pulse_rank], axis=1)\n        kwds = {'weights': weights}\n        \n    # Computing the stardard PCA\n    else:\n        kwds = {}\n        if drop_aux:\n            event_pulses = event_pulses.loc[~event_pulses.auxiliary]\n                \n    try:\n        # Computing the PCA vectors & variance\n        pca = ThisPCA(n_components=1).fit(event_pulses[['x', 'y', 'z']], **kwds)\n        \n        # Calculating the first principal component vector\n        vector = pca.components_[0]\n\n        # Calculating a distance of a point from a plane\n        event_pulses['distance'] = (event_pulses.x * vector[0] +\n                                    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 * np.array(vector)\n\n        vector = np.clip(vector, -1, 1)\n        \n        # Calculating the angles\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 float(zenith), float(azimuth)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:45.470360Z","iopub.execute_input":"2023-03-14T13:13:45.474362Z","iopub.status.idle":"2023-03-14T13:13:45.522267Z","shell.execute_reply.started":"2023-03-14T13:13:45.474307Z","shell.execute_reply":"2023-03-14T13:13:45.517924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_sample = train_meta_batch_125[:10000].reset_index(drop=True).copy()\ntrain_meta_sample.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:45.529908Z","iopub.execute_input":"2023-03-14T13:13:45.531051Z","iopub.status.idle":"2023-03-14T13:13:45.598821Z","shell.execute_reply.started":"2023-03-14T13:13:45.530851Z","shell.execute_reply":"2023-03-14T13:13:45.593633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sample = train_batch_125 \\\n                    .loc[:train_meta_sample.iloc[-1].last_pulse_index.astype(int)].reset_index(drop=True)\n\n# Adding coordinates of each sensor in train_sample\ntrain_sample[['x', 'y', 'z']] = sensor_geometry \\\n                    .loc[train_sample.sensor_id].reset_index()[['x', 'y', 'z']]\n\n# Zero-centering of the coordinates\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-14T13:13:45.604931Z","iopub.execute_input":"2023-03-14T13:13:45.606525Z","iopub.status.idle":"2023-03-14T13:13:46.156987Z","shell.execute_reply.started":"2023-03-14T13:13:45.606354Z","shell.execute_reply":"2023-03-14T13:13:46.155647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:46.158484Z","iopub.execute_input":"2023-03-14T13:13:46.158855Z","iopub.status.idle":"2023-03-14T13:13:46.386824Z","shell.execute_reply.started":"2023-03-14T13:13:46.158816Z","shell.execute_reply":"2023-03-14T13:13:46.385371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\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(drop=True).copy()\n    zen_wpca_non_aux, az_wpca_non_aux = get_angles(event_pulses)\n    zen_wpca_all, az_wpca_all = get_angles(event_pulses, drop_aux=False)\n    zen_pca_non_aux, az_pca_non_aux = get_angles(event_pulses, PCA)\n    zen_pca_all, az_pca_all = get_angles(event_pulses, PCA, drop_aux=False)\n    train_meta_sample.loc[i, ['zen_wpca_non_aux','az_wpca_non_aux']] = \\\n                                                            zen_wpca_non_aux, az_wpca_non_aux\n    train_meta_sample.loc[i, ['zen_wpca_all', 'az_wpca_all']] = zen_wpca_all, az_wpca_all\n    train_meta_sample.loc[i, ['zen_pca_non_aux', 'az_pca_non_aux']] = \\\n                                                            zen_pca_non_aux, az_pca_non_aux\n    train_meta_sample.loc[i, ['zen_pca_all', 'az_pca_all']] = zen_pca_all, az_pca_all","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:13:46.388404Z","iopub.execute_input":"2023-03-14T13:13:46.388832Z","iopub.status.idle":"2023-03-14T13:19:30.402033Z","shell.execute_reply.started":"2023-03-14T13:13:46.388790Z","shell.execute_reply":"2023-03-14T13:19:30.400628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_sample[['zen_wpca_non_aux', 'az_wpca_non_aux', 'zen_wpca_all', 'az_wpca_all',\n                   'zen_pca_non_aux', 'az_pca_non_aux', 'zen_pca_all', 'az_pca_all']].head(3)","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:19:30.404319Z","iopub.execute_input":"2023-03-14T13:19:30.404679Z","iopub.status.idle":"2023-03-14T13:19:30.422803Z","shell.execute_reply.started":"2023-03-14T13:19:30.404645Z","shell.execute_reply":"2023-03-14T13:19:30.421453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_group_columns = {\n    'wpca_non_aux': ['zen_wpca_non_aux', 'az_wpca_non_aux'],\n    'wpca_all': ['zen_wpca_all', 'az_wpca_all'],\n    'pca_non_aux': ['zen_pca_non_aux', 'az_pca_non_aux'],\n    'pca_all': ['zen_pca_all', 'az_pca_all']\n}\n\nerrors = {'wpca_non_aux': [], 'wpca_all': [], 'pca_non_aux': [], 'pca_all': []}\nfor key, values in result_group_columns.items():\n    error = angular_dist_score(train_meta_sample.azimuth.to_list(), \\\n                            train_meta_sample.zenith.to_list(), \\\n                            train_meta_sample[values[1]].astype(float).to_list(), \\\n                            train_meta_sample[values[0]].astype(float).to_list())\n    errors[key].append(error)\n\npd.DataFrame(errors, index=['error'])","metadata":{"execution":{"iopub.status.busy":"2023-03-14T13:19:30.425058Z","iopub.execute_input":"2023-03-14T13:19:30.425541Z","iopub.status.idle":"2023-03-14T13:19:30.502986Z","shell.execute_reply.started":"2023-03-14T13:19:30.425486Z","shell.execute_reply":"2023-03-14T13:19:30.500433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, in the above part of the notebook, we have looked at how PCA can be implemented to predict the diraction of neutrino particle in two different ways: without weightings of input data (the Standard PCA algorithm) and with it (the Weighted one). We have seen that both methods have the same goal of finding the linear sub-space which maximizes the variance of projected data and then using the PC1 vector as a target neutrino path prediction. But these methods produce different resulting reconstructions of the same original data, due to the fact that the Weighted PCA algorithm focuses on maximizing the *weighted* variance explained by each principal component, and we get non-standard principal components because their linear combination is not necessarily the best at explaining the total variance of a dataset.\n\nSo, for comparison in practice, the two algorithms were applied to the same sample of the training data and we've got the following results: in absence of noise data (without aux pulses), the Standard PCA approach performs better for our predictions than the Weighted one. But in the case of noise data (all pulses data with a high Signal/Noise ratio), on the contrary, the Weighted PCA algorithm shows better predictive power.","metadata":{}},{"cell_type":"markdown","source":"Finally, we will test the performance of the Weighted PCA algorithm on a hidden test dataset.","metadata":{}},{"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-14T13:19:30.504931Z","iopub.execute_input":"2023-03-14T13:19:30.505589Z","iopub.status.idle":"2023-03-14T13:19:30.537516Z","shell.execute_reply.started":"2023-03-14T13:19:30.505530Z","shell.execute_reply":"2023-03-14T13:19:30.536085Z"},"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-14T13:19:30.539816Z","iopub.execute_input":"2023-03-14T13:19:30.540306Z","iopub.status.idle":"2023-03-14T13:19:30.548790Z","shell.execute_reply.started":"2023-03-14T13:19:30.540255Z","shell.execute_reply":"2023-03-14T13:19:30.547401Z"},"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-14T13:19:30.550423Z","iopub.execute_input":"2023-03-14T13:19:30.551272Z","iopub.status.idle":"2023-03-14T13:19:30.633002Z","shell.execute_reply.started":"2023-03-14T13:19:30.551227Z","shell.execute_reply":"2023-03-14T13:19:30.631238Z"},"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-14T13:19:30.637593Z","iopub.execute_input":"2023-03-14T13:19:30.638268Z","iopub.status.idle":"2023-03-14T13:19:30.660120Z","shell.execute_reply.started":"2023-03-14T13:19:30.638222Z","shell.execute_reply":"2023-03-14T13:19:30.658712Z"},"trusted":true},"execution_count":null,"outputs":[]}]}