{"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":"- Thanks to [GREYSNOW](https://www.kaggle.com/code/shlomoron/icecube-eda-pca-baseline-cv-1-23-lb-1-218), I conducted EDA and PCA.\n- 理解しながら走り書きしたものです．少しでも参考になれば幸いです．","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport plotly.express as px\nfrom sklearn.decomposition import PCA\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nimport multiprocessing\nimport seaborn as sns","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# train meta","metadata":{}},{"cell_type":"code","source":"train_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train_meta.parquet')","metadata":{"execution":{"iopub.status.busy":"2023-03-30T07:55:53.851966Z","iopub.execute_input":"2023-03-30T07:55:53.852416Z","iopub.status.idle":"2023-03-30T07:56:33.626029Z","shell.execute_reply.started":"2023-03-30T07:55:53.852370Z","shell.execute_reply":"2023-03-30T07:56:33.624950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta","metadata":{"execution":{"iopub.status.busy":"2023-03-30T07:56:33.635447Z","iopub.execute_input":"2023-03-30T07:56:33.635972Z","iopub.status.idle":"2023-03-30T07:56:33.680781Z","shell.execute_reply.started":"2023-03-30T07:56:33.635938Z","shell.execute_reply":"2023-03-30T07:56:33.679828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta.describe()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T07:59:29.367146Z","iopub.execute_input":"2023-03-30T07:59:29.367835Z","iopub.status.idle":"2023-03-30T08:01:17.553554Z","shell.execute_reply.started":"2023-03-30T07:59:29.367794Z","shell.execute_reply":"2023-03-30T08:01:17.552555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta[train_meta['batch_id'] == 1]","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:02:40.702872Z","iopub.execute_input":"2023-03-30T08:02:40.703272Z","iopub.status.idle":"2023-03-30T08:02:41.166742Z","shell.execute_reply.started":"2023-03-30T08:02:40.703237Z","shell.execute_reply":"2023-03-30T08:02:41.165481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"what is first_pulse_index and last_pulse_index?\n- event_idのtrain_batch_{n}のindexに該当する","metadata":{}},{"cell_type":"markdown","source":"# train batch 1","metadata":{}},{"cell_type":"code","source":"train_batch_1 = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/train/batch_1.parquet')\ntrain_batch_1","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:03:25.357665Z","iopub.execute_input":"2023-03-30T08:03:25.358620Z","iopub.status.idle":"2023-03-30T08:03:27.222291Z","shell.execute_reply.started":"2023-03-30T08:03:25.358575Z","shell.execute_reply":"2023-03-30T08:03:27.221119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 各batchごとのevent_idが格納されている\n- auxiliaryは不確か度合いのようなもので，Falseの時は十分確かなデータが取られているが，Trueの時は微妙みたいなニュアンス","metadata":{}},{"cell_type":"markdown","source":"# sensor geometry","metadata":{}},{"cell_type":"code","source":"sensor_geometry = pd.read_csv('/kaggle/input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv')\nsensor_geometry","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:04:31.282336Z","iopub.execute_input":"2023-03-30T08:04:31.282740Z","iopub.status.idle":"2023-03-30T08:04:31.303902Z","shell.execute_reply.started":"2023-03-30T08:04:31.282704Z","shell.execute_reply":"2023-03-30T08:04:31.302761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- センサーは5160個存在する\n- 各センサーの位置が記されている","metadata":{}},{"cell_type":"code","source":"plt.scatter(sensor_geometry.x, sensor_geometry.y)\nplt.xlabel('x')\nplt.ylabel('y')","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:06:11.220809Z","iopub.execute_input":"2023-03-30T08:06:11.221249Z","iopub.status.idle":"2023-03-30T08:06:11.478010Z","shell.execute_reply.started":"2023-03-30T08:06:11.221210Z","shell.execute_reply":"2023-03-30T08:06:11.476684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(sensor_geometry.x, sensor_geometry.z)\nplt.xlabel('x')\nplt.ylabel('z')","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:06:35.848390Z","iopub.execute_input":"2023-03-30T08:06:35.848800Z","iopub.status.idle":"2023-03-30T08:06:36.103842Z","shell.execute_reply.started":"2023-03-30T08:06:35.848766Z","shell.execute_reply":"2023-03-30T08:06:36.102647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 地面に垂直な棒がたくさん並んでいるイメージ","metadata":{}},{"cell_type":"code","source":"fig = px.scatter_3d(sensor_geometry, x='x', y='y', z='z', opacity=0.5)\nfig.update_traces(marker_size=2)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:05:20.829886Z","iopub.execute_input":"2023-03-30T08:05:20.830269Z","iopub.status.idle":"2023-03-30T08:05:22.425268Z","shell.execute_reply.started":"2023-03-30T08:05:20.830234Z","shell.execute_reply":"2023-03-30T08:05:22.424385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test meta","metadata":{}},{"cell_type":"code","source":"test_meta = pd.read_parquet('/kaggle/input/icecube-neutrinos-in-deep-ice/test_meta.parquet')\ntest_meta","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:07:04.942262Z","iopub.execute_input":"2023-03-30T08:07:04.943324Z","iopub.status.idle":"2023-03-30T08:07:04.960081Z","shell.execute_reply.started":"2023-03-30T08:07:04.943263Z","shell.execute_reply":"2023-03-30T08:07:04.959131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 3データしかない（Trainは131953924データ存在）","metadata":{}},{"cell_type":"markdown","source":"# understanding the data","metadata":{}},{"cell_type":"markdown","source":"試しにcase_event_idx=3(event_id=67)を取り出して見てみている","metadata":{}},{"cell_type":"code","source":"case_event_idx = 3\ncase_event_pulses = train_batch_1.iloc[\n    train_meta.iloc[case_event_idx].first_pulse_index.astype(int): \n    train_meta.iloc[case_event_idx].last_pulse_index.astype(int)+1\n].copy()\nprint(len(case_event_pulses))\ncase_event_pulses","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:11:34.970283Z","iopub.execute_input":"2023-03-30T08:11:34.971468Z","iopub.status.idle":"2023-03-30T08:11:34.990904Z","shell.execute_reply.started":"2023-03-30T08:11:34.971424Z","shell.execute_reply":"2023-03-30T08:11:34.989672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- event_id=67を抽出している\n- （first_pulse_indexとかで指定する必要はなくて，event_idで指定するだけで良く無い？）","metadata":{}},{"cell_type":"code","source":"case_event_pulses = case_event_pulses.reset_index()\ncase_event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[case_event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\ncase_event_pulses","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:16:09.274117Z","iopub.execute_input":"2023-03-30T08:16:09.274660Z","iopub.status.idle":"2023-03-30T08:16:09.305824Z","shell.execute_reply.started":"2023-03-30T08:16:09.274604Z","shell.execute_reply":"2023-03-30T08:16:09.304371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"case_event_pulses.describe()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:16:54.885292Z","iopub.execute_input":"2023-03-30T08:16:54.886013Z","iopub.status.idle":"2023-03-30T08:16:54.986008Z","shell.execute_reply.started":"2023-03-30T08:16:54.885946Z","shell.execute_reply":"2023-03-30T08:16:54.984117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.catplot(case_event_pulses[['x', 'y', 'z']], kind='bar')","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:17:36.944117Z","iopub.execute_input":"2023-03-30T08:17:36.945037Z","iopub.status.idle":"2023-03-30T08:17:37.210707Z","shell.execute_reply.started":"2023-03-30T08:17:36.944993Z","shell.execute_reply":"2023-03-30T08:17:37.209552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- かなりマイナス側に寄っている","metadata":{}},{"cell_type":"code","source":"case_event_pulses.x = case_event_pulses.x - case_event_pulses.x.mean()\ncase_event_pulses.y = case_event_pulses.y - case_event_pulses.y.mean()\ncase_event_pulses.z = case_event_pulses.z - case_event_pulses.z.mean()\ncase_event_pulses.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:18:18.263396Z","iopub.execute_input":"2023-03-30T08:18:18.263820Z","iopub.status.idle":"2023-03-30T08:18:18.283457Z","shell.execute_reply.started":"2023-03-30T08:18:18.263781Z","shell.execute_reply":"2023-03-30T08:18:18.282056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.catplot(case_event_pulses[['x', 'y', 'z']], kind='bar')","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:18:21.574048Z","iopub.execute_input":"2023-03-30T08:18:21.575325Z","iopub.status.idle":"2023-03-30T08:18:21.853406Z","shell.execute_reply.started":"2023-03-30T08:18:21.575257Z","shell.execute_reply":"2023-03-30T08:18:21.852400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 平均値で補正","metadata":{}},{"cell_type":"markdown","source":"# calclate the vector","metadata":{}},{"cell_type":"code","source":"zenith_target = train_meta.iloc[case_event_idx].zenith\nazimuth_target = train_meta.iloc[case_event_idx].azimuth\n","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:19:23.086461Z","iopub.execute_input":"2023-03-30T08:19:23.086903Z","iopub.status.idle":"2023-03-30T08:19:23.093776Z","shell.execute_reply.started":"2023-03-30T08:19:23.086866Z","shell.execute_reply":"2023-03-30T08:19:23.092461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zenith_target, azimuth_target","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:34:34.780021Z","iopub.execute_input":"2023-03-30T08:34:34.780530Z","iopub.status.idle":"2023-03-30T08:34:34.788827Z","shell.execute_reply.started":"2023-03-30T08:34:34.780484Z","shell.execute_reply":"2023-03-30T08:34:34.787495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 方位角．\n    - zenithが首の上下の傾き\n    - azimuthが首の左右の回転\n- というイメージ","metadata":{}},{"cell_type":"code","source":"np.cos(zenith_target)","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:34:52.709396Z","iopub.execute_input":"2023-03-30T08:34:52.709872Z","iopub.status.idle":"2023-03-30T08:34:52.719404Z","shell.execute_reply.started":"2023-03-30T08:34:52.709798Z","shell.execute_reply":"2023-03-30T08:34:52.717714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vector_target = [\n    np.cos(azimuth_target) * np.sin(zenith_target), \n    np.sin(azimuth_target) * np.sin(zenith_target),\n    np.cos(zenith_target)]","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:20:38.740674Z","iopub.execute_input":"2023-03-30T08:20:38.741130Z","iopub.status.idle":"2023-03-30T08:20:38.746830Z","shell.execute_reply.started":"2023-03-30T08:20:38.741087Z","shell.execute_reply":"2023-03-30T08:20:38.745540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 方位角を3次元ベクトルに変換","metadata":{}},{"cell_type":"code","source":"vector_target","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:20:42.634938Z","iopub.execute_input":"2023-03-30T08:20:42.635420Z","iopub.status.idle":"2023-03-30T08:20:42.643073Z","shell.execute_reply.started":"2023-03-30T08:20:42.635375Z","shell.execute_reply":"2023-03-30T08:20:42.641827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vector_base = np.array([-500, 500])\nx = vector_base*vector_target[0]\ny = vector_base*vector_target[1]\nz = vector_base*vector_target[2]\nvector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\nvector_target_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:21:04.891183Z","iopub.execute_input":"2023-03-30T08:21:04.892180Z","iopub.status.idle":"2023-03-30T08:21:04.904726Z","shell.execute_reply.started":"2023-03-30T08:21:04.892130Z","shell.execute_reply":"2023-03-30T08:21:04.903900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Ground Truthの三次元ベクトルがもとまった","metadata":{}},{"cell_type":"code","source":"fig_auxiliary = px.scatter_3d(\n    case_event_pulses.loc[case_event_pulses.auxiliary],\n    x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\nfig_non_auxiliary = px.scatter_3d(\n    case_event_pulses.loc[~case_event_pulses.auxiliary],\n    x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\nfig_line = px.line_3d(\n    vector_target_df, x=\"x\", y=\"y\", z=\"z\")\n\nfig_auxiliary.update_traces(marker_size=2)\nfig_non_auxiliary.update_traces(marker_size=2)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:39:05.771532Z","iopub.execute_input":"2023-03-30T08:39:05.772000Z","iopub.status.idle":"2023-03-30T08:39:05.960873Z","shell.execute_reply.started":"2023-03-30T08:39:05.771954Z","shell.execute_reply":"2023-03-30T08:39:05.959556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- red is auxiliary is True, blue is auxiliary is False","metadata":{}},{"cell_type":"markdown","source":"# PCA","metadata":{}},{"cell_type":"code","source":"pca = PCA(n_components=1).fit(case_event_pulses[['x', 'y', 'z']])","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:49:55.630936Z","iopub.execute_input":"2023-03-30T08:49:55.631418Z","iopub.status.idle":"2023-03-30T08:49:55.641896Z","shell.execute_reply.started":"2023-03-30T08:49:55.631377Z","shell.execute_reply":"2023-03-30T08:49:55.640393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pca.components_","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:49:56.771578Z","iopub.execute_input":"2023-03-30T08:49:56.772020Z","iopub.status.idle":"2023-03-30T08:49:56.780948Z","shell.execute_reply.started":"2023-03-30T08:49:56.771973Z","shell.execute_reply":"2023-03-30T08:49:56.779374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- xyzの第一主成分を抽出している","metadata":{}},{"cell_type":"code","source":"vector = pca.components_[0]\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})\n\nfig_auxiliary = px.scatter_3d(case_event_pulses.loc[case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\nfig_non_auxiliary = px.scatter_3d(case_event_pulses.loc[~case_event_pulses.auxiliary],\n                              x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\nfig_line = px.line_3d(vector_target_df, x=\"x\", y=\"y\", z=\"z\")\nfig_fit_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\", color_discrete_sequence=['magenta'])\n\nfig_auxiliary.update_traces(marker_size=2)\nfig_non_auxiliary.update_traces(marker_size=2)\n\nfig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data + fig_fit_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T08:50:18.907784Z","iopub.execute_input":"2023-03-30T08:50:18.908205Z","iopub.status.idle":"2023-03-30T08:50:19.145508Z","shell.execute_reply.started":"2023-03-30T08:50:18.908167Z","shell.execute_reply":"2023-03-30T08:50:19.144165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- The figure demonstrates the alignment of the angles with the 'true' pulses. Unfortunately, many events (probably most) do not have a clear pulse line like this one (hence I chose it for demonstration). Try to change case_event_idx to other numbers, then rerun the last four cells and see for yourself :) An important thing to note is that most non-auxiliary pulses appear as a short vertical line of close dots, not because they come straight from the top but due to the structure of the detector. Watch [this](https://www.youtube.com/watch?v=iv-Rz3-s4BM) video for a better explanation, and in particular, 3:12-3:25.\n\n- 本当か確かめるために１０パターン試してみた","metadata":{}},{"cell_type":"code","source":"for case_event_idx in range(10):\n    case_event_pulses = train_batch_1.iloc[\n        train_meta.iloc[case_event_idx].first_pulse_index.astype(int): \n        train_meta.iloc[case_event_idx].last_pulse_index.astype(int)+1\n    ].copy()\n    print(len(case_event_pulses))\n    case_event_pulses\n\n    case_event_pulses = case_event_pulses.reset_index()\n    case_event_pulses[['x', 'y', 'z']] = sensor_geometry.loc[case_event_pulses.sensor_id].reset_index()[['x', 'y', 'z']]\n    case_event_pulses.x = case_event_pulses.x - case_event_pulses.x.mean()\n    case_event_pulses.y = case_event_pulses.y - case_event_pulses.y.mean()\n    case_event_pulses.z = case_event_pulses.z - case_event_pulses.z.mean()\n    case_event_pulses.head()\n\n    zenith_target = train_meta.iloc[case_event_idx].zenith\n    azimuth_target = train_meta.iloc[case_event_idx].azimuth\n\n    vector_target = [\n        np.cos(azimuth_target) * np.sin(zenith_target), \n        np.sin(azimuth_target) * np.sin(zenith_target),\n        np.cos(zenith_target)]\n\n    vector_base = np.array([-500, 500])\n    x = vector_base*vector_target[0]\n    y = vector_base*vector_target[1]\n    z = vector_base*vector_target[2]\n    vector_target_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n    vector_target_df.head()\n\n    pca = PCA(n_components=1).fit(case_event_pulses[['x', 'y', 'z']])\n    vector = pca.components_[0]\n    vector_base = np.array([-500, 500])\n    x = vector_base*vector[0]\n    y = vector_base*vector[1]\n    z = vector_base*vector[2]\n    vector_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\n\n    fig_auxiliary = px.scatter_3d(case_event_pulses.loc[case_event_pulses.auxiliary],\n                                  x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['red'])\n    fig_non_auxiliary = px.scatter_3d(case_event_pulses.loc[~case_event_pulses.auxiliary],\n                                  x='x', y='y', z='z', opacity=0.5, color_discrete_sequence=['blue'])\n    fig_line = px.line_3d(vector_target_df, x=\"x\", y=\"y\", z=\"z\")\n    fig_fit_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\", color_discrete_sequence=['magenta'])\n\n    fig_auxiliary.update_traces(marker_size=2)\n    fig_non_auxiliary.update_traces(marker_size=2)\n\n    fig = go.Figure(data = fig_auxiliary.data + fig_non_auxiliary.data + fig_line.data + fig_fit_line.data)\n    print('red: predict, blue: true')\n    fig.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-30T09:10:31.297961Z","iopub.execute_input":"2023-03-30T09:10:31.298455Z","iopub.status.idle":"2023-03-30T09:10:34.401921Z","shell.execute_reply.started":"2023-03-30T09:10:31.298410Z","shell.execute_reply":"2023-03-30T09:10:34.400587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- 確かに予測（赤）は正解（青）と比べて全然違う\n    - 一見赤はauxiliary=Falseの様子を十分捉えているような感じに見えるが，正解角度は全く別を示している（例：event_id=43）\n    - 時系列情報を含まない予測なので，時系列を含めたら上手い感じで予測できるのかも？LSTMに期待という気持ちです","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}