{"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":"# 3D line fitting + computed vector direction choice\nIn this notebook I will present a way of imporving 3D line fitting by using time values. Using this way I improved my LB score from `1.26` to `1.193`. With chosing a proper direciton of predicted vector the score could be around `0.8`.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nfrom pathlib import Path\nimport math\nfrom tqdm import tqdm\nfrom random import sample","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-29T13:26:35.980177Z","iopub.execute_input":"2023-01-29T13:26:35.980666Z","iopub.status.idle":"2023-01-29T13:26:36.009655Z","shell.execute_reply.started":"2023-01-29T13:26:35.980565Z","shell.execute_reply":"2023-01-29T13:26:36.008541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subset = 'train'\nlimit = 10_000\n\ndataset_path = Path(r'/kaggle/input/icecube-neutrinos-in-deep-ice')\ntest_meta_data = pd.read_parquet(dataset_path / f'{subset}_meta.parquet')\nif limit is not None:\n    test_meta_data = test_meta_data[:limit]\nsensor_geo = pd.read_csv(dataset_path / 'sensor_geometry.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-29T13:26:36.011379Z","iopub.execute_input":"2023-01-29T13:26:36.012195Z","iopub.status.idle":"2023-01-29T13:27:35.022306Z","shell.execute_reply.started":"2023-01-29T13:26:36.012161Z","shell.execute_reply":"2023-01-29T13:27:35.019228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def angular_dist_score(az_true, zen_true, az_pred, zen_pred):\n    '''\n    calculate the MAE of the angular distance between two directions.\n    The two vectors are first converted to cartesian unit vectors,\n    and then their scalar product is computed, which is equal to\n    the cosine of the angle between the two vectors. The inverse \n    cosine (arccos) thereof is then the angle between the two input vectors\n    \n    Parameters:\n    -----------\n    \n    az_true : float (or array thereof)\n        true azimuth value(s) in radian\n    zen_true : float (or array thereof)\n        true zenith value(s) in radian\n    az_pred : float (or array thereof)\n        predicted azimuth value(s) in radian\n    zen_pred : float (or array thereof)\n        predicted zenith value(s) in radian\n    \n    Returns:\n    --------\n    \n    dist : float\n        mean over the angular distance(s) in radian\n    '''\n    \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    \n    # pre-compute all sine and cosine values\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    \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    \n    # scalar product of the two cartesian vectors (x = sz*ca, y = sz*sa, z = cz)\n    scalar_prod = sz1*sz2*(ca1*ca2 + sa1*sa2) + (cz1*cz2)\n    \n    # scalar product of two unit vectors is always between -1 and 1, this is against nummerical instability\n    # that might otherwise occure from the finite precision of the sine and cosine functions\n    scalar_prod =  np.clip(scalar_prod, -1, 1)\n    \n    # convert back to an angle (in radian)\n    return np.average(np.abs(np.arccos(scalar_prod)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-29T13:27:35.033901Z","iopub.execute_input":"2023-01-29T13:27:35.035128Z","iopub.status.idle":"2023-01-29T13:27:35.053138Z","shell.execute_reply.started":"2023-01-29T13:27:35.035062Z","shell.execute_reply":"2023-01-29T13:27:35.052205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _cartesian_to_spherical(x, y, z):\n    r = math.sqrt(x * x + y * y + z * z)\n    c = -1 if y < 0 else 1\n    azimuth = math.acos(x / math.sqrt(x * x + y * y)) * c\n    zenith = math.acos(z / r)\n    return azimuth, zenith\n\ndef _adjust_spherical(azimuth, zenith):\n    if azimuth < 0 and zenith < 0:\n        azimuth *= -1\n        zenith *= -1\n    if azimuth < 0:\n        azimuth += math.pi * 2\n    elif zenith < 0:\n        zenith += math.pi\n    return azimuth, zenith\n\ndef cartesian_to_spherical(pred_direction: np.ndarray):\n    try:\n        a, z = _adjust_spherical(*_cartesian_to_spherical(*(pred_direction.tolist())))\n        if np.isnan(a):\n            a = 0.0\n        if np.isnan(z):\n            z = 0.0\n        return a, z\n    except:\n        return 0.0, 0.0","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-29T13:27:35.056918Z","iopub.execute_input":"2023-01-29T13:27:35.059173Z","iopub.status.idle":"2023-01-29T13:27:35.076350Z","shell.execute_reply.started":"2023-01-29T13:27:35.059115Z","shell.execute_reply":"2023-01-29T13:27:35.074720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"total_score = {'svd': 0.0, 'svd+t': 0.0, 'best': 0.0}\n\nperc = 0.25\nprev_batch_id = -1\n\nprogress_bar = tqdm(range(len(test_meta_data)))\nfor i in progress_bar:\n    meta_data_row = test_meta_data.iloc[i]\n    batch_id = int(meta_data_row['batch_id'])\n    first_pulse_index = int(meta_data_row['first_pulse_index'])\n    last_pulse_index = int(meta_data_row['last_pulse_index'])\n    \n    if batch_id != prev_batch_id:\n        train_batch = pd.read_parquet(dataset_path / (f'{subset}/batch_{batch_id}.parquet'))\n        prev_batch_id = batch_id\n\n    event_data = train_batch.iloc[first_pulse_index:last_pulse_index + 1]\n    event_data = event_data[~event_data['auxiliary']]\n    if len(event_data) > 800:\n        event_data = event_data.sort_values('charge', ascending=False)[:800]\n\n    event_data_with_pos = pd.merge(event_data, sensor_geo, on='sensor_id', how='left')\n\n    time_sorted = event_data_with_pos.sort_values('time', ascending=True)\n    m = max(1, int(perc*len(time_sorted)))\n    first = time_sorted.iloc[:m]\n    last = time_sorted.iloc[-m:]\n    \n    first = np.array([first['x'].mean(), first['y'].mean(), first['z'].mean()])\n    last = np.array([last['x'].mean(), last['y'].mean(), last['z'].mean()])\n\n    time_dir = last - first\n    time_dir /= np.linalg.norm(time_dir)\n\n    data = np.concatenate(\n        (np.array(event_data_with_pos['x'])[:, np.newaxis], \n         np.array(event_data_with_pos['y'])[:, np.newaxis], \n         np.array(event_data_with_pos['z'])[:, np.newaxis]), \n        axis=1\n    )\n\n    _, _, v = np.linalg.svd(data - data.mean(axis=0))\n\n    pred_direction = -v[0]\n    pred_direction_neg = -pred_direction.copy()\n    \n    # check if predicted direction is aligned with vector computed based on time\n    if np.dot(time_dir, pred_direction) > 0:\n        pred_direction_t = -pred_direction.copy()\n    else:\n        pred_direction_t = pred_direction.copy()\n    \n    azimuth, zenith = cartesian_to_spherical(pred_direction)\n    azimuth_t, zenith_t = cartesian_to_spherical(pred_direction_t)\n    azimuth_neg, zenith_neg = cartesian_to_spherical(pred_direction_neg)\n    \n    if subset == 'train':\n        true_azimuth = meta_data_row['azimuth']\n        true_zenith = meta_data_row['zenith']\n        \n        score = angular_dist_score(\n            true_azimuth, true_zenith,\n            azimuth, zenith\n        )\n        score_adj = angular_dist_score(\n            true_azimuth, true_zenith,\n            azimuth_t, zenith_t\n        )\n        score_neg = angular_dist_score(\n            true_azimuth, true_zenith,\n            azimuth_neg, zenith_neg\n        )\n        total_score['svd'] += score\n        total_score['svd+t'] += score_adj\n        total_score['best'] += min(score, score_neg)\n    \n        progress_bar.set_description(\n            f'svd: {total_score[\"svd\"]/(i + 1):.3f}, '\n            f'svd+t: {total_score[\"svd+t\"]/(i + 1):.3f}, '\n            f'best: {total_score[\"best\"]/(i + 1):.3f} | '\n        )\ntotal_score['svd'] /= len(test_meta_data)\ntotal_score['svd+t'] /= len(test_meta_data)\ntotal_score['best'] /= len(test_meta_data)","metadata":{"execution":{"iopub.status.busy":"2023-01-29T13:27:35.078438Z","iopub.execute_input":"2023-01-29T13:27:35.079317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result = pd.DataFrame(\n    data=total_score.values(),\n    columns=['score'],\n    index=total_score.keys()\n)\nresult.index.name = 'method'","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(result)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}