{"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":"code","source":"# Import Modules\nimport gc\nimport os\nimport multiprocessing\nimport time\nimport numpy as np\nimport pandas as pd\nimport pyarrow.parquet as pq\nimport tensorflow as tf\nfrom tqdm import tqdm","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-20T16:40:03.835155Z","iopub.execute_input":"2023-04-20T16:40:03.836212Z","iopub.status.idle":"2023-04-20T16:40:03.842065Z","shell.execute_reply.started":"2023-04-20T16:40:03.836172Z","shell.execute_reply":"2023-04-20T16:40:03.840979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"home_dir = \"/kaggle/input/icecube-neutrinos-in-deep-ice/\"\n#test_meta_df = pd.read_parquet('/kaggle/input/meta-neutrinos/meta_'+str(353)).reset_index()\ntest_meta_df = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet').reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:03.846765Z","iopub.execute_input":"2023-04-20T16:40:03.847042Z","iopub.status.idle":"2023-04-20T16:40:03.933699Z","shell.execute_reply.started":"2023-04-20T16:40:03.847016Z","shell.execute_reply":"2023-04-20T16:40:03.932716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Directories and constants\nhome_dir = \"/kaggle/input/icecube-neutrinos-in-deep-ice/\"\n#test_format = home_dir + 'train/batch_{batch_id:d}.parquet'\ntest_format = home_dir + 'test/batch_{batch_id:d}.parquet'\nmodel_home = \"/kaggle/input/lstmicecubesdata/\"\n\n# Model(s)\nmodel_names = [\"4347_MAE_1-02076_bin24_pp96_n6_batch2048_epoch29.h5\",\n               \"4347_MAE_1-02039_bin24_pp96_n6_batch2048_epoch25.h5\", \n               \"4346_MAE_1-02020_bin24_pp96_n6_batch2048_epoch27.h5\"]\nmodel_weights = np.array([0.30, \n                          0.30,\n                          0.40])","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:03.935873Z","iopub.execute_input":"2023-04-20T16:40:03.936241Z","iopub.status.idle":"2023-04-20T16:40:03.941696Z","shell.execute_reply.started":"2023-04-20T16:40:03.936203Z","shell.execute_reply":"2023-04-20T16:40:03.940632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Models\nmodels = []\nfor model_name in model_names:\n    print(f'\\n========== Model File: {model_name}')\n    \n    # Load Model\n    model_path = model_home + model_name\n    model = tf.keras.models.load_model(model_path)\n    models.append(model)      \n    \n    # Model summary\n    model.summary()\n    \n# Get Model Parameters\npulse_count = model.inputs[0].shape[1]\nfeature_count = model.inputs[0].shape[2]\noutput_bins = model.layers[-1].weights[0].shape[-1]\nbin_num = int(np.sqrt(output_bins))\n\n# Model Parameter Summary\nprint(\"\\n==== Model Parameters\")\nprint(f\"Bin Numbers: {bin_num}\")\nprint(f\"Maximum Pulse Count: {pulse_count}\")\nprint(f\"Features Count: {feature_count}\")","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:03.943409Z","iopub.execute_input":"2023-04-20T16:40:03.944127Z","iopub.status.idle":"2023-04-20T16:40:26.919047Z","shell.execute_reply.started":"2023-04-20T16:40:03.944090Z","shell.execute_reply":"2023-04-20T16:40:26.918193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()\n# Load sensor_geometry\nsensor_geometry_df = pd.read_csv(home_dir + \"sensor_geometry.csv\")\n\n# Get Sensor Information\nsensor_x = sensor_geometry_df.x\nsensor_y = sensor_geometry_df.y\nsensor_z = sensor_geometry_df.z\n\n# Detector constants\nc_const = 0.299792458  # speed of light [m/ns]\n\n# Sensor Min / Max Coordinates\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(f\"time valid length: {t_valid_length} ns\")","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:26.921743Z","iopub.execute_input":"2023-04-20T16:40:26.922107Z","iopub.status.idle":"2023-04-20T16:40:27.177587Z","shell.execute_reply.started":"2023-04-20T16:40:26.922067Z","shell.execute_reply":"2023-04-20T16:40:27.176484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create Azimuth Edges\nazimuth_edges = np.linspace(0, 2 * np.pi, bin_num + 1)\nprint(azimuth_edges)\n\n# Create Zenith Edges\nzenith_edges = []\nzenith_edges.append(0)\nfor bin_idx in range(1, bin_num):\n    zenith_edges.append(np.arccos(np.cos(zenith_edges[-1]) - 2 / (bin_num)))\nzenith_edges.append(np.pi)\nzenith_edges = np.array(zenith_edges)\nprint(zenith_edges)","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:27.179530Z","iopub.execute_input":"2023-04-20T16:40:27.179907Z","iopub.status.idle":"2023-04-20T16:40:27.192851Z","shell.execute_reply.started":"2023-04-20T16:40:27.179867Z","shell.execute_reply":"2023-04-20T16:40:27.191292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-04-20T16:40:27.194872Z","iopub.execute_input":"2023-04-20T16:40:27.195535Z","iopub.status.idle":"2023-04-20T16:40:27.208585Z","shell.execute_reply.started":"2023-04-20T16:40:27.195497Z","shell.execute_reply":"2023-04-20T16:40:27.207423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Placeholder\nopen_batch_dict = dict()\n\n# Read single event from batch_meta_df\ndef read_event(event_idx, batch_meta_df, pulse_count):\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        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 = [(\"time\", \"float16\"),\n             (\"charge\", \"float16\"),\n             (\"auxiliary\", \"float16\"),\n             (\"x\", \"float16\"),\n             (\"y\", \"float16\"),\n             (\"z\", \"float16\"),\n             (\"rank\", \"short\")]    \n    \n    # Create event_x\n    event_x = np.zeros(last_pulse_index - first_pulse_index + 1, dtype)\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    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    # For long event, pick-up\n    if len(event_x) > 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        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[-pulse_count:]\n\n        # Sort events by time \n        event_x = np.sort(event_x, order = \"time\")\n\n    return event_idx, len(event_x), event_x","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:40:27.210463Z","iopub.execute_input":"2023-04-20T16:40:27.210886Z","iopub.status.idle":"2023-04-20T16:40:27.225321Z","shell.execute_reply.started":"2023-04-20T16:40:27.210850Z","shell.execute_reply":"2023-04-20T16:40:27.224265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read Test Meta data\n#test_meta_df = (pq.read_table(home_dir + 'test_meta.parquet').to_pandas())\n#test_meta_df = pd.read_parquet(home_dir + 'train_meta.parquet')[0:100]\n\ngc.collect()\n\nbatch_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# Support Function\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-04-20T16:40:27.227307Z","iopub.execute_input":"2023-04-20T16:40:27.228544Z","iopub.status.idle":"2023-04-20T16:40:27.414190Z","shell.execute_reply.started":"2023-04-20T16:40:27.228502Z","shell.execute_reply":"2023-04-20T16:40:27.412974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get Batch IDs\ntest_batch_ids = test_meta_df.batch_id.unique()\n\n# Submission Placeholders\ntest_event_id = []\ntest_azimuth = []\ntest_zenith = []\n\npred_raw_0=[]\n\n# Batch Loop\nfor batch_id in test_batch_ids:\n    # Batch Meta DF\n    batch_meta_df = test_meta_df_spliter(batch_id)\n\n    # Set Pulses\n    test_x = np.zeros((len(batch_meta_df), pulse_count, feature_count), dtype = \"float16\")    \n    test_x[:, :, 2] = -1    \n\n    # Read Event Data\n    def read_event_local(event_idx):\n        return read_event(event_idx, batch_meta_df, pulse_count)\n    \n    # Multiprocess Events\n    iterator = range(len(batch_meta_df))\n    with multiprocessing.Pool() as pool:\n        for event_idx, pulsecount, event_x in pool.map(read_event_local, iterator):\n            # Features\n            test_x[event_idx, :pulsecount, 0] = event_x[\"time\"]\n            test_x[event_idx, :pulsecount, 1] = event_x[\"charge\"]\n            test_x[event_idx, :pulsecount, 2] = event_x[\"auxiliary\"]\n            test_x[event_idx, :pulsecount, 3] = event_x[\"x\"]\n            test_x[event_idx, :pulsecount, 4] = event_x[\"y\"]\n            test_x[event_idx, :pulsecount, 5] = event_x[\"z\"]\n    \n    del batch_meta_df\n    \n    # Normalize\n    test_x[:, :, 0] /= 1000  # time\n    test_x[:, :, 1] /= 300  # charge\n    test_x[:, :, 3:] /= 600  # space\n    \n    \n    pred_mod1=models[0].predict(test_x, verbose=0)\n    gc.collect()\n    pred_mod2=models[1].predict(test_x, verbose=0)\n    gc.collect()\n    pred_mod3=models[2].predict(test_x, verbose=0)\n    gc.collect()\n    \n    pred_mod=0.3*pred_mod1+0.3*pred_mod2+0.4*pred_mod3\n\n    ","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:41:35.401889Z","iopub.execute_input":"2023-04-20T16:41:35.402630Z","iopub.status.idle":"2023-04-20T16:41:36.833573Z","shell.execute_reply.started":"2023-04-20T16:41:35.402586Z","shell.execute_reply":"2023-04-20T16:41:36.832335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nimport plotly.offline as offline\noffline.init_notebook_mode(connected=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df=pd.DataFrame(angle_bin_vector[0])\ndf['p']=pd.Series(pred_mod[0].tolist())\ndf['s']=4\n\nfig = px.scatter_3d(df, x=0, y=1, z=2,color=\"p\")\nfig.update_traces(marker_size=10)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:42:16.562778Z","iopub.execute_input":"2023-04-20T16:42:16.563925Z","iopub.status.idle":"2023-04-20T16:42:16.631397Z","shell.execute_reply.started":"2023-04-20T16:42:16.563883Z","shell.execute_reply":"2023-04-20T16:42:16.630146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df=pd.DataFrame(angle_bin_vector[0])\ndf['p']=pd.Series(pred_mod[1].tolist())\ndf['s']=4\n\nfig = px.scatter_3d(df, x=0, y=1, z=2,color=\"p\")\nfig.update_traces(marker_size=10)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:42:34.998496Z","iopub.execute_input":"2023-04-20T16:42:34.999538Z","iopub.status.idle":"2023-04-20T16:42:35.060582Z","shell.execute_reply.started":"2023-04-20T16:42:34.999487Z","shell.execute_reply":"2023-04-20T16:42:35.059433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\ndf=pd.DataFrame(angle_bin_vector[0])\ndf['p']=pd.Series(pred_mod[2].tolist())\ndf['s']=4\n\nfig = px.scatter_3d(df, x=0, y=1, z=2,color=\"p\")\nfig.update_traces(marker_size=10)\nfig.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-20T16:41:44.040064Z","iopub.execute_input":"2023-04-20T16:41:44.040443Z","iopub.status.idle":"2023-04-20T16:41:47.448676Z","shell.execute_reply.started":"2023-04-20T16:41:44.040410Z","shell.execute_reply":"2023-04-20T16:41:47.447646Z"},"trusted":true},"execution_count":null,"outputs":[]}]}