{"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 os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport gc\nimport time\nPATH_DATASET = \"/kaggle/input/icecube-neutrinos-in-deep-ice\"\nDATA_DIR=PATH_DATASET","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-28T01:15:08.601762Z","iopub.execute_input":"2023-01-28T01:15:08.602225Z","iopub.status.idle":"2023-01-28T01:15:08.609200Z","shell.execute_reply.started":"2023-01-28T01:15:08.602191Z","shell.execute_reply":"2023-01-28T01:15:08.607399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Starting with @jirkborovec's excellent notebook, we will modify the direction computation to take into account the time of the sensor detects.\n\nA separate discussion topic (TBD) will give the derivation of the formula, which is as follows.  It is called the Line-fit method.\nWe assume that the particle is following the path $pos = r + v * t$ where $v$ is the velocity and $r$ is the reference position at time 0.\n\nLet $\\langle x \\rangle$ be the average value of the quantity $x$, and \n    $ti$=vector of times of our measurements, and $ri$=vector of sensor positions (x,y,z)\n    \nThen \n    \n$$v_{est} = \\frac{ \\langle ri * ti \\rangle -  \\langle r \\rangle * \\langle t \\rangle}{\\langle t^2 \\rangle - \\langle t \\rangle ^2}$$\n    \nwhere $v_{est}$ is the estimated velocity vector of the charged particle, which we assume has the same\n    trajectory as the neutrino that generated it.\n\nReferences:  \n1. Mirco Hunnefield Masters Thesis, *Online Reconstruction of Muon-Neutrino Events in IceCube using Deep Learning Techniques*\n2. Kai Schatto, PhD Thesis, *Stacked searches for high-energy\nneutrinos from blazars with IceCube*\n","metadata":{}},{"cell_type":"code","source":"geometry = pd.read_csv(os.path.join(DATA_DIR, \"sensor_geometry.csv\"))\ngeometry.set_index(\"sensor_id\", inplace=True)\ngeometry = geometry.apply(np.float32)\ngeometry.info()","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-28T01:15:08.611588Z","iopub.execute_input":"2023-01-28T01:15:08.612934Z","iopub.status.idle":"2023-01-28T01:15:08.643814Z","shell.execute_reply.started":"2023-01-28T01:15:08.612845Z","shell.execute_reply":"2023-01-28T01:15:08.642952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#! pip install -q scikit-spatial --no-index -f /kaggle/input/icecube-neutrino-eda-3d-interactive-viewer/frozen-packages/","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-28T01:15:08.645248Z","iopub.execute_input":"2023-01-28T01:15:08.645818Z","iopub.status.idle":"2023-01-28T01:15:08.649354Z","shell.execute_reply.started":"2023-01-28T01:15:08.645787Z","shell.execute_reply":"2023-01-28T01:15:08.648594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\nimport warnings\nwarnings.filterwarnings('ignore')\n\ndef cartesian_to_sphere(x, y, z):\n    # https://en.wikipedia.org/wiki/Spherical_coordinate_system\n    x2y2 = x**2 + y**2\n    r = math.sqrt(x2y2 + z**2)\n    azimuth = math.acos(x / math.sqrt(x2y2)) * np.sign(y)\n    zenith = math.acos(z / r)\n    return azimuth, zenith\n\ndef adjust_sphere(azimuth, zenith):\n    if zenith < 0:\n        zenith += math.pi\n        azimuth += math.pi\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    azimuth = azimuth % (2 * math.pi)\n    return azimuth, zenith\n\ndef sphere_to_cartesian(azimuth, zenith):\n    # see: https://stackoverflow.com/a/10868220/4521646\n    x = math.cos(azimuth) * math.sin(zenith)\n    y = math.sin(azimuth) * math.sin(zenith)\n    z = math.cos(zenith)\n    return x, y, z","metadata":{"execution":{"iopub.status.busy":"2023-01-28T01:15:08.651385Z","iopub.execute_input":"2023-01-28T01:15:08.651686Z","iopub.status.idle":"2023-01-28T01:15:08.667847Z","shell.execute_reply.started":"2023-01-28T01:15:08.651651Z","shell.execute_reply":"2023-01-28T01:15:08.666848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Here is our function implementing the above Line-fit method**","metadata":{}},{"cell_type":"code","source":"def compute_direction(r, t):\n    \"\"\"compute Line-fit using r vector and time vectors\n    return v_est and r_est \"\"\"\n    def avg(p):\n        if len(p.shape)==1:\n            p=p.reshape(-1,1)\n        return np.mean(p, axis=0)\n    #r + v*t is path of particle\n    q = avg(t**2)-avg(t)**2\n    v_est = (avg(r*t) - avg(r)*avg(t))/q\n    r_est = avg(r) - v_est*avg(t)\n    #v_est is the direction of travel\n    return v_est, r_est\n","metadata":{"execution":{"iopub.status.busy":"2023-01-28T01:15:08.669380Z","iopub.execute_input":"2023-01-28T01:15:08.669965Z","iopub.status.idle":"2023-01-28T01:15:08.690345Z","shell.execute_reply.started":"2023-01-28T01:15:08.669925Z","shell.execute_reply":"2023-01-28T01:15:08.689211Z"},"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":"ssub = pd.read_parquet(os.path.join(PATH_DATASET, \"sample_submission.parquet\"))\nssub.set_index(\"event_id\", inplace=True)\ndisplay(ssub.head())","metadata":{"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-01-28T01:15:08.691792Z","iopub.execute_input":"2023-01-28T01:15:08.692150Z","iopub.status.idle":"2023-01-28T01:15:08.717954Z","shell.execute_reply.started":"2023-01-28T01:15:08.692122Z","shell.execute_reply":"2023-01-28T01:15:08.716577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ssub = ssub.apply(np.float32)\ndisplay(ssub.info())","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-28T01:15:08.719406Z","iopub.execute_input":"2023-01-28T01:15:08.720016Z","iopub.status.idle":"2023-01-28T01:15:08.734275Z","shell.execute_reply.started":"2023-01-28T01:15:08.719977Z","shell.execute_reply":"2023-01-28T01:15:08.733320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fitting to all data","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    df = pd.read_parquet(batch_file)\n    # del df['time'], df['charge']\n    # display(df.head())\n    gc.collect()\n    \n    for eid, dfg in df.groupby(\"event_id\"):\n        dfg = dfg[~dfg['auxiliary']]\n        dfg = dfg.merge(geometry, left_on=\"sensor_id\", right_index=True)\n        # TODO: for memory reason subsample point cloud\n        if len(dfg) > 800:\n            dfg = dfg.sort_values('charge', ascending=False)[:800]\n        # display(dfg)\n        #These are the sensor position coordinates\n        xyz = dfg[['x','y','z']].values\n        #These are the times of detection\n        ti = dfg[['time']].values\n        v_est,r_est = compute_direction(xyz, ti)\n        #we need to switch direction because we want the direction the neutrino came from, but\n        #v_est is the direction it is traveling\n        v_est = -v_est\n        azimuth_, zenith_ = cartesian_to_sphere(*v_est)\n        azimuth_, zenith_ = adjust_sphere(azimuth_, zenith_)\n        ssub.at[eid, \"azimuth\"] = azimuth_\n        ssub.at[eid, \"zenith\"] = zenith_\n        if len(df) < 1e5:\n            print(f\"Estimation event {eid} with azimuth={azimuth_} & zenith={zenith_}\")\n\n    del df, dfg\n    gc.collect()\n    time.sleep(1)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-28T01:15:08.735445Z","iopub.execute_input":"2023-01-28T01:15:08.736331Z","iopub.status.idle":"2023-01-28T01:15:10.045840Z","shell.execute_reply.started":"2023-01-28T01:15:08.736297Z","shell.execute_reply":"2023-01-28T01:15:10.044509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Finalize submission","metadata":{}},{"cell_type":"code","source":"ssub.fillna(0).to_csv('submission.csv', index=True)\n\n!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-01-28T01:15:10.048589Z","iopub.execute_input":"2023-01-28T01:15:10.049009Z","iopub.status.idle":"2023-01-28T01:15:11.153352Z","shell.execute_reply.started":"2023-01-28T01:15:10.048965Z","shell.execute_reply":"2023-01-28T01:15:11.151832Z"},"trusted":true},"execution_count":null,"outputs":[]}]}