{"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":"# EDA for [IceCube - Neutrinos in Deep Ice](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/overview) Competition\n> **Work in progress, I welcome any feedback.**","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfiles = []\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n        \nfor file in sorted(files):\n    print(file)\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-03-27T21:40:06.681494Z","iopub.execute_input":"2023-03-27T21:40:06.682051Z","iopub.status.idle":"2023-03-27T21:40:06.827034Z","shell.execute_reply.started":"2023-03-27T21:40:06.681999Z","shell.execute_reply":"2023-03-27T21:40:06.825689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pyarrow\n!pip install fastparquet\n!pip install seaborn\n!pip install scikit-learn","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-03-27T21:40:06.829960Z","iopub.execute_input":"2023-03-27T21:40:06.830417Z","iopub.status.idle":"2023-03-27T21:40:57.374879Z","shell.execute_reply.started":"2023-03-27T21:40:06.830376Z","shell.execute_reply":"2023-03-27T21:40:57.373335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\npd.set_option('display.max_rows', 200)\nimport numpy as np\n\n%matplotlib inline\nimport matplotlib.pyplot as plt\nfrom matplotlib.lines import Line2D # for legend handle\nfrom matplotlib import cm # color map\nimport seaborn as sns\n\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.linear_model import LinearRegression","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:47:09.090042Z","iopub.execute_input":"2023-03-27T21:47:09.090439Z","iopub.status.idle":"2023-03-27T21:47:09.287622Z","shell.execute_reply.started":"2023-03-27T21:47:09.090406Z","shell.execute_reply":"2023-03-27T21:47:09.286406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ideas\n[Competition Guide](https://storage.googleapis.com/kaggle-forum-message-attachments/1958559/18618/kaggle_webinar_small.pdf)\n\n**Discussion of making a solution attempt**\n> Objective is to predict the direction of the neutrino (i.e. angular coordinates where it will exit the detector)\n- According to the *Competition Guide* (above), fitting a line to the spatial positions of detections is too simplistic.\n- They suggest to fit the Cherenkov light cone geometry (see [this brilliant visualisation](https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/discussion/381166)).\n- Additionally the Ice itself in which the detectors are housed has a scattering coefficient and absorpivity that varies with depth and wavelength. *\"Accounting for the scattering of photons in ice is crucial\"* - Host\n- As such the correct approach is to pre-simulate what the paths and detection patterns of neutrinos coming in all all possible angle and then to match test results to these.\n- Guide also mentions about using CNNs and Graph Neural networks\n\n**Next Steps**\n- Modelling phase and the attempt at a solutions seems daunting for now so focus on EDA to extract as much information as possible.\n- Treat as an EDA exercise at the very least","metadata":{}},{"cell_type":"markdown","source":"# 1. Look at each of the file types","metadata":{}},{"cell_type":"markdown","source":"### 1a) Sensor Geometry\n- just a mapping between `sensor_id` and the cartesian coordinates `x`, `y` and `z`.\n- may wish to augment this with the spherical polar coordinates: $r$, $\\theta$ and $\\phi$.","metadata":{}},{"cell_type":"code","source":"sensor_geometry = pd.read_csv(\"../input/icecube-neutrinos-in-deep-ice/sensor_geometry.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:40:58.266962Z","iopub.execute_input":"2023-03-27T21:40:58.267328Z","iopub.status.idle":"2023-03-27T21:40:58.290596Z","shell.execute_reply.started":"2023-03-27T21:40:58.267292Z","shell.execute_reply":"2023-03-27T21:40:58.289547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:40:58.291729Z","iopub.execute_input":"2023-03-27T21:40:58.292097Z","iopub.status.idle":"2023-03-27T21:40:58.326001Z","shell.execute_reply.started":"2023-03-27T21:40:58.292061Z","shell.execute_reply":"2023-03-27T21:40:58.324938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:40:58.327491Z","iopub.execute_input":"2023-03-27T21:40:58.327874Z","iopub.status.idle":"2023-03-27T21:40:58.335631Z","shell.execute_reply.started":"2023-03-27T21:40:58.327825Z","shell.execute_reply":"2023-03-27T21:40:58.334411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry.describe().T","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:40:58.337313Z","iopub.execute_input":"2023-03-27T21:40:58.337733Z","iopub.status.idle":"2023-03-27T21:40:58.381554Z","shell.execute_reply.started":"2023-03-27T21:40:58.337698Z","shell.execute_reply":"2023-03-27T21:40:58.380663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1b) Train MetaData\n- The target is the angular coordinates `azimuth` and `zenith` corresponding to $\\theta$ and $\\phi$ in the usual spherical polar coordinate system.\n- Note however that the z-axis is defined to point upwards when standing on the south pole.\n- Therefore $\\theta \\, \\epsilon \\, [0, \\pi]$ and $\\phi \\, \\epsilon \\,  [0, 2\\pi]$.","metadata":{}},{"cell_type":"code","source":"%%time\ntrain_meta = pd.read_parquet(\"../input/icecube-neutrinos-in-deep-ice/train_meta.parquet\", engine='fastparquet')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:40:58.382997Z","iopub.execute_input":"2023-03-27T21:40:58.383363Z","iopub.status.idle":"2023-03-27T21:41:40.662798Z","shell.execute_reply.started":"2023-03-27T21:40:58.383328Z","shell.execute_reply":"2023-03-27T21:41:40.661345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:40.664510Z","iopub.execute_input":"2023-03-27T21:41:40.664916Z","iopub.status.idle":"2023-03-27T21:41:40.680749Z","shell.execute_reply.started":"2023-03-27T21:41:40.664878Z","shell.execute_reply":"2023-03-27T21:41:40.679398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:40.687344Z","iopub.execute_input":"2023-03-27T21:41:40.688512Z","iopub.status.idle":"2023-03-27T21:41:40.697234Z","shell.execute_reply.started":"2023-03-27T21:41:40.688461Z","shell.execute_reply":"2023-03-27T21:41:40.696274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1c) Training Batches\n- There are 560 batches\n- Each training batch appears to contain 200,000 events.\n- Each event contains on average 164 sensor indications\n- Therefore in total ~1e8 events; comparable to the length of `train_meta` (1.3e8)","metadata":{}},{"cell_type":"code","source":"batch_number = 1","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:40.698547Z","iopub.execute_input":"2023-03-27T21:41:40.699569Z","iopub.status.idle":"2023-03-27T21:41:40.710975Z","shell.execute_reply.started":"2023-03-27T21:41:40.699529Z","shell.execute_reply":"2023-03-27T21:41:40.709798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch = pd.read_parquet(f\"../input/icecube-neutrinos-in-deep-ice/train/batch_{batch_number}.parquet\", engine='fastparquet')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:40.712544Z","iopub.execute_input":"2023-03-27T21:41:40.712909Z","iopub.status.idle":"2023-03-27T21:41:44.823147Z","shell.execute_reply.started":"2023-03-27T21:41:40.712876Z","shell.execute_reply":"2023-03-27T21:41:44.821631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:44.824818Z","iopub.execute_input":"2023-03-27T21:41:44.825270Z","iopub.status.idle":"2023-03-27T21:41:44.838925Z","shell.execute_reply.started":"2023-03-27T21:41:44.825229Z","shell.execute_reply":"2023-03-27T21:41:44.837823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:44.840681Z","iopub.execute_input":"2023-03-27T21:41:44.841422Z","iopub.status.idle":"2023-03-27T21:41:44.859901Z","shell.execute_reply.started":"2023-03-27T21:41:44.841361Z","shell.execute_reply":"2023-03-27T21:41:44.858336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch.groupby(by='event_id')['sensor_id'].count().mean()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:44.861912Z","iopub.execute_input":"2023-03-27T21:41:44.862826Z","iopub.status.idle":"2023-03-27T21:41:45.932490Z","shell.execute_reply.started":"2023-03-27T21:41:44.862775Z","shell.execute_reply":"2023-03-27T21:41:45.931223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch.index.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:45.934282Z","iopub.execute_input":"2023-03-27T21:41:45.934812Z","iopub.status.idle":"2023-03-27T21:41:46.557638Z","shell.execute_reply.started":"2023-03-27T21:41:45.934757Z","shell.execute_reply":"2023-03-27T21:41:46.556320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_batch.index.unique()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:46.559388Z","iopub.execute_input":"2023-03-27T21:41:46.560021Z","iopub.status.idle":"2023-03-27T21:41:46.788643Z","shell.execute_reply.started":"2023-03-27T21:41:46.559966Z","shell.execute_reply":"2023-03-27T21:41:46.787190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Which sensors record data the most? And with what intensity (`charge`)","metadata":{}},{"cell_type":"markdown","source":"### 2a) First find the number of hits on each detector","metadata":{}},{"cell_type":"code","source":"all_batches = pd.DataFrame()\nnum_batches = 5","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:46.790414Z","iopub.execute_input":"2023-03-27T21:41:46.790915Z","iopub.status.idle":"2023-03-27T21:41:46.797764Z","shell.execute_reply.started":"2023-03-27T21:41:46.790849Z","shell.execute_reply":"2023-03-27T21:41:46.796137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(1,num_batches+1):\n    batch = pd.read_parquet(f\"../input/icecube-neutrinos-in-deep-ice/train/batch_{i}.parquet\",\n                            engine='fastparquet')\n    all_batches = pd.concat((all_batches, batch))","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:41:46.799629Z","iopub.execute_input":"2023-03-27T21:41:46.800193Z","iopub.status.idle":"2023-03-27T21:42:22.671003Z","shell.execute_reply.started":"2023-03-27T21:41:46.800137Z","shell.execute_reply":"2023-03-27T21:42:22.667238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_batches.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:22.676989Z","iopub.execute_input":"2023-03-27T21:42:22.679470Z","iopub.status.idle":"2023-03-27T21:42:22.693529Z","shell.execute_reply.started":"2023-03-27T21:42:22.679365Z","shell.execute_reply":"2023-03-27T21:42:22.691924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_batches.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:22.696227Z","iopub.execute_input":"2023-03-27T21:42:22.697521Z","iopub.status.idle":"2023-03-27T21:42:22.726708Z","shell.execute_reply.started":"2023-03-27T21:42:22.697459Z","shell.execute_reply":"2023-03-27T21:42:22.725426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hits = all_batches['sensor_id'].value_counts().values\nsensor_ids = all_batches['sensor_id'].value_counts().index","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:22.729141Z","iopub.execute_input":"2023-03-27T21:42:22.730114Z","iopub.status.idle":"2023-03-27T21:42:24.924385Z","shell.execute_reply.started":"2023-03-27T21:42:22.730060Z","shell.execute_reply":"2023-03-27T21:42:24.922842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hits_df = pd.DataFrame({'sensor_id': sensor_ids,\n                       'hits': hits})","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:24.925744Z","iopub.execute_input":"2023-03-27T21:42:24.926127Z","iopub.status.idle":"2023-03-27T21:42:24.933291Z","shell.execute_reply.started":"2023-03-27T21:42:24.926089Z","shell.execute_reply":"2023-03-27T21:42:24.932125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hits_df","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:24.935094Z","iopub.execute_input":"2023-03-27T21:42:24.935448Z","iopub.status.idle":"2023-03-27T21:42:24.958708Z","shell.execute_reply.started":"2023-03-27T21:42:24.935414Z","shell.execute_reply":"2023-03-27T21:42:24.957414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\nsns.histplot(data=hits_df, x='hits')\nplt.xlabel('Number of hits')\nplt.title('Distribution of the number of hits on each detector')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:24.961249Z","iopub.execute_input":"2023-03-27T21:42:24.962590Z","iopub.status.idle":"2023-03-27T21:42:25.342586Z","shell.execute_reply.started":"2023-03-27T21:42:24.962526Z","shell.execute_reply":"2023-03-27T21:42:25.341042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"------\n# Aside: Issue with `sensor_geometry` coordinates\n- It appears that after grouping by the `x` coordinate, almost all the x-coordinates had 60 sensors (corresponding to the 60 Digital Optical Modules (DOMs) per string, see [page 10 of the competition guide](https://storage.googleapis.com/kaggle-forum-message-attachments/1958559/18618/kaggle_webinar_small.pdf)).\n- However one of the x-coodinates (~443m) had its 60 sensors smeared between x-values of 443.43 and 444.05.\n- Perhaps the hole that was drilled to house this string wasn't perfectly vertical, or maybe the hole is wide enough to allow the DOMs to move side to side across the above range of x-values.\n- Either way, we should correct this, so our later graphs aren't affected.\n- It turns out that almost the exact same issue occurs in the y-coordinates too.","metadata":{}},{"cell_type":"code","source":"# Only prove this is the case for x-coordinate\nsensor_group_x = sensor_geometry.groupby(by='x').agg(num_sensors=('x', 'count')).sort_index()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.344774Z","iopub.execute_input":"2023-03-27T21:42:25.345333Z","iopub.status.idle":"2023-03-27T21:42:25.374506Z","shell.execute_reply.started":"2023-03-27T21:42:25.345283Z","shell.execute_reply":"2023-03-27T21:42:25.373072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_group_x.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.376096Z","iopub.execute_input":"2023-03-27T21:42:25.376470Z","iopub.status.idle":"2023-03-27T21:42:25.387555Z","shell.execute_reply.started":"2023-03-27T21:42:25.376418Z","shell.execute_reply":"2023-03-27T21:42:25.386523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# there are some x-values that don't contain 60 sensors\nnum_sensor_freq_x = pd.DataFrame({'num_sensors': sensor_group_x['num_sensors'].value_counts().index,\n                              'frequency': sensor_group_x['num_sensors'].value_counts().values}).set_index('num_sensors')\n\nnum_sensor_freq_x","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.395354Z","iopub.execute_input":"2023-03-27T21:42:25.396506Z","iopub.status.idle":"2023-03-27T21:42:25.418044Z","shell.execute_reply.started":"2023-03-27T21:42:25.396462Z","shell.execute_reply":"2023-03-27T21:42:25.415189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# the remaining sensors, that aren't on an x-coodinate that already has 60 sensors, themselves sum to 60 i.e. 1 x-coordinate has been smeared.\nsensor_group_x.loc[sensor_group_x['num_sensors'] != 60]['num_sensors'].sum()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.420150Z","iopub.execute_input":"2023-03-27T21:42:25.420518Z","iopub.status.idle":"2023-03-27T21:42:25.433045Z","shell.execute_reply.started":"2023-03-27T21:42:25.420485Z","shell.execute_reply":"2023-03-27T21:42:25.431748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fix_sensor_geometry(coord: str):\n    sensor_group_coord = sensor_geometry.groupby(by=coord).agg(num_sensors=(coord, 'count')).sort_index()\n    faulty = sensor_group_coord.loc[sensor_group_coord['num_sensors'] != 60].index\n    faulty_mean = round(np.mean(faulty), 2)\n    \n    for i in sensor_geometry[coord].sort_values():\n        if i in faulty:\n            sensor_geometry.loc[sensor_geometry[coord] == i, coord] = faulty_mean\n        \n    return sensor_geometry","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.434990Z","iopub.execute_input":"2023-03-27T21:42:25.435341Z","iopub.status.idle":"2023-03-27T21:42:25.443987Z","shell.execute_reply.started":"2023-03-27T21:42:25.435310Z","shell.execute_reply":"2023-03-27T21:42:25.442407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry = fix_sensor_geometry('x')\nsensor_geometry = fix_sensor_geometry('y')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.446035Z","iopub.execute_input":"2023-03-27T21:42:25.446828Z","iopub.status.idle":"2023-03-27T21:42:25.565161Z","shell.execute_reply.started":"2023-03-27T21:42:25.446787Z","shell.execute_reply":"2023-03-27T21:42:25.563645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# reclculate the groupby in x and check how many detectors each x has, problem solved.\nsensor_group_x = sensor_geometry.groupby(by='x').agg(num_sensors=('x', 'count')).sort_index()\nnum_sensor_freq_x = pd.DataFrame({'num_sensors': sensor_group_x['num_sensors'].value_counts().index,\n                              'frequency': sensor_group_x['num_sensors'].value_counts().values}).set_index('num_sensors')\n\nnum_sensor_freq_x","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.567128Z","iopub.execute_input":"2023-03-27T21:42:25.567610Z","iopub.status.idle":"2023-03-27T21:42:25.589376Z","shell.execute_reply.started":"2023-03-27T21:42:25.567569Z","shell.execute_reply":"2023-03-27T21:42:25.587992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### What about the `z` coordinate?\n- Since the DOMs are placed on the strings which are lowered into the drilled holes, we don't expect there to be planes of z detectors in the x-y plane.\n- We observe that 4790 of the detectors have unique z-value with there being 185 z-values housing 2 sensors.\n- We can treat the z-values then as being continuously distributed and don't have to fix anything.","metadata":{}},{"cell_type":"code","source":"sensor_group_z = sensor_geometry.groupby(by='z').agg(num_sensors=('z', 'count')).sort_index()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.591145Z","iopub.execute_input":"2023-03-27T21:42:25.591555Z","iopub.status.idle":"2023-03-27T21:42:25.606844Z","shell.execute_reply.started":"2023-03-27T21:42:25.591518Z","shell.execute_reply":"2023-03-27T21:42:25.605094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_sensor_freq_z = pd.DataFrame({'num_sensors': sensor_group_z['num_sensors'].value_counts().index,\n                              'frequency': sensor_group_z['num_sensors'].value_counts().values}).set_index('num_sensors')\n\nnum_sensor_freq_z","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.608930Z","iopub.execute_input":"2023-03-27T21:42:25.609458Z","iopub.status.idle":"2023-03-27T21:42:25.629733Z","shell.execute_reply.started":"2023-03-27T21:42:25.609374Z","shell.execute_reply":"2023-03-27T21:42:25.628425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"----","metadata":{}},{"cell_type":"markdown","source":"### 2b) Next try to visualise the variation of hits in sensors spatially in each of the `x`, `y` & `z` directions\n- Now that the `sensor_geometry` coordinates are fixed, we can return to the task of plotting `sensor_hits` vs `position_of_sensor`.","metadata":{}},{"cell_type":"code","source":"# merge the number of hits on each sensor with the coordinates of each sensor\nsensor_geometry_aug = sensor_geometry.merge(hits_df, how='left', on='sensor_id')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.631642Z","iopub.execute_input":"2023-03-27T21:42:25.633083Z","iopub.status.idle":"2023-03-27T21:42:25.654967Z","shell.execute_reply.started":"2023-03-27T21:42:25.633031Z","shell.execute_reply":"2023-03-27T21:42:25.653476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry_aug.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.656584Z","iopub.execute_input":"2023-03-27T21:42:25.657751Z","iopub.status.idle":"2023-03-27T21:42:25.675225Z","shell.execute_reply.started":"2023-03-27T21:42:25.657698Z","shell.execute_reply":"2023-03-27T21:42:25.673526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"intensity_x = sensor_geometry_aug.groupby(by='x').agg(avg_hits=('hits', 'mean')).sort_index()\nintensity_y = sensor_geometry_aug.groupby(by='y').agg(avg_hits=('hits', 'mean')).sort_index()\nintensity_z = sensor_geometry_aug.groupby(by='z').agg(avg_hits=('hits', 'mean')).sort_index()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.677020Z","iopub.execute_input":"2023-03-27T21:42:25.677493Z","iopub.status.idle":"2023-03-27T21:42:25.714084Z","shell.execute_reply.started":"2023-03-27T21:42:25.677434Z","shell.execute_reply":"2023-03-27T21:42:25.712500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"intensity_x.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.715906Z","iopub.execute_input":"2023-03-27T21:42:25.716284Z","iopub.status.idle":"2023-03-27T21:42:25.728610Z","shell.execute_reply.started":"2023-03-27T21:42:25.716251Z","shell.execute_reply":"2023-03-27T21:42:25.727025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2b i) Spatial Distribution of Sensors\n- `x` and `y` have a normal(ish) distribution with obvious peaks at the origins.\n- This is explained by the hexagonal shape of the detector (see page 9 of the [competition guide](https://storage.googleapis.com/kaggle-forum-message-attachments/1958559/18618/kaggle_webinar_small.pdf)).\n- `z` is more or less uniform as expected.","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(1,3,figsize=(12,4))\nsns.histplot(intensity_x.index, bins=10, ax=axes[0])\nsns.histplot(intensity_y.index, bins=10, ax=axes[1])\nsns.histplot(intensity_z.index, bins=10, ax=axes[2])\nplt.suptitle('Distribution of Sensors across the Spatial Coordinates')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:25.730284Z","iopub.execute_input":"2023-03-27T21:42:25.730903Z","iopub.status.idle":"2023-03-27T21:42:26.373321Z","shell.execute_reply.started":"2023-03-27T21:42:25.730750Z","shell.execute_reply":"2023-03-27T21:42:26.371558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2b ii) Average Number of Hits vs Position\n- `x` and `y` have no variation for the most part, bar a few centrally placed sensors recording a greater number of hits (~10000)\n- the `z` coordinate is interesting: the most hits occured at low z, with slightly less at high z but then around the origin, far less hits.\n- Why might this be?","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(1,3, figsize=(15,4))\nsns.scatterplot(x=intensity_x.index, y=intensity_x['avg_hits'], ax=axes[0])\nsns.scatterplot(x=intensity_y.index, y=intensity_y['avg_hits'], ax=axes[1])\nsns.scatterplot(x=intensity_z.index, y=intensity_z['avg_hits'], ax=axes[2])","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:26.375245Z","iopub.execute_input":"2023-03-27T21:42:26.375615Z","iopub.status.idle":"2023-03-27T21:42:27.012521Z","shell.execute_reply.started":"2023-03-27T21:42:26.375580Z","shell.execute_reply":"2023-03-27T21:42:27.011320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2b iii) Average Number of Hits vs Position as a 3D plot\n- The 3D plot captures the variation in `z` well but it is harder to see the variation in `x` and `y`.\n- Looking carefully, you can see the increased intensity around the origin but it is not that clear, we need the interactive version!","metadata":{}},{"cell_type":"code","source":"sensor_aug_x = sensor_geometry_aug['x']\nsensor_aug_y = sensor_geometry_aug['y']\nsensor_aug_z = sensor_geometry_aug['z']\nsensor_aug_hits = sensor_geometry_aug['hits']","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:27.013909Z","iopub.execute_input":"2023-03-27T21:42:27.015320Z","iopub.status.idle":"2023-03-27T21:42:27.022004Z","shell.execute_reply.started":"2023-03-27T21:42:27.015277Z","shell.execute_reply":"2023-03-27T21:42:27.020189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# rescale to match the colour gradient for the plot later\nscaler = MinMaxScaler(feature_range=(0,100))\nsensor_aug_hits = scaler.fit_transform(sensor_geometry_aug['hits'].values.reshape(-1,1))","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:27.024565Z","iopub.execute_input":"2023-03-27T21:42:27.025485Z","iopub.status.idle":"2023-03-27T21:42:27.041044Z","shell.execute_reply.started":"2023-03-27T21:42:27.025426Z","shell.execute_reply":"2023-03-27T21:42:27.039994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(15,4))\ncmap = cm.plasma\n\nax = fig.add_subplot(111,projection='3d')\nax.set_xlabel('x')\nax.set_ylabel('y')\nax.set_zlabel('z')\n# detector = ax.scatter3D(sensor_x, sensor_y, sensor_z, c='darkgray', s=0.5, marker='.')\nhits = ax.scatter3D(sensor_aug_x, sensor_aug_y, sensor_aug_z, s=sensor_aug_hits, marker='o',\n                   c=sensor_aug_hits, cmap=cmap)\nax.set_title('location of hits')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T21:42:27.043139Z","iopub.execute_input":"2023-03-27T21:42:27.043841Z","iopub.status.idle":"2023-03-27T21:42:27.442254Z","shell.execute_reply.started":"2023-03-27T21:42:27.043792Z","shell.execute_reply":"2023-03-27T21:42:27.441353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Plotting an event\n- my attempt to recreate the plot shown in the competition description.\n- Can't figure out how to make interactive in matplotlib since `%matplotlib notebook` appears to be disabled.\n- need to figure out how to regress these datapoints\n- As expected due to the nature of Cherenkov radiation, the scatter plots are rather messy.","metadata":{}},{"cell_type":"code","source":"event_number = 24","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:31.326878Z","iopub.execute_input":"2023-03-27T22:01:31.327432Z","iopub.status.idle":"2023-03-27T22:01:31.334067Z","shell.execute_reply.started":"2023-03-27T22:01:31.327389Z","shell.execute_reply":"2023-03-27T22:01:31.332294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event_true = test_batch[(test_batch.index == event_number)&(test_batch['auxiliary'] == True)]\ntest_event_false = test_batch[(test_batch.index == event_number)&(test_batch['auxiliary'] == False)]","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:31.730101Z","iopub.execute_input":"2023-03-27T22:01:31.731377Z","iopub.status.idle":"2023-03-27T22:01:32.064592Z","shell.execute_reply.started":"2023-03-27T22:01:31.731325Z","shell.execute_reply":"2023-03-27T22:01:32.063188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event_true = test_event_true.merge(sensor_geometry, how='inner', left_on='sensor_id',\n                       right_on='sensor_id')\ntest_event_false = test_event_false.merge(sensor_geometry, how='inner', left_on='sensor_id',\n                       right_on='sensor_id')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:32.066459Z","iopub.execute_input":"2023-03-27T22:01:32.066908Z","iopub.status.idle":"2023-03-27T22:01:32.084511Z","shell.execute_reply.started":"2023-03-27T22:01:32.066844Z","shell.execute_reply":"2023-03-27T22:01:32.082887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_x_false = test_event_false['x']\nbatch_y_false = test_event_false['y']\nbatch_z_false = test_event_false['z']\ncharge_false = test_event_false['charge']*100\ntime_false = test_event_false['time']\nauxiliary_false = test_event_false['auxiliary']\n\nbatch_x_true = test_event_true['x']\nbatch_y_true = test_event_true['y']\nbatch_z_true = test_event_true['z']\ncharge_true = test_event_true['charge']*100\ntime_true = test_event_true['time']\nauxiliary_true = test_event_true['auxiliary']","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:32.183337Z","iopub.execute_input":"2023-03-27T22:01:32.183802Z","iopub.status.idle":"2023-03-27T22:01:32.194623Z","shell.execute_reply.started":"2023-03-27T22:01:32.183766Z","shell.execute_reply":"2023-03-27T22:01:32.192974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_x = sensor_geometry['x']\nsensor_y = sensor_geometry['y']\nsensor_z = sensor_geometry['z']","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:32.596390Z","iopub.execute_input":"2023-03-27T22:01:32.596910Z","iopub.status.idle":"2023-03-27T22:01:32.603177Z","shell.execute_reply.started":"2023-03-27T22:01:32.596854Z","shell.execute_reply":"2023-03-27T22:01:32.601436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(12,4))\ncmap = cm.coolwarm\n\nax1 = fig.add_subplot(121,projection='3d')\nax1.set_xlabel('x')\nax1.set_ylabel('y')\nax1.set_zlabel('z')\ndetector = ax1.scatter3D(sensor_x, sensor_y, sensor_z, c='darkgray', s=0.5, marker='.')\nfalse = ax1.scatter3D(batch_x_false, batch_y_false, batch_z_false, s=charge_false, marker='o', c=time_false, cmap=cmap)\nfig.colorbar(false, shrink=0.5, aspect=10)\nax1.set_title('Auxiliary == False (stronger signals)')\n\nax2 = fig.add_subplot(122,projection='3d')\nax2.set_xlabel('x')\nax2.set_ylabel('y')\nax2.set_zlabel('z')\ndetector = ax2.scatter3D(sensor_x, sensor_y, sensor_z, c='darkgray', s=0.5, marker='.')\ntrue = ax2.scatter3D(batch_x_true, batch_y_true, batch_z_true, s=charge_true, marker='o', c=time_true, cmap=cmap)\nfig.colorbar(true, shrink=0.5, aspect=10)\nax2.set_title('Auxiliary == True (weaker signals)')\n\nplt.suptitle('An Example Event Separated out by Auxiliary')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:33.036962Z","iopub.execute_input":"2023-03-27T22:01:33.037450Z","iopub.status.idle":"2023-03-27T22:01:33.733251Z","shell.execute_reply.started":"2023-03-27T22:01:33.037412Z","shell.execute_reply":"2023-03-27T22:01:33.731712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Regressing an event using a simple linear fit\n- This is not recommended as a solution attempt by the [Competition Guide](https://storage.googleapis.com/kaggle-forum-message-attachments/1958559/18618/kaggle_webinar_small.pdf) but I want to see how it would be done.\n- We will use both `auxiliary` = False and `auxiliary` = True.\n\n### 4a) The process of isolating one event","metadata":{}},{"cell_type":"code","source":"test_event = test_batch.loc[(test_batch.index == event_number)&(test_batch['auxiliary'] == True)]","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:33.801344Z","iopub.execute_input":"2023-03-27T22:01:33.801904Z","iopub.status.idle":"2023-03-27T22:01:33.955611Z","shell.execute_reply.started":"2023-03-27T22:01:33.801835Z","shell.execute_reply":"2023-03-27T22:01:33.954369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:33.968579Z","iopub.execute_input":"2023-03-27T22:01:33.969560Z","iopub.status.idle":"2023-03-27T22:01:33.984238Z","shell.execute_reply.started":"2023-03-27T22:01:33.969500Z","shell.execute_reply":"2023-03-27T22:01:33.982965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:34.104049Z","iopub.execute_input":"2023-03-27T22:01:34.104547Z","iopub.status.idle":"2023-03-27T22:01:34.113776Z","shell.execute_reply.started":"2023-03-27T22:01:34.104508Z","shell.execute_reply":"2023-03-27T22:01:34.112237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event = test_event.merge(sensor_geometry, on='sensor_id', how='inner')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:34.382541Z","iopub.execute_input":"2023-03-27T22:01:34.383460Z","iopub.status.idle":"2023-03-27T22:01:34.396440Z","shell.execute_reply.started":"2023-03-27T22:01:34.383399Z","shell.execute_reply":"2023-03-27T22:01:34.394758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_event.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:01:34.617872Z","iopub.execute_input":"2023-03-27T22:01:34.618373Z","iopub.status.idle":"2023-03-27T22:01:34.636462Z","shell.execute_reply.started":"2023-03-27T22:01:34.618331Z","shell.execute_reply":"2023-03-27T22:01:34.634915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 4b) Now try to find a best fitline for first 6 events\n- [Code Source](https://stackoverflow.com/questions/2298390/fitting-a-line-in-3d)","metadata":{}},{"cell_type":"code","source":"# take first 6 events\nevent_numbers = test_batch.index.unique()[:6]","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:19:05.233599Z","iopub.execute_input":"2023-03-27T22:19:05.234117Z","iopub.status.idle":"2023-03-27T22:19:05.470756Z","shell.execute_reply.started":"2023-03-27T22:19:05.234073Z","shell.execute_reply":"2023-03-27T22:19:05.469526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_data(event_number: int):\n    test_event = test_batch.loc[test_batch.index == event_number]\n    test_event = test_event.merge(sensor_geometry, on='sensor_id', how='inner')\n    \n    data = np.array(test_event[['x', 'y', 'z']])\n    datamean = data.mean(axis=0)\n    \n    # implement a singular value decomposition\n    uu, dd, vv = np.linalg.svd(data - datamean)\n    \n    # vv[0] contains direction vector of the 'best fit' line in the least squares sense\n    linepts = vv[0] * np.mgrid[-1e3:1e3:2j][:, np.newaxis]\n    linepts += datamean\n    \n    return data, linepts","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:19:06.698319Z","iopub.execute_input":"2023-03-27T22:19:06.698771Z","iopub.status.idle":"2023-03-27T22:19:06.707307Z","shell.execute_reply.started":"2023-03-27T22:19:06.698736Z","shell.execute_reply":"2023-03-27T22:19:06.706180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(2,3,figsize=(12,9),subplot_kw=dict(projection=\"3d\"))\n\nfor idx, event_number in enumerate(event_numbers):\n    data_f, linepts_f = fit_data(event_number)\n    axes[idx // 3, idx % 3].scatter3D(sensor_x, sensor_y, sensor_z, c='darkgray', s=0.5, marker='.')\n    axes[idx // 3, idx % 3].scatter3D(*data_f.T, marker='o')\n    axes[idx // 3, idx % 3].plot3D(*linepts_f.T)\n    axes[idx // 3, idx % 3].set_title(f'Batch {batch_number} Event {event_number}')\n    axes[idx // 3, idx % 3].set_xlabel('x')\n    axes[idx // 3, idx % 3].set_ylabel('y')\n    axes[idx // 3, idx % 3].set_zlabel('z')","metadata":{"execution":{"iopub.status.busy":"2023-03-27T22:38:35.144272Z","iopub.execute_input":"2023-03-27T22:38:35.145155Z","iopub.status.idle":"2023-03-27T22:38:37.741546Z","shell.execute_reply.started":"2023-03-27T22:38:35.145108Z","shell.execute_reply":"2023-03-27T22:38:37.740457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}