{"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":"# Simple PCA line prediction\n\nThis notebook uses PCA to predict the line of best fit to the data. It uses a single PCA component to predict the line of best fit to the data and reducing the dimensionality of the data to 1.","metadata":{}},{"cell_type":"code","source":"import time\n\nstart_time = time.time()","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.683317Z","iopub.execute_input":"2023-02-25T13:05:15.684077Z","iopub.status.idle":"2023-02-25T13:05:15.689222Z","shell.execute_reply.started":"2023-02-25T13:05:15.684038Z","shell.execute_reply":"2023-02-25T13:05:15.687886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd \nimport os\nimport random\nfrom typing import Tuple\nimport numpy as np\nimport math\nimport pyarrow.parquet as pq\n","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.799039Z","iopub.execute_input":"2023-02-25T13:05:15.799461Z","iopub.status.idle":"2023-02-25T13:05:15.804990Z","shell.execute_reply.started":"2023-02-25T13:05:15.799423Z","shell.execute_reply":"2023-02-25T13:05:15.804170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with np.errstate(divide='ignore'):\n    np.float64(1.0) / 0.0","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.807120Z","iopub.execute_input":"2023-02-25T13:05:15.807455Z","iopub.status.idle":"2023-02-25T13:05:15.817731Z","shell.execute_reply.started":"2023-02-25T13:05:15.807424Z","shell.execute_reply":"2023-02-25T13:05:15.816900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/icecube-neutrinos-in-deep-ice\"\nBATCH = 1\nEVENT = 24\nMAX_SEQUENCE_LENGTH=100\nEXCLUDE_AUXILIARY= True\nSHOW_DF=False\nIS_TRAINING = False\nSET = 'train' if IS_TRAINING else 'test'","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.819018Z","iopub.execute_input":"2023-02-25T13:05:15.820075Z","iopub.status.idle":"2023-02-25T13:05:15.830536Z","shell.execute_reply.started":"2023-02-25T13:05:15.820040Z","shell.execute_reply":"2023-02-25T13:05:15.829375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def seed_it_all(seed=7):\n    \"\"\" Attempt to be Reproducible \"\"\"\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    random.seed(seed)\n    np.random.seed(seed)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.831648Z","iopub.execute_input":"2023-02-25T13:05:15.832617Z","iopub.status.idle":"2023-02-25T13:05:15.843111Z","shell.execute_reply.started":"2023-02-25T13:05:15.832573Z","shell.execute_reply":"2023-02-25T13:05:15.842190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"seed_it_all(10)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.846380Z","iopub.execute_input":"2023-02-25T13:05:15.847171Z","iopub.status.idle":"2023-02-25T13:05:15.855958Z","shell.execute_reply.started":"2023-02-25T13:05:15.847120Z","shell.execute_reply":"2023-02-25T13:05:15.854896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import importlib\n\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2\n    for col in df.columns:\n        print(f'Optimising col {col}')\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)\n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n\n    end_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))\n    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))\n\n    return df\n\ndef import_data(file: str):\n    \"\"\"create a dataframe and optimize its memory usage\"\"\"\n    function_name = f\"read_{file.split('.')[-1]}\"\n    function = getattr(pd, function_name)\n    df = function(file)\n    df = reduce_mem_usage(df)\n    return df\n\n","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.857925Z","iopub.execute_input":"2023-02-25T13:05:15.858867Z","iopub.status.idle":"2023-02-25T13:05:15.876187Z","shell.execute_reply.started":"2023-02-25T13:05:15.858801Z","shell.execute_reply":"2023-02-25T13:05:15.875300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry = import_data(f'{DATA_DIR}/sensor_geometry.csv')","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.877484Z","iopub.execute_input":"2023-02-25T13:05:15.878185Z","iopub.status.idle":"2023-02-25T13:05:15.905483Z","shell.execute_reply.started":"2023-02-25T13:05:15.878148Z","shell.execute_reply":"2023-02-25T13:05:15.904306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### For Geometry","metadata":{}},{"cell_type":"code","source":"from typing import List, Tuple\nfrom sklearn.decomposition import PCA\n\ndef get_direction(coords: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Get the direction vector from a list of coordinates.\n    \"\"\"\n    pca = PCA(n_components=1)\n    pca.fit(coords) \n    direction_vector = pca.components_#type: ignore\n    return direction_vector","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.906959Z","iopub.execute_input":"2023-02-25T13:05:15.907320Z","iopub.status.idle":"2023-02-25T13:05:15.913735Z","shell.execute_reply.started":"2023-02-25T13:05:15.907285Z","shell.execute_reply":"2023-02-25T13:05:15.912420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_line(origin, vector, extent):\n    below_origin = origin - direction_vector * extent\n    above_origin = origin + direction_vector * extent\n    line = np.vstack((below_origin, above_origin))\n    return line","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.917010Z","iopub.execute_input":"2023-02-25T13:05:15.917770Z","iopub.status.idle":"2023-02-25T13:05:15.928412Z","shell.execute_reply.started":"2023-02-25T13:05:15.917722Z","shell.execute_reply":"2023-02-25T13:05:15.927212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef cartesian_to_sphere(x: float, y: float, z: float) -> Tuple[float, float]:\n    \"\"\"Maps vector cartesian coordinates (x, y, z) from the origin to spherical angles azimuth and zenith.\n    \n    See: https://en.wikipedia.org/wiki/Spherical_coordinate_system\n\n    Args:\n        x (float): The x-coordinate of the point.\n        y (float): The y-coordinate of the point.\n        z (float): The z-coordinate of the point.\n\n    Returns:\n        tuple[float, float]: The azimuth and zenith angles in radians.\n    \"\"\"\n    x2y2 = x**2 + y**2\n    r = math.sqrt(x2y2 + z**2)\n    \n    sqt_x2y2 = math.sqrt(x2y2)\n    x_dv_py = x / sqt_x2y2 if sqt_x2y2 != 0 else 0\n    azimuth = math.acos(x_dv_py) * np.sign(y)\n\n    zenith = math.acos(z / r) if r != 0 else 0\n\n    return azimuth, zenith\n\n\ndef sphere_to_cartesian(azimuth: float, zenith: float) -> Tuple[float, float, float]:\n    \"\"\"Map spherical coordinates to cartesian coordinates.\n    see: https://stackoverflow.com/a/10868220/4521646\n    \n    Args:\n        azimuth (float): The azimuth angle in radians.\n        zenith (float): The zenith angle in radians.\n\n    Returns:\n        tuple: The x, y, z vector cartesian coordinates of the point from the origin.\n    \"\"\"\n    x = math.sin(zenith) * math.cos(azimuth)\n    y = math.sin(zenith) * math.sin(azimuth)\n    z = math.cos(zenith)\n    return x, y, z\n\n\ndef adjust_sphere(azimuth:float, zenith:float) -> Tuple[float, float]:\n    \"\"\"Adjust azimuth and zenith to be within [-pi, pi]\n\n    Args:\n        azimuth (float): The azimuth to adjust\n        zenith (float): The zenith to adjust\n\n    Returns:\n        float: The adjusted azimuth and zenith\n    \"\"\"\n    print('adjust_sphere takes',azimuth, zenith)\n    \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","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.930171Z","iopub.execute_input":"2023-02-25T13:05:15.930638Z","iopub.status.idle":"2023-02-25T13:05:15.945199Z","shell.execute_reply.started":"2023-02-25T13:05:15.930585Z","shell.execute_reply":"2023-02-25T13:05:15.943792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### For scoring","metadata":{}},{"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    The lower the angle, the more similar the two vectors are meaning the score is better.\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 numerical instability\n    # that might otherwise occur 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)))\n","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.947351Z","iopub.execute_input":"2023-02-25T13:05:15.948110Z","iopub.status.idle":"2023-02-25T13:05:15.974994Z","shell.execute_reply.started":"2023-02-25T13:05:15.948060Z","shell.execute_reply":"2023-02-25T13:05:15.973950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Make prediction on test set","metadata":{}},{"cell_type":"code","source":"from dask import dataframe as dd","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.979182Z","iopub.execute_input":"2023-02-25T13:05:15.980191Z","iopub.status.idle":"2023-02-25T13:05:15.992088Z","shell.execute_reply.started":"2023-02-25T13:05:15.980149Z","shell.execute_reply":"2023-02-25T13:05:15.990778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create a dataset object for the Parquet file\n\nmeta_dfd = dd.read_parquet(f'{DATA_DIR}/{SET}_meta.parquet', \n    blocksize=64000000 # = 64 Mb chunks\n)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:15.993964Z","iopub.execute_input":"2023-02-25T13:05:15.994377Z","iopub.status.idle":"2023-02-25T13:05:16.016508Z","shell.execute_reply.started":"2023-02-25T13:05:15.994342Z","shell.execute_reply":"2023-02-25T13:05:16.015317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_ids = meta_dfd['batch_id'].unique().compute().values","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.017967Z","iopub.execute_input":"2023-02-25T13:05:16.018995Z","iopub.status.idle":"2023-02-25T13:05:16.032334Z","shell.execute_reply.started":"2023-02-25T13:05:16.018958Z","shell.execute_reply":"2023-02-25T13:05:16.031201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_event_df(batch_df: dd.DataFrame, sensor_geometry: pd.DataFrame, event_id: str) -> pd.DataFrame:\n    \"\"\"\n    Get a DataFrame for a specific event.\n\n    Parameters:\n    train_batch_df (pandas.DataFrame): The batch DataFrame.\n    sensor_geometry (pandas.DataFrame): The sensor geometry DataFrame.\n    event_id (str): The event identifier.\n\n    Returns:\n    pandas.DataFrame: A DataFrame containing data for the specified event.\n    \"\"\"\n    event_df = batch_df[batch_df['event_id'] == event_id].compute()[:MAX_SEQUENCE_LENGTH]\n    event_df = pd.merge(\n        left=event_df,\n        right=sensor_geometry,\n        how='inner',\n        on='sensor_id'\n    )\n    return event_df","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.033978Z","iopub.execute_input":"2023-02-25T13:05:16.034541Z","iopub.status.idle":"2023-02-25T13:05:16.040710Z","shell.execute_reply.started":"2023-02-25T13:05:16.034504Z","shell.execute_reply":"2023-02-25T13:05:16.039865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(columns=['event_id', 'azimuth', 'zenith'])\nav_batch_time = None\nav_event_time = None\nevent_count = 0\nbatch_count = 0\n\nfor batch_id in batch_ids:\n    batch_start_time = time.time()\n    \n    batch_dfd = dd.read_parquet(f'{DATA_DIR}/{SET}/batch_{batch_id}.parquet', \n        blocksize=64000000 # = 64 Mb chunks,\n    ).reset_index()\n\n    \n    if EXCLUDE_AUXILIARY:\n        batch_dfd = batch_dfd[~batch_dfd['auxiliary']]\n        \n    # Loop through unique event IDs\n    for event_id in batch_dfd['event_id'].unique().compute().values:\n        \n        print('Processing event', event_id, ' in batch ', batch_id)\n\n        event_df = get_event_df(batch_dfd, sensor_geometry, event_id)\n        print('Event df length', len(event_df))\n\n        x, y, z = event_df['x'], event_df['y'], event_df['z']\n        coords = np.array((x,y,z)).T\n        direction_vector = get_direction(coords)\n        origin = np.mean(coords, axis=0)\n        euclidean_distance = np.linalg.norm(coords - origin, axis=1)\n        extent = np.nanmax(np.nan_to_num(euclidean_distance, copy=False, nan=0.0, posinf=0.0, neginf=0.0))\n        line = get_line(origin, direction_vector, extent)\n        x2, y2, z2= line[:,0], line[:,1], line[:,2]\n        x2 = x2.max() - x2.min()\n        y2 = y2.max() - y2.min()\n        z2 = z2.max() - z2.min()\n        az_pred, zen_pred = adjust_sphere(*cartesian_to_sphere(x2,y2,z2))\n        new_row = pd.DataFrame({ 'event_id': [event_id], 'azimuth': az_pred,'zenith': zen_pred})\n        submission = reduce_mem_usage(pd.concat([ new_row, submission.loc[:]]))\n        current_time = time.time() \n        print(\"Total time taken so far : \", round((current_time - start_time)/60), \"minutes\")\n        \n    batch_process_time =  time.time() - batch_start_time\n    print(\"Time for batch: \", round((batch_process_time)/60), \"minutes\")\n    if av_batch_time is None:\n        av_batch_time = batch_process_time\n    else:\n        av_batch_time = (batch_process_time + av_batch_time) / 2\n        \n    print('Average batch time', av_batch_time)\n        \n        ","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.042268Z","iopub.execute_input":"2023-02-25T13:05:16.042699Z","iopub.status.idle":"2023-02-25T13:05:16.147694Z","shell.execute_reply.started":"2023-02-25T13:05:16.042655Z","shell.execute_reply":"2023-02-25T13:05:16.146502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = submission.sort_values([\"event_id\"]).reset_index(drop=True)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.149244Z","iopub.execute_input":"2023-02-25T13:05:16.149678Z","iopub.status.idle":"2023-02-25T13:05:16.168605Z","shell.execute_reply.started":"2023-02-25T13:05:16.149640Z","shell.execute_reply":"2023-02-25T13:05:16.167422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.169926Z","iopub.execute_input":"2023-02-25T13:05:16.170275Z","iopub.status.idle":"2023-02-25T13:05:16.178297Z","shell.execute_reply.started":"2023-02-25T13:05:16.170244Z","shell.execute_reply":"2023-02-25T13:05:16.177027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"end_time = time.time()\ntotal_time = end_time - start_time\nprint(\"Total time taken: \", total_time, \"seconds\")","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.179812Z","iopub.execute_input":"2023-02-25T13:05:16.180187Z","iopub.status.idle":"2023-02-25T13:05:16.191368Z","shell.execute_reply.started":"2023-02-25T13:05:16.180154Z","shell.execute_reply":"2023-02-25T13:05:16.190144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submission.to_csv('submission.csv', index=True)\n!head submission.csv","metadata":{"execution":{"iopub.status.busy":"2023-02-25T13:05:16.194460Z","iopub.execute_input":"2023-02-25T13:05:16.194799Z","iopub.status.idle":"2023-02-25T13:05:17.314261Z","shell.execute_reply.started":"2023-02-25T13:05:16.194767Z","shell.execute_reply":"2023-02-25T13:05:17.312846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}