{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":102335,"databundleVersionId":12518947,"sourceType":"competition"},{"sourceId":12144415,"sourceType":"datasetVersion","datasetId":7648709}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Intro\n","metadata":{}},{"cell_type":"markdown","source":"\n## The goal\nThe goal is is to extract as many features from IMU as possible\n\n## On giants shoulders\nfeature ingeneering inspyred by this projects:\n\n- [CMI3 Pyroch Baseline+model_add, Aug, folds](https://www.kaggle.com/code/myso1987/cmi3-pyroch-baseline-model-add-aug-folds)\n\n- [LB 0.82 | 5fold single bert model](https://www.kaggle.com/code/wasupandceacar/lb-0-82-5fold-single-bert-model)\n\n- [cmi imu-only lstm](https://www.kaggle.com/code/jiazhuang/cmi-imu-only-lstm)","metadata":{}},{"cell_type":"markdown","source":"# Features description","metadata":{}},{"cell_type":"markdown","source":"\nadditional features description:\n- `[acc/gyro]_mag` (2 features): norms (magnitudes) of `linear acceleration` and `orientation` data to calculate only strength of this vectors, ignoring direction\n\n- `linear_acc` & `linear_acc_mag` (4 features): linear acceleration (without gravity) and its magnitude\n\n- `angular_vel_[x/y/z/mag]` (4 features): angular velocity (3 directions) and its magnitude\n\n- `[acc/gyro/linear_acc/angular_vel]_jerk_[x/y/z/mag]` (16 features): approximate derivatives of acceleration, linear acceleration, orientation and angular velocity\n\n- `[acc/gyro/linear_acc/angular_vel]_pow_[x/y/z/mag]` (16 features): squares of according values. Called energy cause $E = \\frac{mv^2}{2}$\n\n- `[acc/gyro/linear_acc/angular_vel]_[lpf/hpf]_[x/y/z/mag]` (32 features): low and hight pass filters to supress data emissinos. `lpf` contains more smoothly data and `hpf` highlights data outliers\n\n- `rot_angle` & `rot_angle_vel` & `angular_distance` (3 features): rotation from quanterion, its jerk and distance\n\n***total:***\n7 (base) + 2 + 4 + 4 + 16 + 16 + 32 + 3 = 84 features\n\n","metadata":{}},{"cell_type":"markdown","source":"# Preparations","metadata":{}},{"cell_type":"markdown","source":"### Imports","metadata":{}},{"cell_type":"code","source":"import os, numpy as np, pandas as pd\nfrom pathlib import Path\nimport warnings \nwarnings.filterwarnings(\"ignore\")\n\nimport torch\nimport torch.nn.functional as F\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:07:24.309417Z","iopub.execute_input":"2025-08-03T11:07:24.309712Z","iopub.status.idle":"2025-08-03T11:07:24.314251Z","shell.execute_reply.started":"2025-08-03T11:07:24.309689Z","shell.execute_reply":"2025-08-03T11:07:24.313307Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Constants","metadata":{}},{"cell_type":"code","source":"# Configuration\n\nPROCESS_DATASET = __name__ == '__main__'\n\nif PROCESS_DATASET:\n\n    #RAW_DIR = Path(\"data\")\n    RAW_DIR = Path(\"/kaggle/input/bfrb-dataset\")\n    #EXPORT_DIR = RAW_DIR\n    EXPORT_DIR = Path(\"/kaggle/working\")\n    FILENAME = 'train.csv'\n\n    assert RAW_DIR.exists(), f\"Directory not found: {RAW_DIR}\"\n    if not os.path.exists(EXPORT_DIR):\n        os.mkdir(EXPORT_DIR)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:14:20.869391Z","iopub.execute_input":"2025-08-03T11:14:20.870572Z","iopub.status.idle":"2025-08-03T11:14:20.878282Z","shell.execute_reply.started":"2025-08-03T11:14:20.870535Z","shell.execute_reply":"2025-08-03T11:14:20.877362Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# FeatureExtractor defenition","metadata":{}},{"cell_type":"markdown","source":"### helpful functions defenitions","metadata":{}},{"cell_type":"code","source":"# Remove Gravity\n\nfrom scipy.spatial.transform import Rotation as R\n\ndef remove_gravity_from_acc(acc_data, rot_data):\n\n    if isinstance(acc_data, pd.DataFrame):\n        acc_values = acc_data[['acc_x', 'acc_y', 'acc_z']].values\n    else:\n        acc_values = acc_data\n\n    if isinstance(rot_data, pd.DataFrame):\n        quat_values = rot_data[['rot_x', 'rot_y', 'rot_z', 'rot_w']].values\n    else:\n        quat_values = rot_data\n\n    num_samples = acc_values.shape[0]\n    linear_accel = np.zeros_like(acc_values)\n    \n    gravity_world = np.array([0, 0, 9.81])\n\n    for i in range(num_samples):\n        if np.all(np.isnan(quat_values[i])) or np.all(np.isclose(quat_values[i], 0)):\n            linear_accel[i, :] = acc_values[i, :] \n            continue\n\n        try:\n            rotation = R.from_quat(quat_values[i])\n            gravity_sensor_frame = rotation.apply(gravity_world, inverse=True)\n            linear_accel[i, :] = acc_values[i, :] - gravity_sensor_frame\n        except ValueError:\n             linear_accel[i, :] = acc_values[i, :]\n             \n    return linear_accel\n\ndef calculate_angular_velocity_from_quat(rot_data, time_delta=1/200): # Assuming 200Hz sampling rate\n    if isinstance(rot_data, pd.DataFrame):\n        quat_values = rot_data[['rot_x', 'rot_y', 'rot_z', 'rot_w']].values\n    else:\n        quat_values = rot_data\n\n    num_samples = quat_values.shape[0]\n    angular_vel = np.zeros((num_samples, 3))\n\n    for i in range(num_samples - 1):\n        q_t = quat_values[i]\n        q_t_plus_dt = quat_values[i+1]\n\n        if np.all(np.isnan(q_t)) or np.all(np.isclose(q_t, 0)) or \\\n           np.all(np.isnan(q_t_plus_dt)) or np.all(np.isclose(q_t_plus_dt, 0)):\n            continue\n\n        try:\n            rot_t = R.from_quat(q_t)\n            rot_t_plus_dt = R.from_quat(q_t_plus_dt)\n\n            # Calculate the relative rotation\n            delta_rot = rot_t.inv() * rot_t_plus_dt\n            \n            # Convert delta rotation to angular velocity vector\n            # The rotation vector (Euler axis * angle) scaled by 1/dt\n            # is a good approximation for small delta_rot\n            angular_vel[i, :] = delta_rot.as_rotvec() / time_delta\n        except ValueError:\n            # If quaternion is invalid, angular velocity remains zero\n            pass\n            \n    return angular_vel\n\ndef calculate_angular_distance(rot_data):\n    if isinstance(rot_data, pd.DataFrame):\n        quat_values = rot_data[['rot_x', 'rot_y', 'rot_z', 'rot_w']].values\n    else:\n        quat_values = rot_data\n\n    num_samples = quat_values.shape[0]\n    angular_dist = np.zeros(num_samples)\n\n    for i in range(num_samples - 1):\n        q1 = quat_values[i]\n        q2 = quat_values[i+1]\n\n        if np.all(np.isnan(q1)) or np.all(np.isclose(q1, 0)) or \\\n           np.all(np.isnan(q2)) or np.all(np.isclose(q2, 0)):\n            angular_dist[i] = 0 # Или np.nan, в зависимости от желаемого поведения\n            continue\n        try:\n            # Преобразование кватернионов в объекты Rotation\n            r1 = R.from_quat(q1)\n            r2 = R.from_quat(q2)\n\n            # Вычисление углового расстояния: 2 * arccos(|real(p * q*)|)\n            # где p* - сопряженный кватернион q\n            # В scipy.spatial.transform.Rotation, r1.inv() * r2 дает относительное вращение.\n            # Угол этого относительного вращения - это и есть угловое расстояние.\n            relative_rotation = r1.inv() * r2\n            \n            # Угол rotation vector соответствует угловому расстоянию\n            # Норма rotation vector - это угол в радианах\n            angle = np.linalg.norm(relative_rotation.as_rotvec())\n            angular_dist[i] = angle\n        except ValueError:\n            angular_dist[i] = 0 # В случае недействительных кватернионов\n            pass\n            \n    return angular_dist\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:07:24.334417Z","iopub.execute_input":"2025-08-03T11:07:24.334710Z","iopub.status.idle":"2025-08-03T11:07:24.622264Z","shell.execute_reply.started":"2025-08-03T11:07:24.334682Z","shell.execute_reply":"2025-08-03T11:07:24.621505Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### `feature_engineering` funtion defenition","metadata":{}},{"cell_type":"code","source":"#################################################################################################################\n####################################              FE               ##############################################\n#################################################################################################################\n\ndef feature_engineering(df: pd.DataFrame):\n    \"\"\"\n    creates new features unsing imu cols in input dataframe,\n    other cols will be ignored\n\n    Args:\n        df (pd.DataFrame): df just like train.csv. *Note:* this \\\n        function will change it.\n\n    Returns:\n        df (pd.DataFrame): input dataframe with extra features in \\\n        propriate columns order: `meta_cols + imu_cols + new_cols \\\n        + thm_cols + tof_cols`\n    \"\"\"\n\n    #######################\n    ### Feature list ######\n    #######################\n    meta_labels = {'gesture', 'gesture_int', 'sequence_type', 'behavior', 'orientation',\n                   'row_id', 'subject', 'phase', 'sequence_id', 'sequence_counter'}\n    meta_cols = [c for c in df.columns if c in meta_labels]\n    \n    feature_cols = [c for c in df.columns if c not in meta_cols]\n\n    imu_cols = [c for c in feature_cols if not (c.startswith('thm_') or c.startswith('tof_'))]\n    thm_cols = [c for c in feature_cols if c.startswith('thm_')]\n    tof_cols = [c for c in feature_cols if c.startswith('tof_')]\n    #print(f\"  IMU {len(imu_cols)} |  THM {len(thm_cols)}|  TOF {len(tof_cols)} | total {len(feature_cols)} features\")\n\n\n\n    #######################\n    #### ace & gyro #######\n    #######################\n    # IMU magnitude\n    df['acc_mag'] = np.sqrt(df['acc_x']**2 + df['acc_y']**2 + df['acc_z']**2)\n    df['rot_mag'] = np.sqrt(df['rot_x']**2 + df['rot_y']**2 + df['rot_z']**2)\n\n\n\n    #######################\n    #### linear_ace #######\n    #######################\n    # Remove gravity\n    def get_linear_accel(df):\n        res = remove_gravity_from_acc(\n            df[['acc_x', 'acc_y', 'acc_z']],\n            df[['rot_x', 'rot_y', 'rot_z', 'rot_w']]\n        )\n        res = pd.DataFrame(res, columns=['linear_acc_x', 'linear_acc_y', 'linear_acc_z'], index=df.index)\n        return res\n    \n    linear_accel_df = df.groupby('sequence_id').apply(get_linear_accel, include_groups=False)\n    linear_accel_df = linear_accel_df.droplevel('sequence_id')\n    df = df.join(linear_accel_df)\n    \n    df['linear_acc_mag'] = np.sqrt(df['linear_acc_x']**2 + df['linear_acc_y']**2 + df['linear_acc_z']**2)\n\n\n\n    #######################\n    #### angular_vel ######\n    #######################\n    # Calc angular velocity\n    def calc_angular_velocity(df):\n        res = calculate_angular_velocity_from_quat( df[['rot_x', 'rot_y', 'rot_z', 'rot_w']] )\n        res = pd.DataFrame(res, columns=['angular_vel_x', 'angular_vel_y', 'angular_vel_z'], index=df.index)\n        return res\n    \n    angular_velocity_df = df.groupby('sequence_id').apply(calc_angular_velocity, include_groups=False)\n    angular_velocity_df = angular_velocity_df.droplevel('sequence_id')\n    df = df.join(angular_velocity_df)\n\n    df['angular_vel_mag'] = np.sqrt(df['angular_vel_x']**2 + df['angular_vel_y']**2 + df['angular_vel_z']**2)\n    \n\n\n    #########################\n    ### apply jerk & pow ####\n    #########################\n    grouped_df = df.groupby('sequence_id')\n    for fe in ('acc', 'rot', 'linear_acc', 'angular_vel'):\n        for dir in ('x', 'y', 'z', 'mag'):\n            # IMU jerks [acc/gyro/linear_acc]_jerk_[x/y/z/mag]\n            df[fe + '_jerk_' + dir] = grouped_df[fe + '_' + dir].diff().fillna(0)\n\n            # IMU Energy (=pow) [acc/gyro/linear_acc]_pow_[x/y/z/mag]\n            df[fe + '_pow_' + dir] = df[fe + '_' + dir] ** 2\n\n\n\n    ###################################\n    ### angle & ang_vel & ang_dist ####\n    ###################################\n    # IMU angle\n    df['rot_angle'] = 2 * np.arccos(df['rot_w'].clip(-1, 1))\n\n    # IMU angular velocity\n    df['rot_angle_vel'] = df.groupby('sequence_id')['rot_angle'].diff().fillna(0)\n    \n    # Calculating angular distance\n    def calc_angular_distance(df):\n        res = calculate_angular_distance(df[['rot_x', 'rot_y', 'rot_z', 'rot_w']])\n        res = pd.DataFrame(res, columns=['angular_distance'], index=df.index)\n        return res\n    \n    angular_distance_df = df.groupby('sequence_id').apply(calc_angular_distance, include_groups=False)\n    angular_distance_df = angular_distance_df.droplevel('sequence_id')\n    df = df.join(angular_distance_df)\n\n\n\n    #########################\n    ####### LPF & HPF #######\n    #########################\n    def create_gaussian_kernel(size: int, channels: int):\n        \"\"\"Create gaussian kernel for smoothing\"\"\"\n        kernel = torch.tensor([np.exp(-(i - size//2)**2/2) for i in range(size)], dtype=torch.float32)\n        kernel = kernel / kernel.sum()\n        kernel = kernel.repeat(channels, 1, 1)  # (out_channels, in_channels/groups, kernel_size)\n        return kernel\n    \n    k = 15\n    grouped = df.groupby('sequence_id')\n    \n    for fe in ('acc', 'rot', 'linear_acc', 'angular_vel'):\n        for dir in ('x', 'y', 'z', 'mag'):\n            col_name = f'{fe}_{dir}'\n            weight = create_gaussian_kernel(k, 1)  # 1 channel, cause process 1 column per iteration\n            \n            # apply convolution to each group\n            lpf_results = []\n            for _, group in grouped:\n                # convert to tensor and add dimentions (batch, channel, length)\n                data = torch.tensor(group[col_name].values, dtype=torch.float32).unsqueeze(0).unsqueeze(0)\n                # apply convolution\n                lpf = F.conv1d(data, weight, padding=k//2)\n                lpf_results.append(lpf.squeeze().numpy())\n            \n            # concatenate results\n            lpf_series = pd.concat([pd.Series(x, index=group.index) for x, (_, group) in zip(lpf_results, grouped)])\n            df[f'{fe}_lpf_{dir}'] = lpf_series\n            df[f'{fe}_hpf_{dir}'] = df[col_name] - df[f'{fe}_lpf_{dir}']\n\n\n\n\n    new_cols = [ c for c in df.columns if (c not in meta_cols) and (c not in feature_cols) ]\n    df[imu_cols + new_cols] = df[imu_cols + new_cols].ffill().bfill().fillna(0).values.astype('float32')\n    \n    return df[meta_cols + imu_cols + new_cols + thm_cols + tof_cols]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:07:24.624208Z","iopub.execute_input":"2025-08-03T11:07:24.624637Z","iopub.status.idle":"2025-08-03T11:07:24.646781Z","shell.execute_reply.started":"2025-08-03T11:07:24.624614Z","shell.execute_reply":"2025-08-03T11:07:24.645806Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Process train dataset","metadata":{}},{"cell_type":"markdown","source":"### load dataset","metadata":{}},{"cell_type":"code","source":"if PROCESS_DATASET:\n    \n    print(\"▶ loading dataset …\")\n    df = pd.read_csv(RAW_DIR / FILENAME)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:07:24.647810Z","iopub.execute_input":"2025-08-03T11:07:24.648141Z","iopub.status.idle":"2025-08-03T11:07:59.075307Z","shell.execute_reply.started":"2025-08-03T11:07:24.648114Z","shell.execute_reply":"2025-08-03T11:07:59.074338Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### implement feature engineering","metadata":{}},{"cell_type":"code","source":"if PROCESS_DATASET:\n    df = feature_engineering(df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:07:59.077321Z","iopub.execute_input":"2025-08-03T11:07:59.077618Z","iopub.status.idle":"2025-08-03T11:13:04.769267Z","shell.execute_reply.started":"2025-08-03T11:07:59.077592Z","shell.execute_reply":"2025-08-03T11:13:04.768245Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if PROCESS_DATASET:\n    meta_cols = {'gesture', 'gesture_int', 'sequence_type', 'behavior', 'orientation',\n                    'row_id', 'subject', 'phase', 'sequence_id', 'sequence_counter'}\n    feature_cols = [c for c in df.columns if c not in meta_cols]\n\n    imu_cols = [c for c in feature_cols if not (c.startswith('thm_') or c.startswith('tof_'))]\n    thm_cols = [c for c in feature_cols if c.startswith('thm_')]\n    tof_cols = [c for c in feature_cols if c.startswith('tof_')]\n    print(f\"  IMU {len(imu_cols)} | THM {len(thm_cols)} | TOF {len(tof_cols)} | total {len(feature_cols)} features\")\n\n\n\nif PROCESS_DATASET:\n    import seaborn as sns\n    import matplotlib.pyplot as plt\n\n    # Матрица корреляций\n    corr_matrix = df.loc[:20000, imu_cols].corr()\n    sns.heatmap(corr_matrix, annot=False)\n    plt.show()\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:13:04.770596Z","iopub.execute_input":"2025-08-03T11:13:04.770948Z","iopub.status.idle":"2025-08-03T11:13:06.865201Z","shell.execute_reply.started":"2025-08-03T11:13:04.770919Z","shell.execute_reply":"2025-08-03T11:13:06.863684Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### save new dataset","metadata":{}},{"cell_type":"code","source":"if PROCESS_DATASET:\n    df.to_csv(EXPORT_DIR / f'featured_{FILENAME}', index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T11:14:28.148295Z","iopub.execute_input":"2025-08-03T11:14:28.148753Z","iopub.status.idle":"2025-08-03T11:18:22.985667Z","shell.execute_reply.started":"2025-08-03T11:14:28.148723Z","shell.execute_reply":"2025-08-03T11:18:22.984876Z"}},"outputs":[],"execution_count":null}]}