{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# IceCube🧊: Neutrino\n\nThe goal of this competition is to predict a neutrino particle’s direction. You will develop a model based on data from the \"IceCube\" detector, which observes the cosmos from deep within the South Pole ice.\n\nThe IceCube Neutrino Observatory is the first detector of its kind, encompassing a cubic kilometer of ice and designed to search for the nearly massless neutrinos. An international group of scientists is responsible for the scientific research that makes up the IceCube Collaboration.\n\nBy making the process faster and more precise, you'll help improve the reconstruction of neutrinos. As a result, we could gain a clearer image of our universe.","metadata":{}},{"cell_type":"markdown","source":"# EDA Phase","metadata":{}},{"cell_type":"code","source":"# %matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:52:11.802042Z","iopub.execute_input":"2023-03-24T17:52:11.802506Z","iopub.status.idle":"2023-03-24T17:52:11.837183Z","shell.execute_reply.started":"2023-03-24T17:52:11.802461Z","shell.execute_reply":"2023-03-24T17:52:11.835913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import numpy as np\n# import pandas as pd\n# import plotly.express as px\n# from sklearn.decomposition import PCA\n# import plotly.graph_objects as go\n# import matplotlib.pyplot as plt\n# import multiprocessing","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:25:17.337692Z","iopub.execute_input":"2023-03-26T19:25:17.338145Z","iopub.status.idle":"2023-03-26T19:25:21.784930Z","shell.execute_reply.started":"2023-03-26T19:25:17.338104Z","shell.execute_reply":"2023-03-26T19:25:21.783164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Install nb_black for autoformating\n# !pip install nb_black --quiet\n# %load_ext lab_black","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:52:16.355968Z","iopub.execute_input":"2023-03-24T17:52:16.356592Z","iopub.status.idle":"2023-03-24T17:52:32.953300Z","shell.execute_reply.started":"2023-03-24T17:52:16.356531Z","shell.execute_reply":"2023-03-24T17:52:32.951662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Browse meta data:\n\n- *[train/test]_meta.parquet*\n- *batch_id (int)*: the ID of the batch the event was placed into.\n- *event_id (int)*: the event ID.\n- *[first/last]_pulse_index (int)*: index of the first/last row in the features dataframe belonging to this event.\n- *[azimuth/zenith] (float32)*: the *[azimuth/zenith]* angle in radians of the neutrino. A value between 0 and 2*pi for the azimuth and 0 and pi for zenith. The target columns. Not provided for the test set. The direction vector represented by zenith and azimuth points to where the neutrino came from.","metadata":{}},{"cell_type":"markdown","source":"# TL;DR\nI perform a PCA on the pulses coordinates and take the first principal component as the neutrino's direction. The PCA is performed only on the non-auxiliary ('true') pulses if there are enough (>15); otherwise, I perform it on all the data. Then, I calculate the average time of the 'true' pulses on the negative and positive sides of the vector and flip its direction if it points to the wrong one. (lower times should be on the positive side since it is supposed to point to the source of the neutrino). This analysis is done only for the 'true' pulses since the auxiliary pulses times are too noisy for this analysis without further refinement. Validation score 1.282, LB score 1.274","metadata":{}},{"cell_type":"markdown","source":"# Train Meta ","metadata":{}},{"cell_type":"code","source":"# BASE_DIR = \"../input/icecube-neutrinos-in-deep-ice/\"\n# train_meta = pd.read_parquet(f\"{BASE_DIR}/train_meta.parquet\")\n# # train_meta.head()\n# print(\"Total number of events: \" + str(len(train_meta)))","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:25:26.683884Z","iopub.execute_input":"2023-03-26T19:25:26.684329Z","iopub.status.idle":"2023-03-26T19:26:07.845077Z","shell.execute_reply.started":"2023-03-26T19:25:26.684289Z","shell.execute_reply":"2023-03-26T19:26:07.843693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_meta.shape  # (131953924, 6)\n# train_meta[\"event_id\"].nunique()  # 131953924","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:53:24.050117Z","iopub.execute_input":"2023-03-24T17:53:24.050509Z","iopub.status.idle":"2023-03-24T17:53:24.056885Z","shell.execute_reply.started":"2023-03-24T17:53:24.050471Z","shell.execute_reply":"2023-03-24T17:53:24.055423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train batch 1","metadata":{}},{"cell_type":"code","source":"# train_batch_1 = pd.read_parquet(f\"{BASE_DIR}/train/batch_1.parquet\")\n# train_batch_1.head()\n# print(\"Number of pulses: \" + str(len(train_batch_1)))\n# print(\"Number of events: \" + str(len(train_batch_1.index.unique())))","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:12.031361Z","iopub.execute_input":"2023-03-26T19:26:12.031892Z","iopub.status.idle":"2023-03-26T19:26:15.939875Z","shell.execute_reply.started":"2023-03-26T19:26:12.031848Z","shell.execute_reply":"2023-03-26T19:26:15.938264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Sensor geometry","metadata":{}},{"cell_type":"code","source":"# sensor_geometry = pd.read_csv(f\"{BASE_DIR}/sensor_geometry.csv\")\n# sensor_geometry.head()\n# plt.scatter(sensor_geometry.x, sensor_geometry.y)\n# # sensor_geometry.shape  # (5160, 4)","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:19.229492Z","iopub.execute_input":"2023-03-26T19:26:19.231169Z","iopub.status.idle":"2023-03-26T19:26:19.610525Z","shell.execute_reply.started":"2023-03-26T19:26:19.231099Z","shell.execute_reply":"2023-03-26T19:26:19.608831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fig = px.scatter_3d(sensor_geometry, x=\"x\", y=\"y\", z=\"z\", opacity=0.5)\n# fig.update_traces(marker_size=2)\n# fig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:22.080754Z","iopub.execute_input":"2023-03-26T19:26:22.081263Z","iopub.status.idle":"2023-03-26T19:26:24.203000Z","shell.execute_reply.started":"2023-03-26T19:26:22.081217Z","shell.execute_reply":"2023-03-26T19:26:24.201857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test Meta","metadata":{}},{"cell_type":"code","source":"# test_meta = pd.read_parquet(f\"{BASE_DIR}/test_meta.parquet\")\n# test_meta.head()\n# print(\"Total number of events: \" + str(len(test_meta)))","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:33.119259Z","iopub.execute_input":"2023-03-26T19:26:33.120093Z","iopub.status.idle":"2023-03-26T19:26:33.139159Z","shell.execute_reply.started":"2023-03-26T19:26:33.120040Z","shell.execute_reply":"2023-03-26T19:26:33.137676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Understanding the data\n- I will explore event 67 (number three in train_meta) as a case event","metadata":{}},{"cell_type":"code","source":"# train_meta","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:38.906025Z","iopub.execute_input":"2023-03-26T19:26:38.906507Z","iopub.status.idle":"2023-03-26T19:26:38.945403Z","shell.execute_reply.started":"2023-03-26T19:26:38.906469Z","shell.execute_reply":"2023-03-26T19:26:38.943546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_idx = 3\n# case_event_pulse = train_batch_1.query(\"event_id == 67\").copy()\n# print(len(case_event_pulse))\n# case_event_pulse.tail()","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:41.316561Z","iopub.execute_input":"2023-03-26T19:26:41.317031Z","iopub.status.idle":"2023-03-26T19:26:41.882628Z","shell.execute_reply.started":"2023-03-26T19:26:41.316988Z","shell.execute_reply":"2023-03-26T19:26:41.881259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse.shape, train_batch_1.query(\"event_id == 67\").shape  # ((142, 4), (142, 4))\n\n# train_meta.iloc[case_event_idx]","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:53:32.133497Z","iopub.execute_input":"2023-03-24T17:53:32.133891Z","iopub.status.idle":"2023-03-24T17:53:32.140439Z","shell.execute_reply.started":"2023-03-24T17:53:32.133852Z","shell.execute_reply":"2023-03-24T17:53:32.138924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Finding the coordinates of each pulse from the sensor_geometry data and mean normalizing","metadata":{}},{"cell_type":"code","source":"# case_event_pulse","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:53:32.142302Z","iopub.execute_input":"2023-03-24T17:53:32.143347Z","iopub.status.idle":"2023-03-24T17:53:32.167898Z","shell.execute_reply.started":"2023-03-24T17:53:32.143300Z","shell.execute_reply":"2023-03-24T17:53:32.166022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse.set_index(\"time\").groupby(\"auxiliary\")[\"charge\"].plot(\n#     style=\".\", figsize=(15, 5), title=\"Event 67 by time and auxiliary\"\n# )\n# plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:45.244463Z","iopub.execute_input":"2023-03-26T19:26:45.244939Z","iopub.status.idle":"2023-03-26T19:26:45.568454Z","shell.execute_reply.started":"2023-03-26T19:26:45.244895Z","shell.execute_reply":"2023-03-26T19:26:45.567053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse.head()\n# sensor_geometry.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-24T17:53:32.527427Z","iopub.execute_input":"2023-03-24T17:53:32.527902Z","iopub.status.idle":"2023-03-24T17:53:32.535134Z","shell.execute_reply.started":"2023-03-24T17:53:32.527858Z","shell.execute_reply":"2023-03-24T17:53:32.533552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse_ = case_event_pulse.merge(\n#     sensor_geometry, on=\"sensor_id\", how=\"left\"\n# ).copy()\n\n# # case_event_pulse = case_event_pulse.reset_index()\n# # case_event_pulse[[\"x\", \"y\", \"z\"]] = sensor_geometry.loc[\n# #     case_event_pulse.sensor_id\n# # ].reset_index()[[\"x\", \"y\", \"z\"]]","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:49.074259Z","iopub.execute_input":"2023-03-26T19:26:49.074662Z","iopub.status.idle":"2023-03-26T19:26:49.091057Z","shell.execute_reply.started":"2023-03-26T19:26:49.074628Z","shell.execute_reply":"2023-03-26T19:26:49.089415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse_.shape, case_event_pulse.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:56.641921Z","iopub.execute_input":"2023-03-26T19:26:56.642368Z","iopub.status.idle":"2023-03-26T19:26:56.652021Z","shell.execute_reply.started":"2023-03-26T19:26:56.642333Z","shell.execute_reply":"2023-03-26T19:26:56.650220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse.x = case_event_pulse.x - case_event_pulse.x.mean()\n# case_event_pulse.y = case_event_pulse.y - case_event_pulse.y.mean()\n# case_event_pulse.z = case_event_pulse.z - case_event_pulse.z.mean()\n# case_event_pulse.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:26:51.241776Z","iopub.execute_input":"2023-03-26T19:26:51.242580Z","iopub.status.idle":"2023-03-26T19:26:51.247015Z","shell.execute_reply.started":"2023-03-26T19:26:51.242531Z","shell.execute_reply":"2023-03-26T19:26:51.245954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# case_event_pulse_.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:27:04.049823Z","iopub.execute_input":"2023-03-26T19:27:04.050275Z","iopub.status.idle":"2023-03-26T19:27:04.067637Z","shell.execute_reply.started":"2023-03-26T19:27:04.050235Z","shell.execute_reply":"2023-03-26T19:27:04.066199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Using Tensorflow LSTM Model Training","metadata":{}},{"cell_type":"code","source":"# Import\nimport numpy as np\nimport os\nimport gc\nimport tensorflow as tf\nimport random\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:28:21.101022Z","iopub.execute_input":"2023-04-02T11:28:21.101598Z","iopub.status.idle":"2023-04-02T11:28:31.919240Z","shell.execute_reply.started":"2023-04-02T11:28:21.101546Z","shell.execute_reply":"2023-04-02T11:28:31.918115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tpu = None\ntry:\n    tpu = tf.distribute.cluster_resolver.TPUClusterResolver()\n    tf.config.experimental_connect_to_cluster(tpu)\n    tf.tpu.experimental.initialize_tpu_system(tpu)\n    strategy = tf.distribute.experimental.TPUStrategy(tpu)\nexcept:\n    strategy = tf.distribute.get_strategy()","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:28:31.921380Z","iopub.execute_input":"2023-04-02T11:28:31.922702Z","iopub.status.idle":"2023-04-02T11:28:31.938256Z","shell.execute_reply.started":"2023-04-02T11:28:31.922658Z","shell.execute_reply":"2023-04-02T11:28:31.937088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Constants","metadata":{}},{"cell_type":"code","source":"# Training\nvalidation_files_amount = 1\ndata_new_load_interval = 6  \ntrain_files_delta = 15  \nepochs = 75  \nbatch_size = 8192  \nlearning_rate = 0.0022  \nverbose = 0\n\n\n# Training Batches\ntrain_batch_id_min = 100\ntrain_batch_id_max = 190\ntrain_batch_ids = list(\n    np.random.permutation(np.arange(train_batch_id_min, train_batch_id_max + 1))\n)\nprint(train_batch_ids)\n\n# Model Parameters\npulse_count = 96\nfeature_count = 6\nlstm_units = 192\nbin_num = 24\n\n# Data\nbase_dir = \"/kaggle/input/lstmicecubesdata/\"\nfile_format = base_dir + \"pp_mpc96_n7_batch_{batch_id:d}.npz\"","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:28:31.939673Z","iopub.execute_input":"2023-04-02T11:28:31.940621Z","iopub.status.idle":"2023-04-02T11:28:31.952966Z","shell.execute_reply.started":"2023-04-02T11:28:31.940576Z","shell.execute_reply":"2023-04-02T11:28:31.951702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set Seed\nseed = 4242\ntf.random.set_seed(seed)\nrandom.seed(seed)\nnp.random.seed(seed)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:28:31.955081Z","iopub.execute_input":"2023-04-02T11:28:31.956099Z","iopub.status.idle":"2023-04-02T11:28:31.960791Z","shell.execute_reply.started":"2023-04-02T11:28:31.956056Z","shell.execute_reply":"2023-04-02T11:28:31.959819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare Metric\n","metadata":{}},{"cell_type":"code","source":"def angular_dist_score(az_true, zen_true, az_pred, zen_pred):\n    \"\"\"\n    Calculate the mean absolute error of the angular distance between two directions.\n    The two vectors are first converted to cartesian unit vectors,\n    and then the Euclidean norm of their difference is computed, which is equal to\n    the sine of half the angle between the two vectors. The inverse\n    sine (arcsin) thereof is then the angle between the two input vectors.\n\n    Parameters:\n    -----------\n    az_true : float or ndarray\n        True azimuth value(s) in radians.\n    zen_true : float or ndarray\n        True zenith value(s) in radians.\n    az_pred : float or ndarray\n        Predicted azimuth value(s) in radians.\n    zen_pred : float or ndarray\n        Predicted zenith value(s) in radians.\n\n    Returns:\n    --------\n    dist : float\n        Mean over the angular distance(s) in radians.\n    \"\"\"\n    az_true, zen_true, az_pred, zen_pred = (\n        np.asarray(az_true),\n        np.asarray(zen_true),\n        np.asarray(az_pred),\n        np.asarray(zen_pred),\n    )\n\n    # Replace NaN and infinity with large finite numbers\n    az_true, zen_true, az_pred, zen_pred = np.nan_to_num(\n        [az_true, zen_true, az_pred, zen_pred], copy=True, nan=1e15, posinf=1e15, neginf=-1e15\n    )\n\n    # Compute cartesian unit vectors\n    true_vec = np.array([np.sin(zen_true) * np.cos(az_true), np.sin(zen_true) * np.sin(az_true), np.cos(zen_true)])\n    pred_vec = np.array([np.sin(zen_pred) * np.cos(az_pred), np.sin(zen_pred) * np.sin(az_pred), np.cos(zen_pred)])\n\n    # Compute sine of half the angle between the vectors\n    sin_half_angle = np.linalg.norm(np.cross(true_vec, pred_vec)) / (np.linalg.norm(true_vec) * np.linalg.norm(pred_vec))\n\n    # Convert back to an angle (in radians)\n    return 2 * np.arcsin(np.clip(sin_half_angle, -1, 1))","metadata":{"execution":{"iopub.status.busy":"2023-04-02T10:46:56.286874Z","iopub.execute_input":"2023-04-02T10:46:56.287450Z","iopub.status.idle":"2023-04-02T10:46:56.300828Z","shell.execute_reply.started":"2023-04-02T10:46:56.287402Z","shell.execute_reply":"2023-04-02T10:46:56.299632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define Azimuth and Zenith Bins¶\n","metadata":{}},{"cell_type":"code","source":"# Create Azimuth Edges\nazimuth_edges = np.linspace(0, 2 * np.pi, bin_num + 1, endpoint=False)\nazimuth_edges[-1] = 2 * np.pi\nprint(azimuth_edges)\n\n# Create Zenith Edges\ncos_zenith_edges = np.linspace(-1, 1, bin_num + 1)\nzenith_edges = np.arccos(cos_zenith_edges)\nprint(zenith_edges)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T10:47:12.570334Z","iopub.execute_input":"2023-04-02T10:47:12.570870Z","iopub.status.idle":"2023-04-02T10:47:12.582463Z","shell.execute_reply.started":"2023-04-02T10:47:12.570821Z","shell.execute_reply":"2023-04-02T10:47:12.580862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Supporting Functions\n","metadata":{}},{"cell_type":"code","source":"angle_bin_zenith0 = zenith_edges[:-1, None]\nangle_bin_zenith1 = zenith_edges[1:, None]\nangle_bin_azimuth0 = azimuth_edges[:-1, None]\nangle_bin_azimuth1 = azimuth_edges[1:, None]\n\narea = (angle_bin_azimuth1 - angle_bin_azimuth0) * (\n    np.cos(angle_bin_zenith0) - np.cos(angle_bin_zenith1)\n)\n\nsin_zenith0 = np.sin(angle_bin_zenith0)\nsin_zenith1 = np.sin(angle_bin_zenith1)\ncos_azimuth0 = np.cos(angle_bin_azimuth0)\ncos_azimuth1 = np.cos(angle_bin_azimuth1)\nsin_azimuth0 = np.sin(angle_bin_azimuth0)\nsin_azimuth1 = np.sin(angle_bin_azimuth1)\nsin_2zenith0 = np.sin(2 * angle_bin_zenith0)\nsin_2zenith1 = np.sin(2 * angle_bin_zenith1)\ncos_2zenith0 = np.cos(2 * angle_bin_zenith0)\ncos_2zenith1 = np.cos(2 * angle_bin_zenith1)\n\nvector_sum_x = (sin_azimuth1 - sin_azimuth0) * (\n    (angle_bin_zenith1 - angle_bin_zenith0) / 2 - (sin_2zenith1 - sin_2zenith0) / 4\n)\nvector_sum_y = (cos_azimuth0 - cos_azimuth1) * (\n    (angle_bin_zenith1 - angle_bin_zenith0) / 2 - (sin_2zenith1 - sin_2zenith0) / 4\n)\nvector_sum_z = (angle_bin_azimuth1 - angle_bin_azimuth0) * (\n    (cos_2zenith0 - cos_2zenith1) / 4\n)\n\narea = np.broadcast_to(area, (bin_num, bin_num, bin_num))\nvector_sum_x = np.broadcast_to(vector_sum_x, (bin_num, bin_num, bin_num))\nvector_sum_y = np.broadcast_to(vector_sum_y, (bin_num, bin_num, bin_num))\nvector_sum_z = np.broadcast_to(vector_sum_z, (bin_num, bin_num, bin_num))\n\nvector_mean_x = np.sum(vector_sum_x, axis=(0, 1)) / np.sum(area, axis=(0, 1))\nvector_mean_y = np.sum(vector_sum_y, axis=(0, 1)) / np.sum(area, axis=(0, 1))\nvector_mean_z = np.sum(vector_sum_z, axis=(0, 1)) / np.sum(area, axis=(0, 1))\n\n\ndef pred_to_angle(pred, epsilon=1e-8):\n    pred_vector = np.einsum('ijk,ik->ij', angle_bin_vector, pred)\n    pred_norm = np.linalg.norm(pred_vector, axis=1)\n    np.clip(pred_norm, epsilon, None, out=pred_norm)\n    pred_vector /= pred_norm[:, None]\n    pred_vector[np.abs(pred_norm) < epsilon] = [1, 0, 0]\n    azimuth = np.arctan2(pred_vector[:, 1], pred_vector[:, 0])\n    azimuth = np.fmod(azimuth + 2*np.pi, 2*np.pi)\n    zenith = np.arccos(pred_vector[:, 2])\n    return azimuth, zenith\n\n\ndef y_to_angle_code(batch_y):\n    azimuth_code = np.digitize(batch_y[:, 0], azimuth_edges[1:])\n    zenith_code = np.digitize(batch_y[:, 1], zenith_edges[1:])\n    angle_code = bin_num * (azimuth_code - 1) + (zenith_code - 1)\n\n    return angle_code","metadata":{"execution":{"iopub.status.busy":"2023-04-02T10:56:48.956047Z","iopub.execute_input":"2023-04-02T10:56:48.957448Z","iopub.status.idle":"2023-04-02T10:56:48.977829Z","shell.execute_reply.started":"2023-04-02T10:56:48.957383Z","shell.execute_reply":"2023-04-02T10:56:48.976465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Loading\n","metadata":{}},{"cell_type":"code","source":"def normalize_data(data):\n    data[:, :, 0] /= 1000  # time\n    data[:, :, 1] /= 300  # charge\n    data[:, :, 3:] /= 600  # space\n    return data\n\n\ndef prep_validation_data(validation_files_amount):\n    print(\"Processing Validation Data...\")\n\n    # Sort train_batch_ids array in ascending order\n    train_batch_ids.sort()\n\n    # Prepare fixed Validation Set\n    val_x, val_y = None, None\n\n    # Loop\n    for batch_id in tqdm(train_batch_ids[:validation_files_amount]):\n        with np.load(file_format.format(batch_id=batch_id)) as val_data_file:\n            x = val_data_file[\"x\"][:, :, [0, 1, 2, 3, 4, 5]]\n            y = val_data_file[\"y\"]\n        val_x = x if val_x is None else np.concatenate([val_x, x], axis=0)\n        val_y = y if val_y is None else np.concatenate([val_y, y], axis=0)\n        del x, y\n\n    # Normalize Data\n    val_x = normalize_data(val_x)\n\n    # Shape Summary\n    print(val_x.shape)\n\n    return val_x, val_y\n\n\ndef prep_training_data(start_batch):\n    print(\"Processing Training Data...\")\n\n    # Placeholders\n    train_x_list, train_y_list = [], []\n\n    # Summary\n    train_ids = sorted(train_batch_ids[start_batch:start_batch + train_files_delta])\n    print(train_ids)\n\n    # Loop\n    for batch_id in tqdm(train_ids):\n        with np.load(file_format.format(batch_id=batch_id)) as train_data_file:\n            train_x = normalize_data(train_data_file[\"x\"][:, :, [0, 1, 2, 3, 4, 5]])\n            train_y = train_data_file[\"y\"]\n            train_y_anglecode = y_to_angle_code(train_y)\n        train_x_list.append(train_x)\n        train_y_list.append(train_y_anglecode)\n        del train_data_file\n\n    # Concatenate arrays\n    train_x = np.concatenate(train_x_list, axis=0)\n    train_y_anglecode = np.concatenate(train_y_list, axis=0)\n\n    # Shape Summary\n    print(train_x.shape)\n\n    return train_x, train_y_anglecode","metadata":{"execution":{"iopub.status.busy":"2023-04-02T10:59:44.292097Z","iopub.execute_input":"2023-04-02T10:59:44.293414Z","iopub.status.idle":"2023-04-02T10:59:44.308921Z","shell.execute_reply.started":"2023-04-02T10:59:44.293344Z","shell.execute_reply":"2023-04-02T10:59:44.307773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prep_validation_data(validation_files_amount)\n# prep_training_data(None)","metadata":{"execution":{"iopub.status.busy":"2023-03-26T19:55:55.375703Z","iopub.execute_input":"2023-03-26T19:55:55.377836Z","iopub.status.idle":"2023-03-26T19:55:55.382911Z","shell.execute_reply.started":"2023-03-26T19:55:55.377686Z","shell.execute_reply":"2023-03-26T19:55:55.380981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras.layers import Input, Masking, Bidirectional, GRU, Dense\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.optimizers import Adam\nimport tensorflow as tf\n\n\ndef create_model(strategy=None):\n    if strategy is not None:\n        with strategy.scope():\n            inputs = Input((pulse_count, feature_count))\n\n            x = Masking(mask_value=0.0, input_shape=(pulse_count, feature_count))(\n                inputs\n            )\n            x = Bidirectional(GRU(lstm_units, return_sequences=True, unroll=True))(x)\n            x = Bidirectional(GRU(lstm_units, return_sequences=True, unroll=True))(x)\n            x = Bidirectional(GRU(lstm_units, unroll=True))(x)\n            x = Dense(256, activation=\"relu\")(x)\n\n            outputs = Dense(bin_num**2, activation=\"softmax\", dtype=\"float32\")(x)\n\n            # Finalize Model\n            model = Model(inputs=inputs, outputs=outputs)\n\n            # Compile model\n            optimizer = Adam(learning_rate=0.001)\n            model.compile(\n                loss=\"sparse_categorical_crossentropy\",\n                optimizer=optimizer,\n                metrics=[\"accuracy\"],\n            )\n\n            # Show Model Summary\n            model.summary()\n\n            return model\n    else:\n        inputs = Input((pulse_count, feature_count))\n\n        x = Masking(mask_value=0.0, input_shape=(pulse_count, feature_count))(inputs)\n        x = Bidirectional(GRU(lstm_units, return_sequences=True, unroll=True))(x)\n        x = Bidirectional(GRU(lstm_units, return_sequences=True, unroll=True))(x)\n        x = Bidirectional(GRU(lstm_units, unroll=True))(x)\n        x = Dense(256, activation=\"relu\")(x)\n\n        outputs = Dense(bin_num**2, activation=\"softmax\", dtype=\"float32\")(x)\n\n        # Finalize Model\n        model = Model(inputs=inputs, outputs=outputs)\n\n        # Compile model\n        optimizer = Adam(learning_rate=0.001)\n        model.compile(\n            loss=\"sparse_categorical_crossentropy\",\n            optimizer=optimizer,\n            metrics=[\"accuracy\"],\n        )\n\n        # Show Model Summary\n        model.summary()\n\n        return model","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:00:51.077476Z","iopub.execute_input":"2023-04-02T11:00:51.078354Z","iopub.status.idle":"2023-04-02T11:00:51.096600Z","shell.execute_reply.started":"2023-04-02T11:00:51.078305Z","shell.execute_reply":"2023-04-02T11:00:51.095485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"create_model(strategy=None)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:00:54.022854Z","iopub.execute_input":"2023-04-02T11:00:54.023592Z","iopub.status.idle":"2023-04-02T11:01:06.333679Z","shell.execute_reply.started":"2023-04-02T11:00:54.023530Z","shell.execute_reply":"2023-04-02T11:01:06.332812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train Model","metadata":{}},{"cell_type":"code","source":"# Define function to load training data\ndef load_train_data(\n    validation_files_amount, data_new_load_interval=None, train_files_delta=None\n):\n    if data_new_load_interval is None and train_files_delta is None:\n        print(\"\\nLoading All Train Data\")\n        start_batch = validation_files_amount\n        end_batch = start_batch + (len(train_batch_ids) - validation_files_amount)\n        trn_x, trn_y_anglecode = prep_training_data(start_batch, end_batch)\n    else:\n        print(f\"\\nLoading Train Data at epoch: {e}\")\n        trn_x, trn_y_anglecode = prep_training_data(validation_files_amount)\n    return trn_x, trn_y_anglecode\n\n\n# Get Fixed Validation Dataset\nval_x, val_y = prep_validation_data(validation_files_amount)\n\n# Create Model\nmodel = create_model()\n\n# Epoch Loop\nfor e in range(epochs):\n    print(f\"=========== EPOCH: {e}\")\n\n    # Load new random batch of training files .. delta wise .. on Kaggle or Colab with limited RAM.\n    if (\n        data_new_load_interval is not None\n        and train_files_delta is not None\n        and e % data_new_load_interval == 0\n    ):\n        if isinstance(train_batch_ids, (list, tuple)) and len(train_batch_ids) > 0:\n            trn_x, trn_y_anglecode = load_train_data(\n                validation_files_amount, data_new_load_interval, train_files_delta\n            )\n        else:\n            print(\"Error: train_batch_ids should be a non-empty list or tuple.\")\n            break\n    else:\n        trn_x, trn_y_anglecode = load_train_data(validation_files_amount)\n\n    # Number of batches\n    batch_count = trn_x.shape[0] // batch_size\n\n    # Random Shuffle each epoch\n    indices = np.arange(trn_x.shape[0])\n    np.random.shuffle(indices)\n    trn_x = trn_x[indices]\n    trn_y_anglecode = trn_y_anglecode[indices]\n\n    # Placeholder\n    losses = []\n    accuracy = []\n\n    # Batch Loop\n    for batch_index in tqdm(range(batch_count), total=batch_count):\n        b_train_x = trn_x[\n            batch_index * batch_size : batch_index * batch_size + batch_size, :\n        ]\n        b_train_y = trn_y_anglecode[\n            batch_index * batch_size : batch_index * batch_size + batch_size\n        ]\n\n        metrics = model.train_on_batch(b_train_x, b_train_y)\n        losses.append(metrics[0])\n        accuracy.append(metrics[1])\n\n    # Save Model\n    model.save(f\"tpu_pp96_n{feature_count}_bin{bin_num}_batch{batch_size}_epoch{e}.h5\")\n\n    # Metrics\n    valid_pred = model.predict(val_x, batch_size=batch_size, verbose=verbose)\n    valid_pred_azimuth, valid_pred_zenith = pred_to_angle(valid_pred)\n    mae = angular_dist_score(\n        val_y[:, 0], val_y[:, 1], valid_pred_azimuth, valid_pred_zenith\n    )\n    print(\n        f\"Total Train Loss: {np.mean(losses):.4f}   Accuracy: {np.mean(accuracy):.4f}  MAE: {mae:.5f}\"\n    )\n\n    # Memory Cleanup\n    gc.collect()\n\n# Save Submission File\nsubmission_df = predict_submission_data()\nsubmission_df.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-04-02T11:12:59.897572Z","iopub.execute_input":"2023-04-02T11:12:59.898101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}