{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# IceCube❄️: baseline with XGBoost regression\n\nThe goal of this competition is to predict a neutrino particle’s direction. You will develop a model based on data from the \"IceCube\" detector, which observes the cosmos from deep within the South Pole ice.\n\nThe IceCube Neutrino Observatory is the first detector of its kind, encompassing a cubic kilometer of ice and designed to search for the nearly massless neutrinos. An international group of scientists is responsible for the scientific research that makes up the IceCube Collaboration.\n\nBy making the process faster and more precise, you'll help improve the reconstruction of neutrinos. As a result, we could gain a clearer image of our universe.\n\nSee EDA: [**IceCube🧊: Neutrino🎆EDA & 3D🔭interactive viewer**](https://www.kaggle.com/code/jirkaborovec/icecube-neutrino-eda-3d-interactive-viewer)\n\n**Tutorials:**\n- https://stackabuse.com/bytes/end-to-end-xgboost-regression-pipeline-with-scikit-learn/\n- https://xgboost.readthedocs.io/en/stable/python/examples/multioutput_regression.html","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nPATH_DATASET = \"/kaggle/input/icecube-neutrinos-in-deep-ice\"","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-23T10:28:11.597850Z","iopub.execute_input":"2023-01-23T10:28:11.598148Z","iopub.status.idle":"2023-01-23T10:28:11.609917Z","shell.execute_reply.started":"2023-01-23T10:28:11.598120Z","shell.execute_reply":"2023-01-23T10:28:11.608850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse meta data\n\n**[train/test]_meta.parquet**\n\n- **batch_id** (int): the ID of the batch the event was placed into.\n- **event_id** (int): the event ID.\n- **[first/last]_pulse_index** (int): index of the first/last row in the features dataframe belonging to this event.\n- **[azimuth/zenith]** (float32): the [azimuth/zenith] angle in radians of the neutrino. A value between 0 and 2*pi for the azimuth and 0 and pi for zenith. The target columns. Not provided for the test set. The direction vector represented by zenith and azimuth points to where the neutrino came from.","metadata":{}},{"cell_type":"code","source":"meta_test = pd.read_parquet(os.path.join(PATH_DATASET, \"test_meta.parquet\"))\nprint(f\"length: {len(meta_test)}\")\nmeta_test.head()\ndel meta_test","metadata":{"execution":{"iopub.status.busy":"2023-01-23T10:28:11.611219Z","iopub.execute_input":"2023-01-23T10:28:11.611676Z","iopub.status.idle":"2023-01-23T10:28:11.749602Z","shell.execute_reply.started":"2023-01-23T10:28:11.611631Z","shell.execute_reply":"2023-01-23T10:28:11.748410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_train = pd.read_parquet(os.path.join(PATH_DATASET, \"train_meta.parquet\"))\nprint(f\"total events: {len(meta_train)}\")\nprint(f\"nb batches: {len(meta_train['batch_id'].unique())}\")\nmeta_train.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T10:28:11.751884Z","iopub.execute_input":"2023-01-23T10:28:11.752427Z","iopub.status.idle":"2023-01-23T10:28:54.618625Z","shell.execute_reply.started":"2023-01-23T10:28:11.752384Z","shell.execute_reply":"2023-01-23T10:28:54.617539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### azimuth and zenith","metadata":{}},{"cell_type":"code","source":"plt.hist2d(meta_train[\"azimuth\"], meta_train[\"zenith\"], bins=(50, 50), cmap=plt.cm.jet)\nplt.xlabel('azimuth'), plt.ylabel('zenith')\nplt.colorbar()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T10:28:54.619773Z","iopub.execute_input":"2023-01-23T10:28:54.620649Z","iopub.status.idle":"2023-01-23T10:29:11.514808Z","shell.execute_reply.started":"2023-01-23T10:28:54.620604Z","shell.execute_reply":"2023-01-23T10:29:11.513749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"meta_train_ = meta_train[meta_train['batch_id'] == 10]\nprint(f\"batch events: {len(meta_train_)}\")\nplt.hist2d(meta_train_[\"azimuth\"], meta_train_[\"zenith\"], bins=(50, 50), cmap=plt.cm.jet)\nplt.xlabel('azimuth'), plt.ylabel('zenith')\nplt.colorbar()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T10:29:11.516116Z","iopub.execute_input":"2023-01-23T10:29:11.516501Z","iopub.status.idle":"2023-01-23T10:29:12.269544Z","shell.execute_reply.started":"2023-01-23T10:29:11.516467Z","shell.execute_reply":"2023-01-23T10:29:12.268536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Browse training/test data\n\n**[train/test]/batch_[n].parquet** Each batch contains tens of thousands of events. Each event may contain thousands of pulses, each of which is the digitized output from a photomultiplier tube and occupies one row.\n\n- **event_id** (int): the event ID. Saved as the index column in parquet.\n- **time** (int): the time of the pulse in nanoseconds in the current event time window. The absolute time of a pulse has no relevance, and only the relative time with respect to other pulses within an event is of relevance.\n- **sensor_id** (int): the ID of which of the 5160 IceCube photomultiplier sensors recorded this pulse.\n- **charge** (float32): An estimate of the amount of light in the pulse, in units of photoelectrons (p.e.). A physical photon does not exactly result in a measurement of 1 p.e. but rather can take values spread around 1 p.e. As an example, a pulse with charge 2.7 p.e. could quite likely be the result of two or three photons hitting the photomultiplier tube around the same time. This data has float16 precision but is stored as float32 due to limitations of the version of pyarrow the data was prepared with.\n- **auxiliary** (bool): If True, the pulse was not fully digitized, is of lower quality, and was more likely to originate from noise. If False, then this pulse was contributed to the trigger decision and the pulse was fully digitized.","metadata":{}},{"cell_type":"code","source":"df_test = pd.read_parquet(os.path.join(PATH_DATASET, \"test/batch_661.parquet\"))\nprint(f\"length: {len(df_test)}\")\nprint(f\"events: {len(df_test.index.unique())}\")\ndisplay(df_test.head())\ndel df_test","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:06:24.719188Z","iopub.execute_input":"2023-01-23T06:06:24.720000Z","iopub.status.idle":"2023-01-23T06:06:24.748644Z","shell.execute_reply.started":"2023-01-23T06:06:24.719894Z","shell.execute_reply":"2023-01-23T06:06:24.747320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SKIP: as this is already part of EDA\n\n# df_train = pd.read_parquet(os.path.join(PATH_DATASET, \"train/batch_1.parquet\"))\n# print(f\"length: {len(df_train)}\")\n# print(f\"events: {len(df_train.index.unique())}\")\n# display(df_train.head())\n# del df_train","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:06:24.750843Z","iopub.execute_input":"2023-01-23T06:06:24.751315Z","iopub.status.idle":"2023-01-23T06:06:24.755874Z","shell.execute_reply.started":"2023-01-23T06:06:24.751268Z","shell.execute_reply":"2023-01-23T06:06:24.755018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Converting training data\n\nconverting data where single row is one event so sensors are columns\n\n**NOTE:** avoid usung pandas as it will blow memory","metadata":{}},{"cell_type":"code","source":"from tqdm.auto import tqdm\n\ndef transform_batch(path_batch_parquet, nb_sensors=5160):\n    df = pd.read_parquet(path_batch_parquet)\n    # df[df['auxiliary']]['charge'] /= 10.\n    data = np.zeros((len(df.index.unique()), nb_sensors), dtype=np.float16)\n    event_ids = []\n    for i, (idx, dfg) in tqdm(enumerate(df.groupby(level=0))):\n        event_ids.append(idx)\n        data[i, dfg['sensor_id']] = dfg['charge']\n    # df = pd.DataFrame(data, index=event_ids, columns=[f\"sid-{i}\" for i in range(nb_sensors)])\n    return event_ids, data","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:06:24.756890Z","iopub.execute_input":"2023-01-23T06:06:24.757197Z","iopub.status.idle":"2023-01-23T06:06:24.896486Z","shell.execute_reply.started":"2023-01-23T06:06:24.757168Z","shell.execute_reply":"2023-01-23T06:06:24.895164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"use just single batch, maybe depending on memory consimptions we can use more...","metadata":{}},{"cell_type":"code","source":"event_ids, data = transform_batch(os.path.join(PATH_DATASET, \"train/batch_10.parquet\"))\nprint(f\"data size: {data.shape}\")\n\nmeta_train_ = meta_train[meta_train['event_id'].isin(event_ids)]\nmeta_train_ = dict(zip(meta_train_['event_id'].values, meta_train_[[\"azimuth\", \"zenith\"]].values.tolist()))\nprint(f\"LUT size: {len(meta_train_)}\")\nangles = np.array([meta_train_[eid] for eid in event_ids], dtype=np.float16)\nprint(f\"angles size: {angles.shape}\")","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:06:24.902139Z","iopub.execute_input":"2023-01-23T06:06:24.902527Z","iopub.status.idle":"2023-01-23T06:07:07.230666Z","shell.execute_reply.started":"2023-01-23T06:06:24.902496Z","shell.execute_reply":"2023-01-23T06:07:07.229585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data preprocessing","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_test, y_train, y_test = train_test_split(data, angles, train_size=0.8)\ndel data, angles","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:07:07.231986Z","iopub.execute_input":"2023-01-23T06:07:07.232383Z","iopub.status.idle":"2023-01-23T06:07:08.569102Z","shell.execute_reply.started":"2023-01-23T06:07:07.232325Z","shell.execute_reply":"2023-01-23T06:07:08.567885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom xgboost import XGBRegressor\n\npreprocess = Pipeline([\n    ('scaler', StandardScaler()),\n    # https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html\n    (\"PCA\", PCA(\n        n_components=510,\n        copy=False,\n    )),\n])\n\nX_train = preprocess.fit_transform(X_train)\nX_test = preprocess.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T06:07:08.570690Z","iopub.execute_input":"2023-01-23T06:07:08.571058Z","iopub.status.idle":"2023-01-23T06:13:00.527653Z","shell.execute_reply.started":"2023-01-23T06:07:08.571023Z","shell.execute_reply":"2023-01-23T06:13:00.525122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Regression baseline\n\nExplore: https://www.kaggle.com/code/prashant111/a-guide-on-xgboost-hyperparameters-tuning/notebook\n\nLearning rate tuning: https://towardsdatascience.com/betaboosting-2cd6c697eb93","metadata":{}},{"cell_type":"markdown","source":"Fit to the data... yes, validation on tests set is not the best but simplest one for now :)","metadata":{}},{"cell_type":"code","source":"# https://xgboost.readthedocs.io/en/stable/parameter.html\nregressor = XGBRegressor(\n    # tree_method=\"hist\",\n    n_estimators=256,\n    learning_rate=0.02,\n    # grow_policy=\"lossguide\",\n    num_target=y_train.shape[1],\n    verbosity=2,  # 0 (silent), 1 (warning), 2 (info), 3 (debug)\n)\n\nregressor.fit(\n    X_train, y_train,\n    # extra arguments: https://stackoverflow.com/a/35634198/4521646\n    eval_set=[(X_test, y_test)], \n    verbose=True,\n    early_stopping_rounds=5,\n)\n\ndel X_train, y_train","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-23T06:13:00.536152Z","iopub.execute_input":"2023-01-23T06:13:00.536962Z","iopub.status.idle":"2023-01-23T07:27:36.130795Z","shell.execute_reply.started":"2023-01-23T06:13:00.536890Z","shell.execute_reply":"2023-01-23T07:27:36.129432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results = regressor.evals_result()\ndf_res = pd.DataFrame(results['validation_0'])\ndf_res.plot()\nplt.xlabel('progress'), plt.grid()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:36.133152Z","iopub.execute_input":"2023-01-23T07:27:36.133525Z","iopub.status.idle":"2023-01-23T07:27:36.439715Z","shell.execute_reply.started":"2023-01-23T07:27:36.133492Z","shell.execute_reply":"2023-01-23T07:27:36.438273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r2 = regressor.score(X_test, y_test)\nprint(f\"XGBoost regressor r2: {r2}\")\ndel X_test, y_test","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:36.441350Z","iopub.execute_input":"2023-01-23T07:27:36.441834Z","iopub.status.idle":"2023-01-23T07:27:36.797428Z","shell.execute_reply.started":"2023-01-23T07:27:36.441800Z","shell.execute_reply":"2023-01-23T07:27:36.796214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare submission\n\nAn example submission with the correct columns and properly ordered event IDs. The sample submission is provided in the parquet format so it can be read quickly but your final submission must be a csv.","metadata":{}},{"cell_type":"code","source":"import gc, time\ngc.collect()\ntime.sleep(9)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T10:29:12.270899Z","iopub.execute_input":"2023-01-23T10:29:12.271201Z","iopub.status.idle":"2023-01-23T10:29:21.385219Z","shell.execute_reply.started":"2023-01-23T10:29:12.271173Z","shell.execute_reply":"2023-01-23T10:29:21.384409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ssub = pd.read_parquet(os.path.join(PATH_DATASET, \"sample_submission.parquet\"))\n# dtypes={\"azimuth\": np.float16, \"zenith\": np.float16}\nprint(f\"length: {len(ssub)}\")\nssub.set_index(\"event_id\", inplace=True)\ndisplay(ssub.head())\nssub = ssub.apply(np.float16)\ndisplay(ssub.info())","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:36.799571Z","iopub.execute_input":"2023-01-23T07:27:36.800471Z","iopub.status.idle":"2023-01-23T07:27:36.857960Z","shell.execute_reply.started":"2023-01-23T07:27:36.800417Z","shell.execute_reply":"2023-01-23T07:27:36.856878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fill defaults","metadata":{}},{"cell_type":"code","source":"# azimuth seems to have random distrib so haveing 1. does not matter\n# ssub['azimuth'] = meta_train['azimuth'].median()\n\n# zenith seesms to have Gausina distribution\nssub['zenith'] = meta_train['zenith'].mean()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:36.859406Z","iopub.execute_input":"2023-01-23T07:27:36.859754Z","iopub.status.idle":"2023-01-23T07:27:38.056734Z","shell.execute_reply.started":"2023-01-23T07:27:36.859723Z","shell.execute_reply":"2023-01-23T07:27:38.055452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Predictions / inference","metadata":{}},{"cell_type":"code","source":"ls = glob.glob(os.path.join(PATH_DATASET, \"test\", \"*.parquet\"))\n\nfor batch_file in ls:\n    print(f\"processing: {batch_file}\")\n    event_ids, data = transform_batch(batch_file)\n    preds = regressor.predict(preprocess.transform(data))\n    del data\n    #print(preds)\n    for eid, (a, z) in zip(event_ids, preds):\n        ssub.at[eid, \"azimuth\"] = a\n        ssub.at[eid, \"zenith\"] = z\n    del event_ids, preds\n    gc.collect()\n    time.sleep(9)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:38.058465Z","iopub.execute_input":"2023-01-23T07:27:38.058945Z","iopub.status.idle":"2023-01-23T07:27:38.197853Z","shell.execute_reply.started":"2023-01-23T07:27:38.058898Z","shell.execute_reply":"2023-01-23T07:27:38.196432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ssub.to_csv('submission.csv', index=True)\n\n!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-01-23T07:27:38.200542Z","iopub.execute_input":"2023-01-23T07:27:38.201432Z","iopub.status.idle":"2023-01-23T07:27:39.575305Z","shell.execute_reply.started":"2023-01-23T07:27:38.201390Z","shell.execute_reply":"2023-01-23T07:27:39.573925Z"},"trusted":true},"execution_count":null,"outputs":[]}]}