{"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\nimport math\nPATH_DATASET = \"/kaggle/input/icecube-neutrinos-in-deep-ice\"\nDATA_DIR=PATH_DATASET\n#MODE='test'\nMODE='test'\nls = glob.glob(os.path.join(PATH_DATASET, \"test\", \"*.parquet\"))\nif len(ls)>1:\n    #make sure we run in test mode if we are running in submission environment\n    MODE='test'\nimport multiprocessing as mp\nimport queue\nmp.cpu_count()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-02-02T22:28:00.608827Z","iopub.execute_input":"2023-02-02T22:28:00.609296Z","iopub.status.idle":"2023-02-02T22:28:00.637698Z","shell.execute_reply.started":"2023-02-02T22:28:00.609260Z","shell.execute_reply":"2023-02-02T22:28:00.636723Z"},"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\nThis original method has been updated with new improvements:\n1. Select points based on a time window - we want the most primary points that fit within a window (4000nsec) but probably could be tuned.\n2. Select proximal points based on @roberthatch's excellent expository notebooks. I use a actual distance in x,y,z instead of a string_id based approach, which allows signals to be compared across sensor types (DeepCore vs regular) with a delta in x,y,z and a delta time.\n3. Use a weighted combination of primary signal, 1 or more nearby sensors, and 2 or more nearby sensors to determine the actual weights to use in a weighted version of the above least squares.  Thanks again to @roberthatch for pointing out the easy weighting concept.\n\nI did a quick search over the various meta parameters to get the best result with the time remaining before the early notebook deadline.\n","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"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-02-02T22:28:00.639444Z","iopub.execute_input":"2023-02-02T22:28:00.639999Z","iopub.status.idle":"2023-02-02T22:28:00.663839Z","shell.execute_reply.started":"2023-02-02T22:28:00.639964Z","shell.execute_reply":"2023-02-02T22:28:00.662654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#PARAMS\nDELTA_T=200\nDELTA_POS=55\nWTS=(.1,.9,0)  # weights for weight formula","metadata":{"execution":{"iopub.status.busy":"2023-02-02T22:28:00.665340Z","iopub.execute_input":"2023-02-02T22:28:00.665707Z","iopub.status.idle":"2023-02-02T22:28:00.672495Z","shell.execute_reply.started":"2023-02-02T22:28:00.665673Z","shell.execute_reply":"2023-02-02T22:28:00.671268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def nearby(df, delta_t=300, delta_dist=55):\n    \"\"\"is there a nearby event in time and space.  Return an array giving counts of nearby sensors\"\"\"\n    times = df['time'].values   # convert to different type? speed issues\n    delta_time = np.abs(times - times[:, None])\n    sensors = df['sensor_id'].values\n    delta_sensor = np.abs(sensors - sensors[:, None])\n    delta_time[(delta_sensor==0)] = 1000    # we ignore from same sensor\n    pos = df[['x','y','z']].values\n    delta_pos = np.abs(pos - pos[:, None, :])\n    #all axes less than 55 meters (17*3 = 51)?\n    close = np.all( delta_pos < delta_dist, axis = 2)\n    #close2 =  delta_pos[:,:,2] < 55 #just z axis\n    xx = (delta_time < delta_t) & (close)\n    ret = np.count_nonzero(xx, axis=1)\n    return ret\n\ndef v_to_polar(x, y, z):\n    \"\"\"Convert x,y,z to polar direction, accounting for arrival reversal\"\"\"\n    x, y, z = -x,-y,-z\n    x2y2 = x**2 + y**2\n    r = math.sqrt(x2y2 + z**2)\n    if x2y2 < 1e-6:\n        x2y2 = 1e-6\n    azimuth = math.acos(x / math.sqrt(x2y2)) * np.sign(y)\n    zenith = math.acos(z / r)  # returns 0-2pi, cannot be negative\n    azimuth, zenith = adjust_polar(azimuth, zenith)\n    return azimuth, zenith\n\ndef adjust_polar(azimuth, zenith):\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    elif zenith < 0:\n        zenith += math.pi\n    azimuth = azimuth % (2 * math.pi)\n    return azimuth, zenith","metadata":{"execution":{"iopub.status.busy":"2023-02-02T22:28:00.675302Z","iopub.execute_input":"2023-02-02T22:28:00.675665Z","iopub.status.idle":"2023-02-02T22:28:00.688550Z","shell.execute_reply.started":"2023-02-02T22:28:00.675632Z","shell.execute_reply":"2023-02-02T22:28:00.687332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def time_window(times, width=4000):\n    \"\"\"find min and max time that gives most points in fixed size window\"\"\"\n    mn = times.min()\n    mx = times.max()\n    if mx-mn < width:\n        #do not need to filter at all\n        return mn, mx\n    #try the center\n    avg = times.mean()\n    step = 25\n    best_count = 0\n    best_center = avg\n    for center in np.arange(avg - 1000, avg + 1000, step):\n        c = np.count_nonzero(np.logical_and(times > center - width/2, times < center + width/2))\n        if c > best_count:\n            best_count = c\n            best_center = center\n    if best_count < 3:\n        return mn, mx\n    m1 = max(best_center - width/2, mn)\n    m2 = min(best_center + width/2, mx)\n    return m1,m2\n","metadata":{"execution":{"iopub.status.busy":"2023-02-02T22:28:00.689734Z","iopub.execute_input":"2023-02-02T22:28:00.690630Z","iopub.status.idle":"2023-02-02T22:28:00.703816Z","shell.execute_reply.started":"2023-02-02T22:28:00.690576Z","shell.execute_reply":"2023-02-02T22:28:00.702777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Here is our function implementing the above Line-fit method, with weights added**","metadata":{}},{"cell_type":"code","source":"def compute_direction(r, t, use_weights=None):\n    \"\"\"compute r vector from time,x,y,z columns: x,y,z,time,charge\"\"\"\n    assert r.shape[1]==3\n    assert t.shape[1]==1\n    if use_weights is None:\n        weight = np.ones_like(t)\n    else:\n        weight = use_weights\n    weight = (weight/np.sum(weight)).reshape(-1,1)            \n    def avg(p):\n        if len(p.shape)==1:\n            p=p.reshape(-1,1)\n        return np.sum(p*weight, axis=0)\n    #r + v*t is path of particle\n    q = avg(t**2)-avg(t)**2\n    #import code\n    #code.interact(local=locals())\n    v_est = (avg(r*t) - avg(r)*avg(t))/q\n    r_est = avg(r) - v_est*avg(t)\n    #v_est is where the direction of travel, but we want where it came from\n    return v_est, r_est\n","metadata":{"execution":{"iopub.status.busy":"2023-02-02T22:28:00.705239Z","iopub.execute_input":"2023-02-02T22:28:00.705891Z","iopub.status.idle":"2023-02-02T22:28:00.719149Z","shell.execute_reply.started":"2023-02-02T22:28:00.705837Z","shell.execute_reply":"2023-02-02T22:28:00.718149Z"},"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-02-02T22:28:00.720882Z","iopub.execute_input":"2023-02-02T22:28:00.721466Z","iopub.status.idle":"2023-02-02T22:28:00.747505Z","shell.execute_reply.started":"2023-02-02T22:28:00.721422Z","shell.execute_reply":"2023-02-02T22:28:00.746427Z"},"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-02-02T22:28:00.795278Z","iopub.execute_input":"2023-02-02T22:28:00.795914Z","iopub.status.idle":"2023-02-02T22:28:00.810144Z","shell.execute_reply.started":"2023-02-02T22:28:00.795878Z","shell.execute_reply":"2023-02-02T22:28:00.809017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fitting to all data","metadata":{}},{"cell_type":"code","source":"def do_one_batch(batch_file, que):\n    \"\"\"Do one batch and put the results on the que\"\"\"\n    print(f\"processing: {batch_file}\")\n    result = []\n    start = time.time()\n    df_big = pd.read_parquet(batch_file)\n    # del df['time'], df['charge']\n    # display(df.head())\n    gc.collect()\n    event_count=0\n    df_big = df_big.merge(geometry, left_on=\"sensor_id\", right_index=True)\n    for eid, df in df_big.groupby(\"event_id\"):\n        t0,t1 = time_window(df[~df['auxiliary']].time, width=4000)\n        df = df[ (df.time>=t0) & (df.time <= t1) ]        \n        #take a max of 500 events, sorted properly\n        if len(df) > 500:\n            df = df.sort_values(['auxiliary','charge'], ascending=[True, False])[:500]\n        #primary sensors only\n        df_pri = df[~df['auxiliary']]  \n        #get counts of nearby sensor pulses\n        ccc = nearby(df, delta_t=DELTA_T, delta_dist=DELTA_POS)\n        if np.sum(ccc) < 5:\n            #just use primary\n            weight=None\n            use_df = df_pri\n        else:\n            weight = WTS[0]*(~df['auxiliary']).values + WTS[1]*(ccc>0) + WTS[2]*(ccc>1)\n            use_df = df\n        xyz = use_df[['x','y','z']].values\n        ti = use_df['time'].values.reshape(-1,1)\n        v_est, r_est = compute_direction(xyz, ti, use_weights=weight)\n        azimuth, zenith = v_to_polar(*v_est)\n        if len(df_big) < 1e5:\n            print(f\"Estimation event {eid} with azimuth={azimuth} & zenith={zenith}\")\n            \n        if event_count%1000 == 0 and MODE == 'train':\n            print(f'{time.time()-start:8.1f} sec {event_count}')\n        result.append((eid, azimuth, zenith))\n        if MODE=='train' and len(result)>1000:\n            break\n        event_count+=1\n    del df_big, df, df_pri\n    print(f'returning {len(result)} for {batch_file}')\n    que.put(result)\n    ","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-02-02T22:28:00.812194Z","iopub.execute_input":"2023-02-02T22:28:00.812582Z","iopub.status.idle":"2023-02-02T22:28:00.832933Z","shell.execute_reply.started":"2023-02-02T22:28:00.812548Z","shell.execute_reply":"2023-02-02T22:28:00.831440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ls = glob.glob(os.path.join(PATH_DATASET, MODE, \"*.parquet\"))\nque = mp.Queue()    \nCPUS=4\nwhile len(ls)>0 or len(mp.active_children()) > 0:\n    while len(mp.active_children()) < CPUS and len(ls)>0:\n        file=ls.pop()\n        if MODE == 'train':\n            print(f'starting batch {file}')\n        p = mp.Process(target=do_one_batch, args=(file,que))\n        p.start()\n    try:\n        result = que.get(False, 1)\n        print(f'got {len(result)} {result[:2]}')\n        for eid, azimuth, zenith in result:\n            ssub.at[eid, \"azimuth\"] = azimuth\n            ssub.at[eid, \"zenith\"] = zenith\n    except queue.Empty:\n        time.sleep(1)\ntime.sleep(5)\n#Get any remaining que items\ntry:\n    result = que.get(False, 1)\n    print(f'got {len(result)} {result[:2]}')\n    for eid, azimuth, zenith in result:\n        ssub.at[eid, \"azimuth\"] = azimuth\n        ssub.at[eid, \"zenith\"] = zenith\nexcept queue.Empty:\n    time.sleep(1)\n                ","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-02-02T22:28:00.834728Z","iopub.execute_input":"2023-02-02T22:28:00.836017Z","iopub.status.idle":"2023-02-02T22:28:06.874234Z","shell.execute_reply.started":"2023-02-02T22:28:00.835966Z","shell.execute_reply":"2023-02-02T22:28:06.872656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Finalize submission","metadata":{}},{"cell_type":"code","source":"\nssub.fillna(0).to_csv('submission.csv', index=True)\n\n!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-02-02T22:28:06.877119Z","iopub.execute_input":"2023-02-02T22:28:06.877793Z","iopub.status.idle":"2023-02-02T22:28:08.004974Z","shell.execute_reply.started":"2023-02-02T22:28:06.877754Z","shell.execute_reply":"2023-02-02T22:28:08.003439Z"},"trusted":true},"execution_count":null,"outputs":[]}]}