{"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":{"papermill":{"duration":0.01442,"end_time":"2023-03-10T09:03:51.425482","exception":false,"start_time":"2023-03-10T09:03:51.411062","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"This notebook is based on existing interesting research (<a href=\"https://www.kaggle.com/code/averkovanika/standard-vs-weighted-pca-who-s-won\">Notebook link</a>)\n\nAccording existing research We assign ranks to pulses using the following ranking system:\n\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\n\nIn this notebook I've used different weights in order to decrease the error. This aprouch ","metadata":{}},{"cell_type":"code","source":"!pip install wpca","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:43:41.804123Z","iopub.execute_input":"2023-04-03T18:43:41.804702Z","iopub.status.idle":"2023-04-03T18:43:57.023616Z","shell.execute_reply.started":"2023-04-03T18:43:41.804656Z","shell.execute_reply":"2023-04-03T18:43:57.022582Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"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":{"papermill":{"duration":10.11877,"end_time":"2023-03-10T09:04:01.588006","exception":false,"start_time":"2023-03-10T09:03:51.469236","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:43:57.025397Z","iopub.execute_input":"2023-04-03T18:43:57.025762Z","iopub.status.idle":"2023-04-03T18:43:57.031661Z","shell.execute_reply.started":"2023-04-03T18:43:57.025722Z","shell.execute_reply":"2023-04-03T18:43:57.030395Z"},"_kg_hide-input":true,"_kg_hide-output":true,"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":{"papermill":{"duration":2.986815,"end_time":"2023-03-10T09:04:04.588534","exception":false,"start_time":"2023-03-10T09:04:01.601719","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:43:57.032987Z","iopub.execute_input":"2023-04-03T18:43:57.033277Z","iopub.status.idle":"2023-04-03T18:43:59.441270Z","shell.execute_reply.started":"2023-04-03T18:43:57.033248Z","shell.execute_reply":"2023-04-03T18:43:59.440214Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH_DATASETS = '/kaggle/input/icecube-neutrinos-in-deep-ice/'","metadata":{"papermill":{"duration":0.023033,"end_time":"2023-03-10T09:04:04.624318","exception":false,"start_time":"2023-03-10T09:04:04.601285","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:43:59.443536Z","iopub.execute_input":"2023-04-03T18:43:59.443853Z","iopub.status.idle":"2023-04-03T18:43:59.451711Z","shell.execute_reply.started":"2023-04-03T18:43:59.443820Z","shell.execute_reply":"2023-04-03T18:43:59.450629Z"},"_kg_hide-input":true,"_kg_hide-output":true,"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":false,"papermill":{"duration":0.02887,"end_time":"2023-03-10T09:04:04.665810","exception":false,"start_time":"2023-03-10T09:04:04.636940","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:43:59.453911Z","iopub.execute_input":"2023-04-03T18:43:59.454332Z","iopub.status.idle":"2023-04-03T18:43:59.468672Z","shell.execute_reply.started":"2023-04-03T18:43:59.454288Z","shell.execute_reply":"2023-04-03T18:43:59.467520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data loading","metadata":{"papermill":{"duration":0.012555,"end_time":"2023-03-10T09:04:04.692072","exception":false,"start_time":"2023-03-10T09:04:04.679517","status":"completed"},"tags":[]}},{"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-04-03T18:43:59.470086Z","iopub.execute_input":"2023-04-03T18:43:59.470402Z","iopub.status.idle":"2023-04-03T18:44:40.699390Z","shell.execute_reply.started":"2023-04-03T18:43:59.470371Z","shell.execute_reply":"2023-04-03T18:44:40.698199Z"},"_kg_hide-output":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train_meta\n","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:40.701035Z","iopub.execute_input":"2023-04-03T18:44:40.701558Z","iopub.status.idle":"2023-04-03T18:44:40.711508Z","shell.execute_reply.started":"2023-04-03T18:44:40.701520Z","shell.execute_reply":"2023-04-03T18:44:40.710270Z"},"_kg_hide-input":true,"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-04-03T18:44:40.712914Z","iopub.execute_input":"2023-04-03T18:44:40.713272Z","iopub.status.idle":"2023-04-03T18:44:40.744880Z","shell.execute_reply.started":"2023-04-03T18:44:40.713221Z","shell.execute_reply":"2023-04-03T18:44:40.743802Z"},"_kg_hide-input":true,"_kg_hide-output":true,"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-04-03T18:44:40.746053Z","iopub.execute_input":"2023-04-03T18:44:40.746353Z","iopub.status.idle":"2023-04-03T18:44:44.925437Z","shell.execute_reply.started":"2023-04-03T18:44:40.746324Z","shell.execute_reply":"2023-04-03T18:44:44.924014Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's begin with:","metadata":{}},{"cell_type":"code","source":"event_id_1 = 403533054","metadata":{"papermill":{"duration":0.020328,"end_time":"2023-03-10T09:04:46.056315","exception":false,"start_time":"2023-03-10T09:04:46.035987","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:44.929218Z","iopub.execute_input":"2023-04-03T18:44:44.929657Z","iopub.status.idle":"2023-04-03T18:44:44.934409Z","shell.execute_reply.started":"2023-04-03T18:44:44.929622Z","shell.execute_reply":"2023-04-03T18:44:44.933120Z"},"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()","metadata":{"papermill":{"duration":0.283649,"end_time":"2023-03-10T09:04:46.353910","exception":false,"start_time":"2023-03-10T09:04:46.070261","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:44.935792Z","iopub.execute_input":"2023-04-03T18:44:44.936592Z","iopub.status.idle":"2023-04-03T18:44:45.494193Z","shell.execute_reply.started":"2023-04-03T18:44:44.936545Z","shell.execute_reply":"2023-04-03T18:44:45.492881Z"},"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Adding zero-centered coordinates of each sensor of the event","metadata":{"papermill":{"duration":0.013455,"end_time":"2023-03-10T09:04:46.380775","exception":false,"start_time":"2023-03-10T09:04:46.367320","status":"completed"},"tags":[]}},{"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.auxiliary.value_counts()","metadata":{"papermill":{"duration":0.03733,"end_time":"2023-03-10T09:04:46.430792","exception":false,"start_time":"2023-03-10T09:04:46.393462","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.495983Z","iopub.execute_input":"2023-04-03T18:44:45.497092Z","iopub.status.idle":"2023-04-03T18:44:45.518728Z","shell.execute_reply.started":"2023-04-03T18:44:45.497051Z","shell.execute_reply":"2023-04-03T18:44:45.517369Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Calculating the true vector of the event","metadata":{"papermill":{"duration":0.012916,"end_time":"2023-03-10T09:04:46.458018","exception":false,"start_time":"2023-03-10T09:04:46.445102","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":0.035608,"end_time":"2023-03-10T09:04:46.507456","exception":false,"start_time":"2023-03-10T09:04:46.471848","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.520089Z","iopub.execute_input":"2023-04-03T18:44:45.521248Z","iopub.status.idle":"2023-04-03T18:44:45.541268Z","shell.execute_reply.started":"2023-04-03T18:44:45.521196Z","shell.execute_reply":"2023-04-03T18:44:45.540019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The Weighted PCA","metadata":{"papermill":{"duration":0.014027,"end_time":"2023-03-10T09:04:47.442983","exception":false,"start_time":"2023-03-10T09:04:47.428956","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"According this research (<a href=\"https://www.kaggle.com/code/averkovanika/standard-vs-weighted-pca-who-s-won\">Notebook link</a>) We assign ranks to pulses using the following ranking system:\n\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.","metadata":{"papermill":{"duration":0.014641,"end_time":"2023-03-10T09:04:47.472041","exception":false,"start_time":"2023-03-10T09:04:47.457400","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"And one more important thing that we need to consider: \n- The 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\".\n ","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.01435,"end_time":"2023-03-10T09:04:47.501029","exception":false,"start_time":"2023-03-10T09:04:47.486679","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"### So let's calculate the error using these conditions","metadata":{"execution":{"iopub.status.busy":"2023-03-19T10:48:56.918729Z","iopub.execute_input":"2023-03-19T10:48:56.919214Z","iopub.status.idle":"2023-03-19T10:48:56.924785Z","shell.execute_reply.started":"2023-03-19T10:48:56.919152Z","shell.execute_reply":"2023-03-19T10:48:56.923624Z"}}},{"cell_type":"markdown","source":"* Calculating the value of this threshold transit time for neutrinos - 6200 ns.","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.014162,"end_time":"2023-03-10T09:04:47.529893","exception":false,"start_time":"2023-03-10T09:04:47.515731","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":0.028546,"end_time":"2023-03-10T09:04:47.572762","exception":false,"start_time":"2023-03-10T09:04:47.544216","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.542666Z","iopub.execute_input":"2023-04-03T18:44:45.543005Z","iopub.status.idle":"2023-04-03T18:44:45.555007Z","shell.execute_reply.started":"2023-04-03T18:44:45.542972Z","shell.execute_reply":"2023-04-03T18:44:45.553931Z"},"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":{"papermill":{"duration":0.014333,"end_time":"2023-03-10T09:04:47.601886","exception":false,"start_time":"2023-03-10T09:04:47.587553","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":0.025927,"end_time":"2023-03-10T09:04:47.643638","exception":false,"start_time":"2023-03-10T09:04:47.617711","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.556295Z","iopub.execute_input":"2023-04-03T18:44:45.556633Z","iopub.status.idle":"2023-04-03T18:44:45.563876Z","shell.execute_reply.started":"2023-04-03T18:44:45.556602Z","shell.execute_reply":"2023-04-03T18:44:45.562700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Applying the function to our event.","metadata":{"papermill":{"duration":0.015564,"end_time":"2023-03-10T09:04:47.674520","exception":false,"start_time":"2023-03-10T09:04:47.658956","status":"completed"},"tags":[]}},{"cell_type":"code","source":"event_pulses = get_ranked_pulses(event_pulses)\nevent_pulses.pulse_rank.value_counts()","metadata":{"papermill":{"duration":0.037136,"end_time":"2023-03-10T09:04:47.728284","exception":false,"start_time":"2023-03-10T09:04:47.691148","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.565909Z","iopub.execute_input":"2023-04-03T18:44:45.566320Z","iopub.status.idle":"2023-04-03T18:44:45.585042Z","shell.execute_reply.started":"2023-04-03T18:44:45.566288Z","shell.execute_reply":"2023-04-03T18:44:45.583731Z"},"_kg_hide-output":true,"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":{"papermill":{"duration":0.01494,"end_time":"2023-03-10T09:04:47.758180","exception":false,"start_time":"2023-03-10T09:04:47.743240","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":0.043892,"end_time":"2023-03-10T09:04:47.817167","exception":false,"start_time":"2023-03-10T09:04:47.773275","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.586556Z","iopub.execute_input":"2023-04-03T18:44:45.587033Z","iopub.status.idle":"2023-04-03T18:44:45.620462Z","shell.execute_reply.started":"2023-04-03T18:44:45.586997Z","shell.execute_reply":"2023-04-03T18:44:45.618723Z"},"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":{"papermill":{"duration":0.029703,"end_time":"2023-03-10T09:04:47.866175","exception":false,"start_time":"2023-03-10T09:04:47.836472","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.622986Z","iopub.execute_input":"2023-04-03T18:44:45.623991Z","iopub.status.idle":"2023-04-03T18:44:45.634140Z","shell.execute_reply.started":"2023-04-03T18:44:45.623919Z","shell.execute_reply":"2023-04-03T18:44:45.631870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Okey. So the value of error is 1.5519125875455784.","metadata":{}},{"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":{"papermill":{"duration":0.036302,"end_time":"2023-03-10T09:04:47.921417","exception":false,"start_time":"2023-03-10T09:04:47.885115","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.636876Z","iopub.execute_input":"2023-04-03T18:44:45.638070Z","iopub.status.idle":"2023-04-03T18:44:45.674746Z","shell.execute_reply.started":"2023-04-03T18:44:45.637876Z","shell.execute_reply":"2023-04-03T18:44:45.672929Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### But why did we decide that the weights of different types of pulses should differ by one?\n\nFor example, let's make the weight for pulses:\n\n- non-auxiliary pulses within the valid time window -- 5 insted 3\n- auxiliary pulses within the valid time window -- 1 insted 2\n- other pulses -- 0\n\n#### So in other words we will slightly change the ratio of the weights of these pulses\n\nMaking new function that assigns new weights to all pulses in the event.","metadata":{}},{"cell_type":"code","source":"def get_weight_all(group):\n\n    if group == 2:\n        return 1\n    if group == 3:\n        return 5\n    else:\n        return 0","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.682464Z","iopub.execute_input":"2023-04-03T18:44:45.687588Z","iopub.status.idle":"2023-04-03T18:44:45.699660Z","shell.execute_reply.started":"2023-04-03T18:44:45.687491Z","shell.execute_reply":"2023-04-03T18:44:45.697955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Applying the function to our event.","metadata":{}},{"cell_type":"code","source":"event_pulses['weight'] = event_pulses['pulse_rank'].apply(get_weight_all)\nevent_pulses['weight'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.707903Z","iopub.execute_input":"2023-04-03T18:44:45.712698Z","iopub.status.idle":"2023-04-03T18:44:45.731176Z","shell.execute_reply.started":"2023-04-03T18:44:45.712620Z","shell.execute_reply":"2023-04-03T18:44:45.729856Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Let's calculate the error again.","metadata":{}},{"cell_type":"code","source":"#Making the weight matrix\nweights_new = np.stack([event_pulses.weight, event_pulses.weight, event_pulses.weight], 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_new})\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-04-03T18:44:45.733066Z","iopub.execute_input":"2023-04-03T18:44:45.733838Z","iopub.status.idle":"2023-04-03T18:44:45.748924Z","shell.execute_reply.started":"2023-04-03T18:44:45.733800Z","shell.execute_reply":"2023-04-03T18:44:45.747887Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating the value of error\nwpca_all_error_new = angular_dist_score(azimuth_true, zenith_true, az_wpca_all, zen_wpca_all)\nprint(f'The value of error: {wpca_all_error_new}')","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.750194Z","iopub.execute_input":"2023-04-03T18:44:45.751187Z","iopub.status.idle":"2023-04-03T18:44:45.761553Z","shell.execute_reply.started":"2023-04-03T18:44:45.751147Z","shell.execute_reply":"2023-04-03T18:44:45.760405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### There's a difference, isn't there?","metadata":{}},{"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_new = pd.DataFrame({'x_pca': x_wpca_all, 'y_pca': y_wpca_all, 'z_pca': z_wpca_all})\nvector_wpca_all_df_new","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.763061Z","iopub.execute_input":"2023-04-03T18:44:45.765149Z","iopub.status.idle":"2023-04-03T18:44:45.783477Z","shell.execute_reply.started":"2023-04-03T18:44:45.765100Z","shell.execute_reply":"2023-04-03T18:44:45.782589Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Let's check for non-aux pulses.","metadata":{}},{"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":{"papermill":{"duration":0.034025,"end_time":"2023-03-10T09:04:47.971540","exception":false,"start_time":"2023-03-10T09:04:47.937515","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.786525Z","iopub.execute_input":"2023-04-03T18:44:45.787565Z","iopub.status.idle":"2023-04-03T18:44:45.804367Z","shell.execute_reply.started":"2023-04-03T18:44:45.787486Z","shell.execute_reply":"2023-04-03T18:44:45.803005Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"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":{"papermill":{"duration":0.024795,"end_time":"2023-03-10T09:04:48.042233","exception":false,"start_time":"2023-03-10T09:04:48.017438","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.806158Z","iopub.execute_input":"2023-04-03T18:44:45.806945Z","iopub.status.idle":"2023-04-03T18:44:45.815827Z","shell.execute_reply.started":"2023-04-03T18:44:45.806891Z","shell.execute_reply":"2023-04-03T18:44:45.814640Z"},"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":{"papermill":{"duration":0.039896,"end_time":"2023-03-10T09:04:48.171155","exception":false,"start_time":"2023-03-10T09:04:48.131259","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.816971Z","iopub.execute_input":"2023-04-03T18:44:45.817288Z","iopub.status.idle":"2023-04-03T18:44:45.834188Z","shell.execute_reply.started":"2023-04-03T18:44:45.817257Z","shell.execute_reply":"2023-04-03T18:44:45.833101Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses.pulse_rank.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.839609Z","iopub.execute_input":"2023-04-03T18:44:45.840575Z","iopub.status.idle":"2023-04-03T18:44:45.852522Z","shell.execute_reply.started":"2023-04-03T18:44:45.840531Z","shell.execute_reply":"2023-04-03T18:44:45.851207Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### There are less pulses. It was to be expected. \n\nLets make another one function with weights for test.","metadata":{}},{"cell_type":"code","source":"def get_weight_non(group):\n    if group == event_pulses.pulse_rank.value_counts().index[0]:\n        return 1\n    if group == event_pulses.pulse_rank.value_counts().index[1]:       \n        return 0.1\n    else:\n        return 0.001","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.853837Z","iopub.execute_input":"2023-04-03T18:44:45.854178Z","iopub.status.idle":"2023-04-03T18:44:45.860312Z","shell.execute_reply.started":"2023-04-03T18:44:45.854145Z","shell.execute_reply":"2023-04-03T18:44:45.859185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Applying the function to our event.","metadata":{}},{"cell_type":"code","source":"event_pulses['weight'] = event_pulses['pulse_rank'].apply(get_weight_non)\nevent_pulses.head()","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.861845Z","iopub.execute_input":"2023-04-03T18:44:45.862165Z","iopub.status.idle":"2023-04-03T18:44:45.890835Z","shell.execute_reply.started":"2023-04-03T18:44:45.862133Z","shell.execute_reply":"2023-04-03T18:44:45.889665Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Let's calculate the error with new weights.","metadata":{}},{"cell_type":"code","source":"#Making the weight matrix\nweights_new = np.stack([event_pulses.weight, event_pulses.weight, event_pulses.weight], 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_new})\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-04-03T18:44:45.892194Z","iopub.execute_input":"2023-04-03T18:44:45.893067Z","iopub.status.idle":"2023-04-03T18:44:45.910035Z","shell.execute_reply.started":"2023-04-03T18:44:45.893031Z","shell.execute_reply":"2023-04-03T18:44:45.908857Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculating the value of error\nwpca_non_error_new = angular_dist_score(azimuth_true, zenith_true, az_wpca_non, zen_wpca_non)\nprint(f'The value of error: {wpca_non_error_new}')","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.911362Z","iopub.execute_input":"2023-04-03T18:44:45.914706Z","iopub.status.idle":"2023-04-03T18:44:45.919670Z","shell.execute_reply.started":"2023-04-03T18:44:45.914663Z","shell.execute_reply":"2023-04-03T18:44:45.918786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The same situation! The error has become less.","metadata":{}},{"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_new = pd.DataFrame({'x_pca': x_wpca_non, 'y_pca': y_wpca_non, 'z_pca': z_wpca_non})\nvector_wpca_non_aux_df_new","metadata":{"execution":{"iopub.status.busy":"2023-04-03T18:44:45.920981Z","iopub.execute_input":"2023-04-03T18:44:45.921577Z","iopub.status.idle":"2023-04-03T18:44:45.943792Z","shell.execute_reply.started":"2023-04-03T18:44:45.921540Z","shell.execute_reply":"2023-04-03T18:44:45.942662Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The Weighted PCA visualization with old weights","metadata":{"papermill":{"duration":0.014957,"end_time":"2023-03-10T09:04:48.205134","exception":false,"start_time":"2023-03-10T09:04:48.190177","status":"completed"},"tags":[]}},{"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: {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 with pulse_rank as weights\", 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,"papermill":{"duration":0.066021,"end_time":"2023-03-10T09:04:48.285951","exception":false,"start_time":"2023-03-10T09:04:48.219930","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:45.945409Z","iopub.execute_input":"2023-04-03T18:44:45.945804Z","iopub.status.idle":"2023-04-03T18:44:46.419804Z","shell.execute_reply.started":"2023-04-03T18:44:45.945768Z","shell.execute_reply":"2023-04-03T18:44:46.418719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The Weighted PCA visualization new weights","metadata":{"execution":{"iopub.status.busy":"2023-03-19T12:19:29.463912Z","iopub.execute_input":"2023-03-19T12:19:29.464593Z","iopub.status.idle":"2023-03-19T12:19:29.471075Z","shell.execute_reply.started":"2023-03-19T12:19:29.464529Z","shell.execute_reply":"2023-03-19T12:19:29.469905Z"}}},{"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_new:.4f})',\n                        f'Non-auxiliary pulses fitting (error: {wpca_non_error_new:.4f})')\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_new.x_pca,\n            y=vector_wpca_all_df_new.y_pca,\n            z=vector_wpca_all_df_new.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_new.x_pca,\n            y=vector_wpca_non_aux_df_new.y_pca,\n            z=vector_wpca_non_aux_df_new.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 with new weights\", 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":{"execution":{"iopub.status.busy":"2023-04-03T18:44:46.421712Z","iopub.execute_input":"2023-04-03T18:44:46.422174Z","iopub.status.idle":"2023-04-03T18:44:46.480716Z","shell.execute_reply.started":"2023-04-03T18:44:46.422128Z","shell.execute_reply":"2023-04-03T18:44:46.479511Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fitting to a sample of the training data","metadata":{"papermill":{"duration":0.014928,"end_time":"2023-03-10T09:04:48.316463","exception":false,"start_time":"2023-03-10T09:04:48.301535","status":"completed"},"tags":[],"_kg_hide-input":true,"_kg_hide-output":true}},{"cell_type":"markdown","source":"Making a function that calculates the target angles of the event.","metadata":{"papermill":{"duration":0.015197,"end_time":"2023-03-10T09:04:48.346746","exception":false,"start_time":"2023-03-10T09:04:48.331549","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def get_angles(event_pulses, wpca_type, drop_aux=True):\n\n    event_pulses = get_ranked_pulses(event_pulses)\n\n    \n    # Computing the weighted PCA\n    if wpca_type == 1:\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    elif wpca_type == 2:\n        #if it is necessary to throw out auxiliary pulses\n        if drop_aux:\n            event_pulses = event_pulses.loc[~event_pulses.auxiliary]\n            event_pulses['weight'] = event_pulses['pulse_rank'].apply(get_weight_non) # было get_weight_non\n            weights = np.stack([event_pulses.weight, event_pulses.weight, \\\n                                event_pulses.weight], axis=1)\n        else:\n            event_pulses['weight'] = event_pulses['pulse_rank'].apply(get_weight_all)\n            weights = np.stack([event_pulses.weight, event_pulses.weight, \\\n                                event_pulses.weight], axis=1)\n        kwds = {'weights': weights}\n                 \n    try:\n        # Computing the weighted PCA vectors & variance\n        wpca = WPCA(n_components=1).fit(event_pulses[['x', 'y', 'z']], **kwds)\n        \n        # Calculating the first principal component vector\n        vector = wpca.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":{"papermill":{"duration":0.030524,"end_time":"2023-03-10T09:04:48.392671","exception":false,"start_time":"2023-03-10T09:04:48.362147","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:46.482226Z","iopub.execute_input":"2023-04-03T18:44:46.482550Z","iopub.status.idle":"2023-04-03T18:44:46.496757Z","shell.execute_reply.started":"2023-04-03T18:44:46.482520Z","shell.execute_reply":"2023-04-03T18:44:46.495930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_sample = train_meta_batch_125[:10000].reset_index(drop=True).copy()  # 10000\ntrain_meta_sample.head(3)","metadata":{"papermill":{"duration":0.032999,"end_time":"2023-03-10T09:04:48.441691","exception":false,"start_time":"2023-03-10T09:04:48.408692","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:46.498289Z","iopub.execute_input":"2023-04-03T18:44:46.498760Z","iopub.status.idle":"2023-04-03T18:44:46.523159Z","shell.execute_reply.started":"2023-04-03T18:44:46.498723Z","shell.execute_reply":"2023-04-03T18:44:46.521983Z"},"_kg_hide-output":true,"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":{"papermill":{"duration":0.273943,"end_time":"2023-03-10T09:04:48.731203","exception":false,"start_time":"2023-03-10T09:04:48.457260","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:46.524538Z","iopub.execute_input":"2023-04-03T18:44:46.524853Z","iopub.status.idle":"2023-04-03T18:44:46.979363Z","shell.execute_reply.started":"2023-04-03T18:44:46.524822Z","shell.execute_reply":"2023-04-03T18:44:46.978132Z"},"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()","metadata":{"papermill":{"duration":0.154408,"end_time":"2023-03-10T09:04:48.901557","exception":false,"start_time":"2023-03-10T09:04:48.747149","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:46.981214Z","iopub.execute_input":"2023-04-03T18:44:46.981825Z","iopub.status.idle":"2023-04-03T18:44:47.122659Z","shell.execute_reply.started":"2023-04-03T18:44:46.981787Z","shell.execute_reply":"2023-04-03T18:44:47.121277Z"},"_kg_hide-input":true,"_kg_hide-output":true,"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, 1)\n    zen_wpca_all, az_wpca_all = get_angles(event_pulses, 1, drop_aux=False)\n    zen_wpca_non_aux_new_weight, az_wpca_non_aux_new_weight = get_angles(event_pulses, 2)\n    zen_wpca_all_new_weight, az_wpca_all_new_weight = get_angles(event_pulses, 2, 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_wpca_non_aux_new_weight', 'az_wpca_non_aux_new_weight']] = \\\n                                                            zen_wpca_non_aux_new_weight, az_wpca_non_aux_new_weight\n    train_meta_sample.loc[i, ['zen_wpca_all_new_weight', 'az_wpca_all_new_weight']] = \\\n                                                            zen_wpca_all_new_weight, az_wpca_all_new_weight","metadata":{"papermill":{"duration":265.564626,"end_time":"2023-03-10T09:09:14.483180","exception":false,"start_time":"2023-03-10T09:04:48.918554","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T18:44:47.124093Z","iopub.execute_input":"2023-04-03T18:44:47.125162Z","iopub.status.idle":"2023-04-03T19:02:15.725256Z","shell.execute_reply.started":"2023-04-03T18:44:47.125125Z","shell.execute_reply":"2023-04-03T19:02:15.723563Z"},"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_wpca_non_aux_new_weight', 'az_wpca_non_aux_new_weight', \n                   'zen_wpca_all_new_weight', 'az_wpca_all_new_weight']].head(3)","metadata":{"papermill":{"duration":0.036137,"end_time":"2023-03-10T09:09:14.537636","exception":false,"start_time":"2023-03-10T09:09:14.501499","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T19:02:15.726999Z","iopub.execute_input":"2023-04-03T19:02:15.727405Z","iopub.status.idle":"2023-04-03T19:02:15.747289Z","shell.execute_reply.started":"2023-04-03T19:02:15.727368Z","shell.execute_reply":"2023-04-03T19:02:15.745668Z"},"_kg_hide-output":true,"_kg_hide-input":true,"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    'wpca_non_aux_new_weight': ['zen_wpca_non_aux_new_weight', 'az_wpca_non_aux_new_weight'],\n    'wpca_all_new_weight': ['zen_wpca_all_new_weight', 'az_wpca_all_new_weight']\n}\n\nerrors = {'wpca_non_aux': [], 'wpca_non_aux_new_weight': [], 'wpca_all': [], 'wpca_all_new_weight': []}\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)","metadata":{"papermill":{"duration":0.078664,"end_time":"2023-03-10T09:09:14.632932","exception":false,"start_time":"2023-03-10T09:09:14.554268","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-04-03T19:02:15.749025Z","iopub.execute_input":"2023-04-03T19:02:15.749524Z","iopub.status.idle":"2023-04-03T19:02:15.821663Z","shell.execute_reply.started":"2023-04-03T19:02:15.749445Z","shell.execute_reply":"2023-04-03T19:02:15.820682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### As a result of our research the errors have been less especially for noisy data. If we have chosen the right weights for the pulses, ranking  by category is reasonable.","metadata":{}}]}