{"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":"# Ensemble of 3 LSTMs\n\n- with Data point picking preprocess and Azimuth-shifted model\n\n---\n\n**ABSTRACT**\n\nThe [*IceCube - Neutrinos in Deep Ice*](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/overview) is a competition to predict the neutrino's incidence direction. The competition hosts provided their simulated data of neutrinos (with some background muons) with true incidence angle, azimuth and zenith. As many participates already observed, a [linear fitting](https://www.kaggle.com/code/shlomoron/icecube-eda-pca-baseline-cv-1-28-lb-1-274) could not get enough score. Muon contaminated events, neutrinos passing detector edge, or curved secondary particles are easily get confused in linear fitting method. And there must be more complicated types of events.\n\nTo handle that many types of events, the LSTM architecture was tried. As the host allowed to [submit with pretrained model](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/discussion/379463), I have splited my code into few notebooks.\n\n1. Data Preprocessing: [LSTM Preprocessing Point Picker](https://www.kaggle.com/code/seungmoklee/lstm-preprocessing-point-picker)\n  - Preprocessed Dataset: [IceCubeData](https://www.kaggle.com/datasets/seungmoklee/icecubedata)\n2. LSTM Training\n  - Initial Training: [LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/lstm-w-gpu-w-npz)\n  - Initial Training with Azimuth Angle Shifted Label: [LSTM AzShift w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/lstm-azshift-w-gpu-w-npz)\n  - Additional Training: [cont' LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/cont-lstm-w-gpu-w-npz)\n  - Additional Training with Azimuth Angle Shifted Label: [cont' AzShift LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/cont-azshift-lstm-w-gpu-w-npz)\n  - Trained models can be found in this dataset: [IceCubeModels](https://www.kaggle.com/datasets/seungmoklee/icecubemodels)\n3. Ensembling models: [LSTM Ensemble - 2 Structure and 2 Phase](https://www.kaggle.com/code/seungmoklee/lstm-ensemble-2-structure-and-2-phase)\n4. Submit: *This one you are reading right now!*\n\n---\n\n*You may have to adjust some hyperparameters to reproduce my score. I've written all the details as far as I can. Please leave comment if you are missing anything, including notebooks, datasets, hyperparameters, settings and explanations. I'd reply.*\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 1. Data Preprocessing\n\n- Notebook link: [LSTM Preprocessing Point Picker](https://www.kaggle.com/code/seungmoklee/lstm-preprocessing-point-picker)\n- Dataset link: [IceCubeData](https://www.kaggle.com/datasets/seungmoklee/icecubedata)\n\nI decided to use maximaly 128 data points for each event.\nEach data points contain its time (relative to the start of each event), charge, auxiliary, position, position resolution (roughly estimated) and rank (to be explained below).\nI assigned the resolution to be the gap of the detector.\nYou can find the relavant figure also in this notebook.\n\nAs some events had too much number of data points, I assigned the rank for each data points considering its importance.\nThe non-auxiliary pulse with highest charge would be the most important one.\nFrom that point, I defined the `valid time window`.\nAs neutrinos travel with speed of light, it must take less than 6200 ns to transverse the detector.\nSo the pulses further than 6200 ns may not be in our interest.\nFrom this intuition, I took non-aux pulses within that `valid time window` first, taking from stronger (high charge) to weaker pulses.\nThe following groups were 2. non-aux pulses out of the `valid time window`, 3. aux pulses in the `valid time window` and 4. aux pulses out of the `valid time window`.\nIf the event had pulses less than 128, zero padding was applied.\n\nI've processed many batches and uploaded them into the dataset linked above.\nI used batches from 400 to 419 for training, and from 101 to 105 for validation and ensemble optimization.\n\nMultiprocessing technique made the code faster about 3 times.\nI learned it from this [notebook](https://www.kaggle.com/code/shlomoron/icecube-eda-pca-baseline-cv-1-28-lb-1-274)\nThe notebook also contained so much intuitions for starting the competetion.\nThanks [greySnow](https://www.kaggle.com/shlomoron).\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 2. LSTM Training\n\n- Notebook links\n  - [LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/lstm-w-gpu-w-npz)\n  - [LSTM AzShift w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/lstm-azshift-w-gpu-w-npz)\n  - [cont' LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/cont-lstm-w-gpu-w-npz)\n  - [cont' AzShift LSTM w/ GPU w/ npz](https://www.kaggle.com/code/seungmoklee/cont-azshift-lstm-w-gpu-w-npz)\n- Dataset link: [IceCubeModels](https://www.kaggle.com/datasets/seungmoklee/icecubemodels)\n\nTraining angle directly is quite ambiguous, as it is cyclic variable.\nDespite training by regression, I trained my models using classification.\nThe idea was adopted from [this paper](http://lps3.doi.org.libproxy.snu.ac.kr/10.1785/0220180311).\nOne-hot encoding labels were gaven to each 2D-angles with fine binning, and trained LSTM models.\nTwo LSTM models were trained, one had 1 LSTM layer with 128 nodes followed by 1 Dense layer with 64 nodes.\nSome additional training was conducted with other data batches, with lower learning rate.\nThe trained models can be found in the dataset above.\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 3. Ensembling models\n\n- Notebook link: [LSTM Ensemble - 2 Structure and 2 Phase](https://www.kaggle.com/code/seungmoklee/lstm-ensemble-2-structure-and-2-phase)\n\nEnsembling various models always helps to squeeze the score.\nI developed optimal combination from this notebook.\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 4. Submit\n\nFinally, submition was done using this notebook!\n\n---","metadata":{}},{"cell_type":"markdown","source":"# START-OF-CODE","metadata":{}},{"cell_type":"code","source":"model_names = [\n    \"PointPicker_mpc128bin16_LSTM128DENSE64_10epc1e-4_10epc1e-5\",\n    \"PointPicker_mpc128bin16_LSTM160DENSE0_10epc1e-4_10epc1e-5\",\n    \"PPAS_mpc128bin16_LSTM128DENSE64_10epc1e-4\",\n]\nmodel_shifted = [False, False, True]\nweights = [0.39, 0.22, 0.39]","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:45:33.445594Z","iopub.execute_input":"2023-02-02T18:45:33.446792Z","iopub.status.idle":"2023-02-02T18:45:33.495951Z","shell.execute_reply.started":"2023-02-02T18:45:33.446622Z","shell.execute_reply":"2023-02-02T18:45:33.494452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set-up\n- Import packages\n- Set hyperparameters","metadata":{}},{"cell_type":"code","source":"# Data I/O and preprocessing\nimport numpy as np\nimport pandas as pd\nimport pyarrow.parquet as pq\n\n# System\nimport time\nimport os\nimport gc\n\n# Graphic\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\n\nfrom tqdm import tqdm\n\n# multiprocessing\nimport multiprocessing\n\n# TENSORFLOW\nimport tensorflow as tf","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-02-02T18:45:33.498369Z","iopub.execute_input":"2023-02-02T18:45:33.499056Z","iopub.status.idle":"2023-02-02T18:45:43.415800Z","shell.execute_reply.started":"2023-02-02T18:45:33.499014Z","shell.execute_reply":"2023-02-02T18:45:43.414504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# directory\nhome_dir = \"/kaggle/input/icecube-neutrinos-in-deep-ice/\"\ntrain_format = home_dir + 'train/batch_{batch_id:d}.parquet'\ntest_format = home_dir + 'test/batch_{batch_id:d}.parquet'\n\nmodel_home = \"/kaggle/input/icecubemodels/\"\n\nweights = np.array(weights)","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:45:43.417390Z","iopub.execute_input":"2023-02-02T18:45:43.418161Z","iopub.status.idle":"2023-02-02T18:45:43.424221Z","shell.execute_reply.started":"2023-02-02T18:45:43.418122Z","shell.execute_reply":"2023-02-02T18:45:43.422913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Model","metadata":{}},{"cell_type":"code","source":"models = list()\nfor model_name in model_names:\n    print(model_name)\n    \n    model_path = model_home + model_name\n    model = tf.keras.models.load_model(model_path)\n    model.summary()\n    \n    models.append(model)","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:45:43.428545Z","iopub.execute_input":"2023-02-02T18:45:43.429030Z","iopub.status.idle":"2023-02-02T18:46:05.653457Z","shell.execute_reply.started":"2023-02-02T18:45:43.428994Z","shell.execute_reply":"2023-02-02T18:46:05.651865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max_pulse_count = model.inputs[0].shape[1]\nn_features = model.inputs[0].shape[2]\noutput_bins = model.layers[-1].weights[0].shape[-1]\n\nbin_num = int(np.sqrt(output_bins))\n\nprint(\"    bin_num    : \", bin_num)\nprint(\"max_pulse_count: \", max_pulse_count)\nprint(\"   n_features  : \", n_features)","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:05.687022Z","iopub.execute_input":"2023-02-02T18:46:05.687538Z","iopub.status.idle":"2023-02-02T18:46:05.696314Z","shell.execute_reply.started":"2023-02-02T18:46:05.687501Z","shell.execute_reply":"2023-02-02T18:46:05.694926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set Detector Geometry","metadata":{}},{"cell_type":"code","source":"%%time\n\n# sensor_geometry\nsensor_geometry_df = pd.read_csv(home_dir + \"sensor_geometry.csv\")\n\n# counts\ndoms_per_string = 60\nstring_num = 86\n\n# index\nouter_long_strings = np.concatenate([np.arange(0, 25), np.arange(27, 34), np.arange(37, 44), np.arange(46, 78)])\ninner_long_strings = np.array([25, 26, 34, 35, 36, 44, 45])\ninner_short_strings = np.array([78, 79, 80, 81, 82, 83, 84, 85])\n\n# known specs\nouter_xy_resolution = 125. / 2\ninner_xy_resolution = 70. / 2\nlong_z_resolution = 17. / 2\nshort_z_resolution = 7. / 2\n\n# evaluate error\nsensor_x = sensor_geometry_df.x\nsensor_y = sensor_geometry_df.y\nsensor_z = sensor_geometry_df.z\nsensor_r_err = np.ones(doms_per_string * string_num)\nsensor_z_err = np.ones(doms_per_string * string_num)\n\nfor string_id in outer_long_strings:\n    sensor_r_err[string_id * doms_per_string:(string_id + 1) * doms_per_string] *= outer_xy_resolution\nfor string_id in np.concatenate([inner_long_strings, inner_short_strings]):\n    sensor_r_err[string_id * doms_per_string:(string_id + 1) * doms_per_string] *= inner_xy_resolution\n\nfor string_id in outer_long_strings:\n    sensor_z_err[string_id * doms_per_string:(string_id + 1) * doms_per_string] *= long_z_resolution\nfor string_id in np.concatenate([inner_long_strings, inner_short_strings]):\n    for dom_id in range(doms_per_string):\n        z = sensor_z[string_id * doms_per_string + dom_id]\n        if (z < -156.) or (z > 95.5 and z < 191.5):\n            sensor_z_err[string_id * doms_per_string + dom_id] *= short_z_resolution\n# register\nsensor_geometry_df[\"r_err\"] = sensor_r_err\nsensor_geometry_df[\"z_err\"] = sensor_z_err","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:05.697439Z","iopub.execute_input":"2023-02-02T18:46:05.697776Z","iopub.status.idle":"2023-02-02T18:46:05.742772Z","shell.execute_reply.started":"2023-02-02T18:46:05.697746Z","shell.execute_reply":"2023-02-02T18:46:05.741735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"theta = np.linspace(0, 2 * np.pi, 361)\n\nfig = plt.figure(figsize=(20, 5))\n\n# 3D plot\nax = fig.add_subplot(141, projection='3d')\n\ns = ax.scatter(sensor_x, sensor_y, sensor_z, s=0.5, c=np.arange(len(sensor_x)), alpha=0.5)\n\nax.set_xlabel(\"X [m]\")\nax.set_ylabel(\"Y [m]\")\nax.set_zlabel(\"Z [m]\")\n\nfig.colorbar(s, ax=ax)\n\n# X-Y plot\nax = fig.add_subplot(142)\n\nfor string_id in outer_long_strings:\n    x = sensor_x[string_id * doms_per_string]\n    y = sensor_y[string_id * doms_per_string]\n    r_err = sensor_r_err[string_id * doms_per_string]\n    scatter_outer_long = ax.scatter(x, y, color=\"blue\", label=\"outer long string\")\n    ax.plot(x + r_err * np.cos(theta), y + r_err * np.sin(theta), color=\"gray\", alpha=0.5)\n    \nfor string_id in inner_long_strings:\n    x = sensor_x[string_id * doms_per_string]\n    y = sensor_y[string_id * doms_per_string]\n    r_err = sensor_r_err[string_id * doms_per_string]\n    scatter_inner_long = ax.scatter(x, y, color=\"orange\", label=\"inner long string\")\n    ax.plot(x + r_err * np.cos(theta), y + r_err * np.sin(theta), color=\"gray\", alpha=0.5)\n    \nfor string_id in inner_short_strings:\n    x = sensor_x[string_id * doms_per_string]\n    y = sensor_y[string_id * doms_per_string]\n    r_err = sensor_r_err[string_id * doms_per_string]\n    scatter_inner_short = ax.scatter(x, y, color=\"red\", label=\"inner short string\")\n    ax.plot(x + r_err * np.cos(theta), y + r_err * np.sin(theta), color=\"gray\", alpha=0.5)\n\nax.set_xlabel(\"X [m]\")\nax.set_ylabel(\"Y [m]\")\nax.legend(handles=[scatter_outer_long, scatter_inner_long, scatter_inner_short])\n\n# X-Z plot\nax = fig.add_subplot(143)\n\nfor string_id in outer_long_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_outer_long = ax.scatter(x, z, s=0.5, color=\"blue\", label=\"outer long string\")\n    \nfor string_id in inner_long_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_inner_long = ax.scatter(x, z, s=0.5, color=\"orange\", label=\"inner long string\")\n\nfor string_id in inner_short_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_inner_short = ax.scatter(x, z, s=0.5, color=\"red\", label=\"inner short string\")\n\nax.set_xlabel(\"X [m]\")\nax.set_ylabel(\"Z [m]\")\nax.legend(handles=[scatter_outer_long, scatter_inner_long, scatter_inner_short])\n\n# X-Z plot zoom\nax = fig.add_subplot(144)\n\nfor string_id in outer_long_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_outer_long = ax.scatter(x, z, s=0.5, color=\"blue\", label=\"outer long string\")\n    \nfor string_id in inner_long_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_inner_long = ax.scatter(x, z, s=0.5, color=\"orange\", label=\"inner long string\")\n\nfor string_id in inner_short_strings:\n    x = sensor_x[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    z = sensor_z[string_id * doms_per_string:string_id * doms_per_string + doms_per_string]\n    scatter_inner_short = ax.scatter(x, z, s=0.5, color=\"red\", label=\"inner short string\")\n\nfor sensor_id in range(doms_per_string * string_num):\n    x = sensor_x[sensor_id]\n    z = sensor_z[sensor_id]\n    z_err = sensor_z_err[sensor_id]\n    if (x > -150 and x < 50) and (z > -200 and z < 250):\n        ax.plot(x + z_err * np.cos(theta), z + z_err * np.sin(theta), color=\"gray\", alpha=0.5)\n\nax.set_xlabel(\"X [m]\")\nax.set_ylabel(\"Z [m]\")\nax.set_xlim(-150, 50)\nax.set_ylim(-200, 250)\nax.legend(handles=[scatter_outer_long, scatter_inner_long, scatter_inner_short])\n\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:05.744302Z","iopub.execute_input":"2023-02-02T18:46:05.744910Z","iopub.status.idle":"2023-02-02T18:46:14.257114Z","shell.execute_reply.started":"2023-02-02T18:46:05.744863Z","shell.execute_reply":"2023-02-02T18:46:14.255938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# detector constants\nc_const = 0.299792458  # speed of light [m/ns]\n\nx_min = sensor_x.min()\nx_max = sensor_x.max()\ny_min = sensor_y.min()\ny_max = sensor_y.max()\nz_min = sensor_z.min()\nz_max = sensor_z.max()\n\ndetector_length = np.sqrt((x_max - x_min)**2 + (y_max - y_min)**2 + (z_max - z_min)**2)\nt_valid_length = detector_length / c_const\n\nprint(\"t_valid_length: \", t_valid_length, \" ns\")","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:14.258778Z","iopub.execute_input":"2023-02-02T18:46:14.259513Z","iopub.status.idle":"2023-02-02T18:46:14.269268Z","shell.execute_reply.started":"2023-02-02T18:46:14.259467Z","shell.execute_reply":"2023-02-02T18:46:14.268014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data I/O Helper","metadata":{}},{"cell_type":"markdown","source":"## Angle one-hot encoding edges\n\n- It is efficient to train the model by classification task, initially.\n- azimuth and zenith are independent\n- azimuth distribution is flat and zenith distribution is sinusoidal.\n  - Flat on the spherical surface\n  - $\\phi > \\pi$ events are a little bit rarer than $\\phi < \\pi$ events, (maybe) because of the neutrino attenuation by earth.\n- So, the uniform bin is used for azimuth, and $\\left| \\cos \\right|$ bin is used for zenith","metadata":{}},{"cell_type":"code","source":"azimuth_edges = np.linspace(0, 2 * np.pi, bin_num + 1)\nazimuth_shift = (azimuth_edges[1] - azimuth_edges[0]) / 2.\n\nzenith_edges_flat = np.linspace(0, np.pi, bin_num + 1)\nzenith_edges = list()\nzenith_edges.append(0)\nfor bin_idx in range(1, bin_num):\n    # cos(zen_before) - cos(zen_now) = 2 / bin_num\n    zen_now = np.arccos(np.cos(zenith_edges[-1]) - 2 / (bin_num))\n    zenith_edges.append(zen_now)\nzenith_edges.append(np.pi)\nzenith_edges = np.array(zenith_edges)","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:14.271938Z","iopub.execute_input":"2023-02-02T18:46:14.272643Z","iopub.status.idle":"2023-02-02T18:46:14.328701Z","shell.execute_reply.started":"2023-02-02T18:46:14.272597Z","shell.execute_reply":"2023-02-02T18:46:14.327026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def y_to_onehot(batch_y):\n    # evaluate bin code\n    azimuth_code = (batch_y[:, 0] > azimuth_edges[1:].reshape((-1, 1))).sum(axis=0)\n    zenith_code = (batch_y[:, 1] > zenith_edges[1:].reshape((-1, 1))).sum(axis=0)\n    angle_code = bin_num * azimuth_code + zenith_code\n\n    # one-hot\n    batch_y_onehot = np.zeros((angle_code.size, bin_num * bin_num))\n    batch_y_onehot[np.arange(angle_code.size), angle_code] = 1\n    \n    return batch_y_onehot","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:14.330625Z","iopub.execute_input":"2023-02-02T18:46:14.332351Z","iopub.status.idle":"2023-02-02T18:46:14.344665Z","shell.execute_reply.started":"2023-02-02T18:46:14.332287Z","shell.execute_reply":"2023-02-02T18:46:14.343365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define a function converts from prediction to angles\n\n- Calculation of the mean-vector in a bin $\\theta \\in ( \\theta_0, \\theta_1 )$ and $\\phi \\in ( \\phi_0, \\phi_1 )$\n  - $\\vec{r} \\left( \\theta, ~ \\phi \\right) = \\left< \\sin \\theta \\cos \\phi, ~ \\sin \\theta \\sin \\phi, ~ \\cos \\theta \\right>$\n  - $\\bar{\\vec{r}} = \\frac{ \\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} \\vec{r} \\left( \\theta, ~ \\phi \\right) \\sin \\theta \\,d\\phi \\,d\\theta }{ \\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} 1 \\sin \\theta \\,d\\phi \\,d\\theta }$\n  - $ \\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} 1 \\sin \\theta \\,d\\phi \\,d\\theta = \\left( \\phi_1 - \\phi_0 \\right) \\left( \\cos \\theta_0 - \\cos \\theta_1 \\right)$\n  - $\n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} {r}_{x} \\left( \\theta, ~ \\phi \\right) \\sin \\theta \\,d\\phi \\,d\\theta = \n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} \\sin^2 \\theta \\cos \\phi \\,d\\phi \\,d\\theta = \n\\left( \\sin \\phi_1 - \\sin \\phi_0 \\right) \\left( \\frac{\\theta_1 - \\theta_0}{2} - \\frac{\\sin 2 \\theta_1 - \\sin 2 \\theta_0}{4} \\right)\n$\n  - $\n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} {r}_{y} \\left( \\theta, ~ \\phi \\right) \\sin \\theta \\,d\\phi \\,d\\theta = \n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} \\sin^2 \\theta \\sin \\phi \\,d\\phi \\,d\\theta = \n\\left( \\cos \\phi_0 - \\cos \\phi_1 \\right) \\left( \\frac{\\theta_1 - \\theta_0}{2} - \\frac{\\sin 2 \\theta_1 - \\sin 2 \\theta_0}{4} \\right)\n$\n  - $\n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} {r}_{z} \\left( \\theta, ~ \\phi \\right) \\sin \\theta \\,d\\phi \\,d\\theta = \n\\int_{\\theta_{0}}^{\\theta_{1}} \\int_{\\phi_0}^{\\phi_1} \\sin \\theta \\cos \\theta \\,d\\phi \\,d\\theta = \n\\left( \\phi_1 - \\phi_0 \\right) \\left( \\frac{\\cos 2 \\theta_0 - \\cos 2 \\theta_1}{4} \\right)\n$","metadata":{}},{"cell_type":"code","source":"angle_bin_zenith0 = np.tile(zenith_edges[:-1], bin_num)\nangle_bin_zenith1 = np.tile(zenith_edges[1:], bin_num)\nangle_bin_azimuth0 = np.repeat(azimuth_edges[:-1], bin_num)\nangle_bin_azimuth1 = np.repeat(azimuth_edges[1:], bin_num)\n\nangle_bin_area = (angle_bin_azimuth1 - angle_bin_azimuth0) * (np.cos(angle_bin_zenith0) - np.cos(angle_bin_zenith1))\nangle_bin_vector_sum_x = (np.sin(angle_bin_azimuth1) - np.sin(angle_bin_azimuth0)) * ((angle_bin_zenith1 - angle_bin_zenith0) / 2 - (np.sin(2 * angle_bin_zenith1) - np.sin(2 * angle_bin_zenith0)) / 4)\nangle_bin_vector_sum_y = (np.cos(angle_bin_azimuth0) - np.cos(angle_bin_azimuth1)) * ((angle_bin_zenith1 - angle_bin_zenith0) / 2 - (np.sin(2 * angle_bin_zenith1) - np.sin(2 * angle_bin_zenith0)) / 4)\nangle_bin_vector_sum_z = (angle_bin_azimuth1 - angle_bin_azimuth0) * ((np.cos(2 * angle_bin_zenith0) - np.cos(2 * angle_bin_zenith1)) / 4)\n\nangle_bin_vector_mean_x = angle_bin_vector_sum_x / angle_bin_area\nangle_bin_vector_mean_y = angle_bin_vector_sum_y / angle_bin_area\nangle_bin_vector_mean_z = angle_bin_vector_sum_z / angle_bin_area\n\nangle_bin_vector = np.zeros((1, bin_num * bin_num, 3))\nangle_bin_vector[:, :, 0] = angle_bin_vector_mean_x\nangle_bin_vector[:, :, 1] = angle_bin_vector_mean_y\nangle_bin_vector[:, :, 2] = angle_bin_vector_mean_z\n\nangle_bin_vector_unit = angle_bin_vector[0].copy()\nangle_bin_vector_unit /= np.sqrt((angle_bin_vector_unit**2).sum(axis=1).reshape((-1, 1)))","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:29.121456Z","iopub.execute_input":"2023-02-02T18:46:29.121979Z","iopub.status.idle":"2023-02-02T18:46:29.135632Z","shell.execute_reply.started":"2023-02-02T18:46:29.121940Z","shell.execute_reply":"2023-02-02T18:46:29.134259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def pred_to_angle(pred, epsilon=1e-8):\n    # convert prediction to vector\n    pred_vector = (pred.reshape((-1, bin_num * bin_num, 1)) * angle_bin_vector).sum(axis=1)\n    \n    # normalize\n    pred_vector_norm = np.sqrt((pred_vector**2).sum(axis=1))\n    mask = pred_vector_norm < epsilon\n    pred_vector_norm[mask] = 1\n    \n    # assign <1, 0, 0> to very small vectors (badly predicted)\n    pred_vector /= pred_vector_norm.reshape((-1, 1))\n    pred_vector[mask] = np.array([1., 0., 0.])\n    \n    # convert to angle\n    azimuth = np.arctan2(pred_vector[:, 1], pred_vector[:, 0])\n    azimuth[azimuth < 0] += 2 * np.pi\n    zenith = np.arccos(pred_vector[:, 2])\n    \n    # mask bad norm predictions as 0, 0\n    azimuth[mask] = 0.\n    zenith[mask] = 0.\n    \n    return azimuth, zenith\n\n\ndef pred_to_angle_azshift(pred, epsilon=1e-8):\n    # convert prediction to vector\n    pred_vector = (pred.reshape((-1, bin_num * bin_num, 1)) * angle_bin_vector).sum(axis=1)\n    \n    # normalize\n    pred_vector_norm = np.sqrt((pred_vector**2).sum(axis=1))\n    mask = pred_vector_norm < epsilon\n    pred_vector_norm[mask] = 1\n    \n    # assign <1, 0, 0> to very small vectors (badly predicted)\n    pred_vector /= pred_vector_norm.reshape((-1, 1))\n    pred_vector[mask] = np.array([1., 0., 0.])\n    \n    # convert to angle\n    azimuth = np.arctan2(pred_vector[:, 1], pred_vector[:, 0])\n    azimuth[azimuth < 0] += 2 * np.pi\n    zenith = np.arccos(pred_vector[:, 2])\n    \n    # shift\n    azimuth -= azimuth_shift\n    azimuth[azimuth < 0] += 2 * np.pi\n    \n    # mask bad norm predictions as 0, 0\n    azimuth[mask] = 0.\n    zenith[mask] = 0.\n    \n    return azimuth, zenith","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:29.907715Z","iopub.execute_input":"2023-02-02T18:46:29.908610Z","iopub.status.idle":"2023-02-02T18:46:29.923123Z","shell.execute_reply.started":"2023-02-02T18:46:29.908559Z","shell.execute_reply":"2023-02-02T18:46:29.921676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def pred_to_angle_argmax(pred):\n    # get the highest score codes\n    pred_code = pred.argmax(axis=1)\n    \n    # get the bin vector\n    pred_vector = angle_bin_vector_unit[pred_code, :]\n    \n    # convert to angle\n    azimuth = np.arctan2(pred_vector[:, 1], pred_vector[:, 0])\n    azimuth[azimuth < 0] += 2 * np.pi\n    zenith = np.arccos(pred_vector[:, 2])\n    \n    return azimuth, zenith","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:30.973835Z","iopub.execute_input":"2023-02-02T18:46:30.975081Z","iopub.status.idle":"2023-02-02T18:46:30.984238Z","shell.execute_reply.started":"2023-02-02T18:46:30.975021Z","shell.execute_reply":"2023-02-02T18:46:30.982613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Weighted-Vector Ensemble","metadata":{}},{"cell_type":"code","source":"def weighted_vector_ensemble(angles, weight):\n    # Convert angle to vector\n    vec_models = list()\n    for angle in angles:\n        az, zen = angle\n\n        sa = np.sin(az)\n        ca = np.cos(az)\n        sz = np.sin(zen)\n        cz = np.cos(zen)\n\n        vec = np.stack([sz * ca, sz * sa, cz], axis=1)\n        vec_models.append(vec)\n    vec_models = np.array(vec_models)\n\n    # Weighted-mean\n    vec_mean = (weight.reshape((-1, 1, 1)) * vec_models).sum(axis=0) / weight.sum()\n    vec_mean /= np.sqrt((vec_mean**2).sum(axis=1)).reshape((-1, 1))\n\n    # Convert vector to angle\n    zenith = np.arccos(vec_mean[:, 2])\n    azimuth = np.arctan2(vec_mean[:, 1], vec_mean[:, 0])\n    azimuth[azimuth < 0] += 2 * np.pi\n    \n    return azimuth, zenith","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:50:41.983976Z","iopub.execute_input":"2023-02-02T18:50:41.984527Z","iopub.status.idle":"2023-02-02T18:50:41.994966Z","shell.execute_reply.started":"2023-02-02T18:50:41.984487Z","shell.execute_reply":"2023-02-02T18:50:41.993853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Single event reader function\n\n- Pick-up important data points first\n    - Rank 3 (First)\n        - not aux, in valid time window\n    - Rank 2\n        - not aux, out of valid time window\n    - Rank 1\n        - aux, in valid time window\n    - Rank 0 (Last)\n        - aux, out of valid time window\n    - In each ranks, take pulses from highest charge","metadata":{}},{"cell_type":"code","source":"open_batch_dict = dict()\n\n\n# read single event from batch_meta_df\ndef read_event(event_idx, batch_meta_df, max_pulse_count, train=True):\n    # read metadata\n    batch_id, first_pulse_index, last_pulse_index = batch_meta_df.iloc[event_idx][[\"batch_id\", \"first_pulse_index\", \"last_pulse_index\"]].astype(\"int\")\n\n    # close past batch df\n    if batch_id - 1 in open_batch_dict.keys():\n        del open_batch_dict[batch_id - 1]\n\n    # open current batch df\n    if batch_id not in open_batch_dict.keys():\n        if train:\n            open_batch_dict.update({batch_id: pd.read_parquet(train_format.format(batch_id=batch_id))})\n        else:\n            open_batch_dict.update({batch_id: pd.read_parquet(test_format.format(batch_id=batch_id))})\n    \n    batch_df = open_batch_dict[batch_id]\n    \n    # read event\n    event_feature = batch_df[first_pulse_index:last_pulse_index + 1]\n    sensor_id = event_feature.sensor_id\n    \n    # merge features into single structured array\n    dtype = [\n        (\"time\", \"float16\"),\n        (\"charge\", \"float16\"),\n        (\"auxiliary\", \"float16\"),\n        (\"x\", \"float16\"),\n        (\"y\", \"float16\"),\n        (\"z\", \"float16\"),\n        (\"r_err\", \"float16\"),\n        (\"z_err\", \"float16\"),\n        (\"rank\", \"short\"),\n    ]\n    event_x = np.zeros(last_pulse_index - first_pulse_index + 1, dtype)\n\n    event_x[\"time\"] = event_feature.time.values - event_feature.time.min()\n    event_x[\"charge\"] = event_feature.charge.values\n    event_x[\"auxiliary\"] = event_feature.auxiliary.values\n\n    event_x[\"x\"] = sensor_geometry_df.x[sensor_id].values\n    event_x[\"y\"] = sensor_geometry_df.y[sensor_id].values\n    event_x[\"z\"] = sensor_geometry_df.z[sensor_id].values\n\n    event_x[\"r_err\"] = sensor_geometry_df.r_err[sensor_id].values\n    event_x[\"z_err\"] = sensor_geometry_df.z_err[sensor_id].values\n    \n    # For long event, pick-up\n    if len(event_x) > max_pulse_count:\n        # Find valid time window\n        t_peak = event_x[\"time\"][event_x[\"charge\"].argmax()]\n        t_valid_min = t_peak - t_valid_length\n        t_valid_max = t_peak + t_valid_length\n\n        t_valid = (event_x[\"time\"] > t_valid_min) * (event_x[\"time\"] < t_valid_max)\n\n        # rank\n        event_x[\"rank\"] = 2 * (1 - event_x[\"auxiliary\"]) + (t_valid)\n\n        # sort by rank and charge (important goes to backward)\n        event_x = np.sort(event_x, order=[\"rank\", \"charge\"])\n\n        # pick-up from backward\n        event_x = event_x[-max_pulse_count:]\n\n        # resort by time\n        event_x = np.sort(event_x, order=\"time\")\n\n    # for train data, give angles together\n    if train:\n        azimuth, zenith = batch_meta_df.iloc[event_idx][[\"azimuth\", \"zenith\"]].astype(\"float16\")\n        event_y = np.array([azimuth, zenith], dtype=\"float16\")\n        \n        return event_idx, len(event_x), event_x, event_y\n    \n    # for test data, just give feature \n    else:\n        return event_idx, len(event_x), event_x","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:36.298510Z","iopub.execute_input":"2023-02-02T18:46:36.299113Z","iopub.status.idle":"2023-02-02T18:46:36.321774Z","shell.execute_reply.started":"2023-02-02T18:46:36.299063Z","shell.execute_reply":"2023-02-02T18:46:36.320118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data I/O (for CPU) & Normalization\n\n- Read data\n- Concatenate and convert\n- Normalize time, charge and position variables","metadata":{}},{"cell_type":"markdown","source":"## Read test metadata and define spliter (for CPU)","metadata":{}},{"cell_type":"code","source":"test_meta_df = pq.read_table(home_dir + 'test_meta.parquet').to_pandas()\ntest_meta_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:38.772817Z","iopub.execute_input":"2023-02-02T18:46:38.773329Z","iopub.status.idle":"2023-02-02T18:46:38.913920Z","shell.execute_reply.started":"2023-02-02T18:46:38.773288Z","shell.execute_reply":"2023-02-02T18:46:38.912675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_counts = test_meta_df.batch_id.value_counts().sort_index()\n\nbatch_max_index = batch_counts.cumsum()\nbatch_max_index[test_meta_df.batch_id.min() - 1] = 0\nbatch_max_index = batch_max_index.sort_index()\n\n\ndef test_meta_df_spliter(batch_id):\n    return test_meta_df.loc[batch_max_index[batch_id - 1]:batch_max_index[batch_id] - 1]","metadata":{"execution":{"iopub.status.busy":"2023-02-02T18:46:39.897196Z","iopub.execute_input":"2023-02-02T18:46:39.897671Z","iopub.status.idle":"2023-02-02T18:46:39.911231Z","shell.execute_reply.started":"2023-02-02T18:46:39.897630Z","shell.execute_reply":"2023-02-02T18:46:39.909409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Read test data and predict batch-by-batch","metadata":{}},{"cell_type":"code","source":"test_batch_ids = test_meta_df.batch_id.unique()\n\ntest_event_id = list()\ntest_azimuth = list()\ntest_zenith = list()\n\nfor batch_id in test_batch_ids:\n    print(batch_id)\n    # READ ONE BATCH OF TEST DATA\n    # get batch meta data\n    batch_meta_df = test_meta_df_spliter(batch_id)\n\n    # register pulses\n    test_x = np.zeros((len(batch_meta_df), max_pulse_count, n_features), dtype=\"float16\")    \n    test_x[:, :, 2] = -1\n    \n\n    def read_event_local(event_idx):\n        return read_event(event_idx, batch_meta_df, max_pulse_count, train=False)\n\n    \n    # scan events\n    iterator = range(len(batch_meta_df))\n    with multiprocessing.Pool() as pool:\n        for event_idx, pulse_count, event_x in pool.map(read_event_local, iterator):\n            # feature\n            test_x[event_idx, :pulse_count, 0] = event_x[\"time\"]\n            test_x[event_idx, :pulse_count, 1] = event_x[\"charge\"]\n            test_x[event_idx, :pulse_count, 2] = event_x[\"auxiliary\"]\n            test_x[event_idx, :pulse_count, 3] = event_x[\"x\"]\n            test_x[event_idx, :pulse_count, 4] = event_x[\"y\"]\n            test_x[event_idx, :pulse_count, 5] = event_x[\"z\"]\n            test_x[event_idx, :pulse_count, 6] = event_x[\"r_err\"]\n            test_x[event_idx, :pulse_count, 7] = event_x[\"z_err\"]\n    \n    del batch_meta_df\n    \n    # CONVERT\n    test_x[:, :, 0] /= 1000  # time\n    test_x[:, :, 1] /= 300  # charge\n    test_x[:, :, 3:] /= 600  # space\n    \n    # PREDICT\n    pred_angles = list()\n    for model, shifted in zip(models, model_shifted):\n        pred_model = model.predict(test_x, verbose=0)\n        \n        if shifted:\n            az_model, zen_model = pred_to_angle_azshift(pred_model)\n        else:\n            az_model, zen_model = pred_to_angle(pred_model)\n    \n        pred_angles.append((az_model, zen_model))\n    \n    pred_azimuth, pred_zenith = weighted_vector_ensemble(pred_angles, weights)\n    \n    event_ids = test_meta_df.event_id[test_meta_df.batch_id == batch_id].values\n    \n    for event_id, azimuth, zenith in zip(event_ids, pred_azimuth, pred_zenith):\n        if np.isfinite(azimuth) and np.isfinite(zenith):\n            test_event_id.append(int(event_id))\n            test_azimuth.append(azimuth)\n            test_zenith.append(zenith)\n        else:\n            test_event_id.append(int(event_id))\n            test_azimuth.append(0.)\n            test_zenith.append(0.)","metadata":{"execution":{"iopub.status.busy":"2023-02-01T17:02:06.971788Z","iopub.execute_input":"2023-02-01T17:02:06.972919Z","iopub.status.idle":"2023-02-01T17:02:08.301916Z","shell.execute_reply.started":"2023-02-01T17:02:06.972861Z","shell.execute_reply":"2023-02-01T17:02:08.300346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submit","metadata":{}},{"cell_type":"code","source":"test_result_dict = {\n    \"event_id\": test_event_id,\n    \"azimuth\": test_azimuth,\n    \"zenith\": test_zenith,\n}\n\ntest_result_df = pd.DataFrame(test_result_dict)\ntest_result_df = test_result_df.sort_values(by=['event_id'])\n\ntest_result_df.to_csv(\"submission.csv\", index=False)\ntest_result_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-01T17:02:08.303614Z","iopub.execute_input":"2023-02-01T17:02:08.304002Z","iopub.status.idle":"2023-02-01T17:02:08.327245Z","shell.execute_reply.started":"2023-02-01T17:02:08.303968Z","shell.execute_reply":"2023-02-01T17:02:08.325899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# END-OF-NOTE","metadata":{}}]}