{"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":"**The goal of this competition is to predict a neutrino particle’s direction. I will develop a model based on data from the \"IceCube\" detector, which observes the cosmos from deep within the South Pole ice.**\n\n**My work could help scientists better understand exploding stars, gamma-ray bursts, and cataclysmic phenomena involving black holes, neutron stars and the fundamental properties of the neutrino itself.**\n\n**One of the most abundant particles in the universe is the neutrino. While similar to an electron, the nearly massless and electrically neutral neutrinos have fundamental properties that make them difficult to detect. Yet, to gather enough information to probe the most violent astrophysical sources, scientists must estimate the direction of neutrino events. If algorithms could be made considerably faster and more accurate, it would allow for more neutrino events to be analyzed, possibly even in real-time and dramatically increase the chance to identify cosmic neutrino sources. Rapid detection could enable networks of telescopes worldwide to search for more transient phenomena.**\n\n**Researchers have developed multiple approaches over the past ten years to reconstruct neutrino events. However, problems arise as existing solutions are far from perfect. They're either fast but inaccurate or more accurate at the price of huge computational costs.**\n\n**The 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\n**By 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.**","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport math\nimport time\nimport gc\nfrom tqdm import tqdm\ntqdm.pandas()","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:17:13.726055Z","iopub.execute_input":"2023-01-31T11:17:13.726891Z","iopub.status.idle":"2023-01-31T11:17:13.756932Z","shell.execute_reply.started":"2023-01-31T11:17:13.726793Z","shell.execute_reply":"2023-01-31T11:17:13.755625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MODE = 'test'\nUSE_POLARS = True\nTRAIN_MAX_EVENTS = None\nTRAIN_BATCH_START = 1\nTRAIN_N_BATCHES = 1\nUSE_UNOPTIMIZED = False\nFIND_BEST_POINTS = True\nMIN_PRIMARY_DATAPOINTS = 2\nMAX_Z = 3\nMAX_DEEP_Z = 1\nMAX_T = 350\nMAX_DEEP_T = 180\nUSE_ENSEMBLE = True\nWEIGHTS = [0.58, 0.42]\nif not USE_ENSEMBLE and not USE_POLARS:\n    ALGORITHM = 'center of charge'\n    if ALGORITHM == 'least squares':\n        USE_WEIGHTED_LEAST_SQUARES = False\n\n        \nINPUT_DIR = '/kaggle/input/icecube-neutrinos-in-deep-ice'\n\nif MODE == 'test':\n    TRAIN_MAX_EVENTS = None\nif USE_ENSEMBLE:\n    USE_WEIGHTED_LEAST_SQUARES = False","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:17:13.758786Z","iopub.execute_input":"2023-01-31T11:17:13.759164Z","iopub.status.idle":"2023-01-31T11:17:13.774973Z","shell.execute_reply.started":"2023-01-31T11:17:13.759127Z","shell.execute_reply":"2023-01-31T11:17:13.773536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if USE_POLARS:\n    try:\n        import polars as pl\n    except:\n        print('Installing polars, please wait about 35 seconds...')\n        !pip install /kaggle/input/polars01516/polars-0.15.16-cp37-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl\n        import polars as pl","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:17:13.776346Z","iopub.execute_input":"2023-01-31T11:17:13.776685Z","iopub.status.idle":"2023-01-31T11:17:28.960174Z","shell.execute_reply.started":"2023-01-31T11:17:13.776656Z","shell.execute_reply":"2023-01-31T11:17:28.958780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def angular_dist_score(az_true, zen_true, az_pred, zen_pred):\n    if not (np.all(np.isfinite(az_true)) and\n            np.all(np.isfinite(zen_true)) and\n            np.all(np.isfinite(az_pred)) and\n            np.all(np.isfinite(zen_pred))):\n        raise ValueError(\"All arguments must be finite\")\n    sa1 = np.sin(az_true)\n    ca1 = np.cos(az_true)\n    sz1 = np.sin(zen_true)\n    cz1 = np.cos(zen_true)\n    sa2 = np.sin(az_pred)\n    ca2 = np.cos(az_pred)\n    sz2 = np.sin(zen_pred)\n    cz2 = np.cos(zen_pred)\n    scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n    scalar_prod =  np.clip(scalar_prod, -1, 1)\n    return np.average(np.abs(np.arccos(scalar_prod)))","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:17:28.962385Z","iopub.execute_input":"2023-01-31T11:17:28.962769Z","iopub.status.idle":"2023-01-31T11:17:28.973614Z","shell.execute_reply.started":"2023-01-31T11:17:28.962733Z","shell.execute_reply":"2023-01-31T11:17:28.971579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def angles_from_vectors(vectors):\n    v_squared = np.square(vectors)\n    xy_sq = np.sum(v_squared[:, 0:2], axis=1)\n    xy_d = np.sqrt(xy_sq)[:, None]\n    np.seterr(divide='ignore', invalid='ignore') \n    vectors[:, 0:2] = np.where(xy_d == 0, xy_d, vectors[:, 0:2]/xy_d)\n    d = np.sqrt(xy_sq + v_squared[:, 2])\n    vectors[:, 2] = np.where(d == 0, d, vectors[:, 2]/d)\n    np.seterr(divide='warn', invalid='warn')\n    vectors =  np.clip(vectors, -1, 1)\n    azimuth = np.arccos(vectors[:, 0])\n    azimuth = np.where(vectors[:, 1] >= 0, azimuth, 2*math.pi - azimuth)\n    azimuth = np.where(np.isfinite(azimuth), azimuth, 0.0)\n    zenith = np.arccos(vectors[:, 2])\n    zenith = np.where(np.isfinite(zenith), zenith, math.pi/2)\n    return np.stack([azimuth, zenith], axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:23:58.476468Z","iopub.execute_input":"2023-01-31T11:23:58.476979Z","iopub.status.idle":"2023-01-31T11:23:58.488887Z","shell.execute_reply.started":"2023-01-31T11:23:58.476941Z","shell.execute_reply":"2023-01-31T11:23:58.486765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def average_angles(az_list, zen_list, weights=None):\n    assert(len(az_list) == len(zen_list))\n    total = az_list[0].shape[0]\n    x = np.zeros(total)\n    y = np.zeros(total)\n    z = np.zeros(total)\n    for i in range(len(az_list)):\n        w = 1\n        if weights is not None:\n            w = weights[i]\n        az = az_list[i]\n        zen = zen_list[i]\n        assert(az.shape[0] == total)\n        assert(zen.shape[0] == total)\n        if not (np.all(np.isfinite(az)) and\n                np.all(np.isfinite(zen))):\n            raise ValueError(\"All arguments must be finite\")\n        sz = np.sin(zen)\n        x += w*np.cos(az)*sz\n        y += w*np.sin(az)*sz\n        z += w*np.cos(zen)\n    tot_w = len(az_list)\n    if weights is not None:\n        tot_w = sum(weights)\n    x = x / tot_w\n    y = y / tot_w\n    z = z / tot_w\n    d = np.sqrt(np.square(x) + np.square(y) + np.square(z))\n    x = x / d\n    y = y / d\n    z = z / d\n    return angles_from_vectors(np.stack([x, y, z], axis=1))","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:24:26.693427Z","iopub.execute_input":"2023-01-31T11:24:26.693875Z","iopub.status.idle":"2023-01-31T11:24:26.707114Z","shell.execute_reply.started":"2023-01-31T11:24:26.693837Z","shell.execute_reply":"2023-01-31T11:24:26.705246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def center_of_charge(batch):\n    batch['ev_t_min'] = batch.groupby('event_id')['time'].transform('min')\n    batch['ev_t_max'] = batch.groupby('event_id')['time'].transform('max')\n    batch['w1'] = batch.charge * (batch.time - batch.ev_t_min) / (batch.ev_t_max - batch.ev_t_min)\n    batch['w0'] = batch.charge - batch.w1\n    batch['wx0'] = batch.x * batch.w0\n    batch['wy0'] = batch.y * batch.w0\n    batch['wz0'] = batch.z * batch.w0\n    batch['wx1'] = batch.x * batch.w1\n    batch['wy1'] = batch.y * batch.w1\n    batch['wz1'] = batch.z * batch.w1\n    df = batch[['w0', 'w1', 'wx0', 'wy0', 'wz0', 'wx1', 'wy1', 'wz1']]\n    df = df.groupby('event_id').sum()\n    df[['wx0', 'wy0', 'wz0']] = df[['wx0', 'wy0', 'wz0']].div(df.w0, axis=0)\n    df[['wx1', 'wy1', 'wz1']] = df[['wx1', 'wy1', 'wz1']].div(df.w1, axis=0)\n    df[['x', 'y', 'z']] = df[['wx0', 'wy0', 'wz0']].values - df[['wx1', 'wy1', 'wz1']].values\n    df = df[['x', 'y', 'z']]\n    df[['azimuth', 'zenith']] = angles_from_vectors(df.values)\n    return(df[['azimuth', 'zenith']])","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:26:29.288230Z","iopub.execute_input":"2023-01-31T11:26:29.288711Z","iopub.status.idle":"2023-01-31T11:26:29.301283Z","shell.execute_reply.started":"2023-01-31T11:26:29.288676Z","shell.execute_reply":"2023-01-31T11:26:29.299569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def least_squares(batch, weighted=False):\n    batch['xt'] = batch.x * batch.time\n    batch['yt'] = batch.y * batch.time\n    batch['zt'] = batch.z * batch.time\n    batch['tt'] = batch.time * batch.time\n    if weighted:\n        df = batch[['x', 'y', 'z', 'time', 'xt', 'yt', 'zt', 'tt']] * batch.charge.values[:, None]\n        df['charge'] = batch.charge\n        df = df.groupby('event_id').sum()\n        df = df.div(df.charge, axis=0)\n    else:\n        df = batch[['x', 'y', 'z', 'time', 'xt', 'yt', 'zt', 'tt']]\n        df = df.groupby('event_id').mean()\n    df[['x', 'y', 'z']] = (\n                              (df[['xt', 'yt', 'zt']].values - (df[['x', 'y', 'z']].values * df['time'].values[:, None]))\n                            / (df['tt'].values - (df.time.values * df.time.values))[:, None]\n                          )\n    df = -df[['x', 'y', 'z']]\n    df[['azimuth', 'zenith']] = angles_from_vectors(df.values)\n    return df[['azimuth', 'zenith']]","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:27:12.084127Z","iopub.execute_input":"2023-01-31T11:27:12.084610Z","iopub.status.idle":"2023-01-31T11:27:12.096370Z","shell.execute_reply.started":"2023-01-31T11:27:12.084574Z","shell.execute_reply":"2023-01-31T11:27:12.095034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def process_batch(batch_id, sensor, max_events=None):\n    print('load batch...')\n    batch = pd.read_parquet(f'{INPUT_DIR}/{MODE}/batch_{batch_id}.parquet')\n    if max_events is not None:\n        batch_i = batch.reset_index()\n        event_ids = batch_i.event_id.drop_duplicates()\n        end_index = event_ids.index[max_events]\n        batch = batch_i[:end_index].set_index('event_id')\n    print(batch.shape)\n    \n    batch = batch.reset_index().merge(sensor, how='left', on='sensor_id', left_index=False).set_index('event_id')\n\n    if FIND_BEST_POINTS:\n        batch = find_best_points_pandas(batch)\n     \n    batch['primary_count'] = batch.groupby('event_id')['auxiliary'].transform('count') - batch.groupby('event_id')['auxiliary'].transform('sum')\n    batch.loc[batch.primary_count < MIN_PRIMARY_DATAPOINTS, 'auxiliary'] = False\n    batch = batch[batch.auxiliary == False]\n    print(batch.shape)\n\n    if USE_ENSEMBLE or ALGORITHM == 'center of charge':\n        df = center_of_charge(batch)\n        if USE_ENSEMBLE:\n            df1 = df\n    if USE_ENSEMBLE or ALGORITHM == 'least squares':\n        df = least_squares(batch, weighted=USE_WEIGHTED_LEAST_SQUARES)\n    if USE_ENSEMBLE:\n        df[['azimuth', 'zenith']] = average_angles([df1.azimuth.values, df.azimuth.values],\n                                                   [df1.zenith.values, df.zenith.values], \n                                                   weights=WEIGHTS)\n    return df","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:28:47.592993Z","iopub.execute_input":"2023-01-31T11:28:47.593386Z","iopub.status.idle":"2023-01-31T11:28:47.605605Z","shell.execute_reply.started":"2023-01-31T11:28:47.593355Z","shell.execute_reply":"2023-01-31T11:28:47.604797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def proximity(e, col):\n    if e.shape[0] > 2000:\n        return e[:, col['auxiliary']]\n\n    deltas = np.abs(e[:, [col['string_id'], col['depth_id'], col['time']]] - e[:, None, [col['string_id'], col['depth_id'], col['time']]])\n    dz = deltas[:, :, 1]\n\n    dz[(dz == 0) | (deltas[:, :, 0] != 0)] = MAX_Z + MAX_DEEP_Z + 1\n\n    mask = (e[:, col['sensor_id']] < 4680)\n    mask = np.broadcast_to(mask, (mask.shape[0], mask.shape[0])).T\n    dz[mask & (deltas[:, :, 2] > MAX_T)] = MAX_Z + MAX_DEEP_Z + 1\n    mask = (e[:, col['sensor_id']] >= 4680)\n    mask = np.broadcast_to(mask, (mask.shape[0], mask.shape[0])).T\n    dz[mask & (deltas[:, :, 2] > MAX_DEEP_T)] = MAX_Z + MAX_DEEP_Z + 1\n    dz = dz.min(axis=1)\n    e[:, col['auxiliary']] = True\n    e[((e[:, col['sensor_id']] < 4680) & (dz <= MAX_Z)) | ((e[:, col['sensor_id']] >= 4680) & (dz <= MAX_DEEP_Z)), col['auxiliary']] = False\n    return e[:, col['auxiliary']]","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:36:38.407032Z","iopub.execute_input":"2023-01-31T11:36:38.407536Z","iopub.status.idle":"2023-01-31T11:36:38.421895Z","shell.execute_reply.started":"2023-01-31T11:36:38.407497Z","shell.execute_reply":"2023-01-31T11:36:38.419893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def proximity_unoptimized(df):\n    if df.shape[0] > 2000:\n        return df\n    deltas = np.abs(df[['string_id', 'depth_id', 'time']].values - df[['string_id', 'depth_id', 'time']].values[:, None, :])\n    mask = (df.sensor_id < 4680)\n    dz = deltas[:, :, 1]\n\n    dz[(dz == 0) | (deltas[:, :, 0] != 0)] = MAX_Z + MAX_DEEP_Z + 1\n\n    mask = (df.sensor_id < 4680)\n    mask = np.broadcast_to(mask, (mask.shape[0], mask.shape[0])).T\n    dz[mask & (deltas[:, :, 2] > MAX_T)] = MAX_Z + MAX_DEEP_Z + 1\n    mask = (df.sensor_id >= 4680)\n    mask = np.broadcast_to(mask, (mask.shape[0], mask.shape[0])).T\n    dz[mask & (deltas[:, :, 2] > MAX_DEEP_T)] = MAX_Z + MAX_DEEP_Z + 1\n\n    dz = dz.min(axis=1)\n    df.auxiliary = True\n    df.loc[((df.sensor_id < 4680) & (dz <= MAX_Z)) | ((df.sensor_id >= 4680) & (dz <= MAX_DEEP_Z)), 'auxiliary'] = False\n    return df","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:37:03.263499Z","iopub.execute_input":"2023-01-31T11:37:03.263970Z","iopub.status.idle":"2023-01-31T11:37:03.277527Z","shell.execute_reply.started":"2023-01-31T11:37:03.263934Z","shell.execute_reply":"2023-01-31T11:37:03.275576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def find_best_points_pandas(batch):\n    if USE_UNOPTIMIZED:\n        cols = ['sensor_id', 'time', 'auxiliary', 'string_id', 'depth_id']\n        batch[cols] = batch[cols].groupby('event_id').progress_apply(proximity_unoptimized)\n        return batch\n\n    column_to_index = { k:v for v,k in enumerate(batch.columns)}\n    events = np.split(batch.values.astype('float32'), np.unique(batch.index.values, return_index=True)[1][1:])\n\n    batch.auxiliary = np.concatenate([proximity(e, column_to_index) for e in tqdm(events)])\n    return batch","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:37:46.918417Z","iopub.execute_input":"2023-01-31T11:37:46.918864Z","iopub.status.idle":"2023-01-31T11:37:46.928545Z","shell.execute_reply.started":"2023-01-31T11:37:46.918831Z","shell.execute_reply":"2023-01-31T11:37:46.926857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if not USE_POLARS:\n    print(MODE)\n    meta = pd.read_parquet(f'{INPUT_DIR}/{MODE}_meta.parquet')\n    print(meta.shape)\n\n    print('load sensor data...')\n    sensor = pd.read_csv(f'{INPUT_DIR}/sensor_geometry.csv')\n    sensor['string_id'] = sensor.sensor_id // 60\n    sensor['depth_id'] = sensor.sensor_id % 60\n    print(sensor.shape)\n\n    if MODE == 'train':\n        batch_id_start = TRAIN_BATCH_START\n        batch_id_end = batch_id_start + TRAIN_N_BATCHES\n    else:\n        batch_id_start = meta.batch_id.values[0]\n        batch_id_end = meta.batch_id.values[-1] + 1\n\n    sub = []\n    for batch_id in range(batch_id_start,batch_id_end):\n        print(batch_id)\n        t = time.time()\n        df = process_batch(batch_id, sensor, max_events=TRAIN_MAX_EVENTS)\n        print(f'Time: {time.time() - t:0.2f}s')\n        if MODE == 'test':\n            sub.append(df)\n        else:\n            if TRAIN_MAX_EVENTS is not None:\n                print(angular_dist_score(meta[meta.batch_id == batch_id].azimuth.values[:TRAIN_MAX_EVENTS],\n                                         meta[meta.batch_id == batch_id].zenith.values[:TRAIN_MAX_EVENTS],\n                                         df.azimuth.values, df.zenith.values))\n            else:\n                print(angular_dist_score(meta[meta.batch_id == batch_id].azimuth.values, meta[meta.batch_id == batch_id].zenith.values,\n                                         df.azimuth.values, df.zenith.values))\n\n    if MODE == 'test':\n        sub = pd.concat(sub, axis=0)\n        sub.to_csv('submission.csv', index=True)\n        print(sub)\n\n    print('done')","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:38:14.993828Z","iopub.execute_input":"2023-01-31T11:38:14.994364Z","iopub.status.idle":"2023-01-31T11:38:15.007846Z","shell.execute_reply.started":"2023-01-31T11:38:14.994320Z","shell.execute_reply":"2023-01-31T11:38:15.006419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def time_weighted_centering(batch, charge_weighted=True):\n    batch = batch.with_columns([pl.col('time').min().over('event_id').alias('ev_t_min'),\n                                pl.col('time').max().over('event_id').alias('ev_t_max')])\n    if charge_weighted:\n        batch = batch.with_columns((pl.col('charge') * (pl.col('time') - pl.col('ev_t_min'))\n                                    / (pl.col('ev_t_max') - pl.col('ev_t_min'))).alias('w1'))\n        batch = batch.with_columns((pl.col('charge') - pl.col('w1')).alias('w0'))\n    else:\n        batch = batch.with_columns(((pl.col('time') - pl.col('ev_t_min'))\n                                    / (pl.col('ev_t_max') - pl.col('ev_t_min'))).alias('w1'))\n        batch = batch.with_columns((pl.lit(1) - pl.col('w1')).alias('w0'))\n\n    batch = batch.select(\n        [\n            pl.col('event_id'),\n            pl.col('w0'),\n            pl.col('w1'),\n            (pl.col('x') * pl.col('w0')).alias('wx0'),\n            (pl.col('y') * pl.col('w0')).alias('wy0'),\n            (pl.col('z') * pl.col('w0')).alias('wz0'),\n            (pl.col('x') * pl.col('w1')).alias('wx1'),\n            (pl.col('y') * pl.col('w1')).alias('wy1'),\n            (pl.col('z') * pl.col('w1')).alias('wz1'),\n        ]\n    ).collect().groupby('event_id', maintain_order=True).sum()\n\n    batch_values = batch.select(\n        [\n            ((pl.col('wx0') / pl.col('w0')) - (pl.col('wx1') / pl.col('w1'))).alias('x'),\n            ((pl.col('wy0') / pl.col('w0')) - (pl.col('wy1') / pl.col('w1'))).alias('y'),\n            ((pl.col('wz0') / pl.col('w0')) - (pl.col('wz1') / pl.col('w1'))).alias('z'),\n        ]\n    ).to_numpy()\n    return angles_from_vectors(batch_values), batch","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:38:59.027678Z","iopub.execute_input":"2023-01-31T11:38:59.028133Z","iopub.status.idle":"2023-01-31T11:38:59.044150Z","shell.execute_reply.started":"2023-01-31T11:38:59.028098Z","shell.execute_reply":"2023-01-31T11:38:59.042522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def find_best_points_polars(batch):\n    df = batch.select([pl.col('sensor_id'), pl.col('time'), pl.col('auxiliary'), pl.col('string_id'), pl.col('depth_id')])\n    column_to_index = { k:v for v,k in enumerate(df.columns)}\n    batch_values = df.to_numpy().astype('float32')\n\n    events = np.split(batch_values, np.unique(batch.select(pl.col('event_id')), return_index=True)[1][1:])\n\n    batch = batch.with_columns(pl.Series(np.concatenate([proximity(e, column_to_index) for e in tqdm(events)]).astype('bool')).alias('auxiliary'))\n\n    return batch","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:39:37.217840Z","iopub.execute_input":"2023-01-31T11:39:37.218855Z","iopub.status.idle":"2023-01-31T11:39:37.228457Z","shell.execute_reply.started":"2023-01-31T11:39:37.218808Z","shell.execute_reply":"2023-01-31T11:39:37.227057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if USE_POLARS:\n    print(MODE)\n    meta = pl.scan_parquet(f'{INPUT_DIR}/{MODE}_meta.parquet')\n\n    print('load sensor data...')\n    sensor = (pl.scan_csv(f'{INPUT_DIR}/sensor_geometry.csv')\n                .with_columns([\n                    pl.col('sensor_id').cast(pl.Int16),\n                    (pl.col('sensor_id') // 60).alias('string_id'),\n                    (pl.col('sensor_id') % 60).alias('depth_id'),                \n                ])\n             )\n\n    print(sensor)\n\n    if MODE == 'train':\n        batch_id_start = TRAIN_BATCH_START\n        batch_id_end = batch_id_start + TRAIN_N_BATCHES\n    else:\n        batch_id_start = meta.select(pl.col('batch_id')).collect()[0, 0]\n        batch_id_end = meta.select(pl.col('batch_id')).collect()[-1, 0] + 1\n\n    print(batch_id_start, batch_id_end)","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:40:00.960338Z","iopub.execute_input":"2023-01-31T11:40:00.960757Z","iopub.status.idle":"2023-01-31T11:40:00.986863Z","shell.execute_reply.started":"2023-01-31T11:40:00.960710Z","shell.execute_reply":"2023-01-31T11:40:00.985971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if USE_POLARS:\n    sub = []\n    for batch_id in range(batch_id_start,batch_id_end):\n        print(batch_id)\n        t = time.time()\n        max_events=TRAIN_MAX_EVENTS\n        print('load batch...')\n        batch = pl.scan_parquet(f'{INPUT_DIR}/{MODE}/batch_{batch_id}.parquet')\n\n        if max_events is not None:\n            batch = batch.collect()\n            last_event_id = batch.select(pl.col('event_id')).unique()[TRAIN_MAX_EVENTS-1, 0]\n            batch = batch.lazy().filter(pl.col('event_id') <= last_event_id)\n\n        batch = batch.join(sensor, on='sensor_id', how='left').collect()\n \n        if FIND_BEST_POINTS:\n            batch = find_best_points_polars(batch)\n\n        batch = batch.lazy().with_columns((pl.col('auxiliary').count().over('event_id') - pl.col('auxiliary').sum().over('event_id')).alias('primary_count'))\n        batch = batch.filter((~pl.col('auxiliary')) | (pl.col('primary_count') < MIN_PRIMARY_DATAPOINTS))\n\n        preds, events = time_weighted_centering(batch)\n        if USE_ENSEMBLE:\n            preds1 = preds\n            preds, events = time_weighted_centering(batch, charge_weighted=False)\n            preds = average_angles([preds1[:, 0], preds[:, 0]], [preds1[:, 1], preds[:, 1]], weights=WEIGHTS)\n\n        if MODE == 'test':\n            sub.append(events.select([pl.col('event_id'), pl.Series(preds[:, 0]).alias('azimuth'),\n                                                         pl.Series(preds[:, 1]).alias('zenith')]))\n        else:\n            meta = meta.filter((pl.col('batch_id') >= batch_id_start) & (pl.col('batch_id') < batch_id_end))\n            if isinstance(meta, pl.LazyFrame):\n                meta = meta.collect()\n            meta_values = meta.filter(pl.col('batch_id') == batch_id).select([pl.col('azimuth'), pl.col('zenith')]).to_numpy()\n            if TRAIN_MAX_EVENTS is not None:\n                print(angular_dist_score(meta_values[:TRAIN_MAX_EVENTS, 0], meta_values[:TRAIN_MAX_EVENTS, 1], preds[:, 0], preds[:, 1]))\n            else:\n                print(angular_dist_score(meta_values[:, 0], meta_values[:, 1], preds[:, 0], preds[:, 1]))\n        print(f'Time: {time.time() - t:0.2f}s')\n\n    if MODE == 'test':\n        sub = pl.concat(sub)\n        sub.write_csv('submission.csv')\n        print(sub)\n\n    print('done')","metadata":{"execution":{"iopub.status.busy":"2023-01-31T11:40:46.179609Z","iopub.execute_input":"2023-01-31T11:40:46.180860Z","iopub.status.idle":"2023-01-31T11:40:46.256999Z","shell.execute_reply.started":"2023-01-31T11:40:46.180804Z","shell.execute_reply":"2023-01-31T11:40:46.256030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**In this competition, I'm identifying which direction neutrinos detected by the IceCube neutrino observatory came from. When detection events can be localized quickly enough, traditional telescopes are recruited to investigate short-lived neutrino sources such as supernovae or gamma ray bursts. Because the sky is huge better localization will not only associate neutrinos with sources but also to help partner observatories limit their search space. With an average of three thousand events per second to process, it's difficult to keep up with the stream of data using traditional methods.**","metadata":{}}]}