{"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"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":102335,"databundleVersionId":12518947,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"%%writefile cmi_utils.py\nimport numpy as np\nimport pandas as pd\nfrom scipy.spatial.transform import Rotation as R\nfrom sklearn.metrics import f1_score\n\n\n############ Quaternion Transformation #############\n\ndef quaternion_to_euler(quat):\n    # [N x,y,z,w]\n    rotation = R.from_quat(quat)\n    # 将旋转对象转换为欧拉角（以弧度为单位）\n    euler_angles = rotation.as_euler('zyx', degrees=False) \n    euler_z,euler_y,euler_x = euler_angles[:, 0], euler_angles[:, 1], euler_angles[:, 2]\n    return euler_x, euler_y, euler_z  \n\n\ndef euler_to_quaternion(euler_x, euler_y, euler_z):\n    \"\"\"将欧拉角转换为四元数\n    参数:\n        roll (float): X轴旋转角度（弧度）\n        pitch (float): Y轴旋转角度（弧度）\n        yaw (float): Z轴旋转角度（弧度）\n        order (str): 旋转顺序，如 'zyx', 'xyz' 等\n    返回:\n        np.array: 四元数 [w, x, y, z]\n    \"\"\"\n    # 使用SciPy的Rotation类（默认顺序为 'xyz'，注意角度需为弧度）\n    rotation = R.from_euler('zyx', np.array([euler_z, euler_y, euler_x]).T, degrees=False)\n    q = rotation.as_quat()  # 返回 [x, y, z, w]\n    return np.array([q[:,3],q[:,0],q[:,1],q[:,2]]).T  # 返回 [w, x,y, z]\n\n\ndef quaternion_multiply(q2, q1):\n    '''\n    params:  q1,q2: [N,[w,x,y,z]]\n    return:  [N,[w,x,y,z]]\n    '''\n    if len(q2.shape)==1:\n        w2, x2, y2, z2 = q2[0],q2[1],q2[2],q2[3]\n    else:\n        w2, x2, y2, z2 = q2[:,0],q2[:,1],q2[:,2],q2[:,3]\n\n    if len(q1.shape)==1:\n        w1, x1, y1, z1 = q1[0],q1[1],q1[2],q1[3]\n    else:\n        w1, x1, y1, z1 = q1[:,0],q1[:,1],q1[:,2],q1[:,3]\n    return np.array([\n        w2*w1 - x2*x1 - y2*y1 - z2*z1,\n        w2*x1 + x2*w1 + y2*z1 - z2*y1,\n        w2*y1 - x2*z1 + y2*w1 + z2*x1,\n        w2*z1 + x2*y1 - y2*x1 + z2*w1\n    ]).T\n\n\ndef quaternion_conjugate(q):\n    if len(q.shape)==1:\n        w, x, y, z = q[0],q[1],q[2],q[3]\n    else:\n        w, x, y, z = q[:,0],q[:,1],q[:,2],q[:,3]\n    return np.array([w, -x, -y, -z]).T\n\n\ndef relative_rotation_quaternion(q1, q2, local=True):\n    # q1,q2: [N,[w,x,y,z]]\n    q1_conjugate = quaternion_conjugate(q1)\n    if local:\n        q_rel = quaternion_multiply(q1_conjugate,q2)\n    else:\n        q_rel = quaternion_multiply(q2, q1_conjugate)\n    return q_rel\n\n\ndef flip_tof(x):\n    x = x.reshape((len(x),5,8,8))\n    x[:,[0,1,3]] =  x[:,[0,1,3],:,::-1]\n    x[:,[2,4]] =  x[:,[2,4],::-1]\n    x2 = x[:,2].copy()\n    x[:,2],x[:,4] =  x[:,4],x2\n    return x.reshape((len(x),-1))\n\n############ Data Augmentation #############\n\n\ndef DataAugmentation(df, n_aug, y_grid, z_grid):\n    quat_rot_y = [np.array([np.cos(np.deg2rad(d) / 2), 0, np.sin(np.deg2rad(d) / 2), 0]) \n                  for d in y_grid]\n    quat_rot_z = [np.array([np.cos(np.deg2rad(d)/2),0,0,np.sin(np.deg2rad(d)/2)]) \n                  for d in z_grid]\n    acc_cols = ['acc_x','acc_y','acc_z']\n    cols_to_convert = df.filter(regex='^(acc|rot|tof|thm)').columns\n    df[cols_to_convert] = df[cols_to_convert].astype(np.float32)\n    df_new = pd.concat([df]*len(quat_rot_z), ignore_index=True)\n    q1 = df[['rot_w','rot_x','rot_y','rot_z']].to_numpy(np.float32)\n    for i,(q2_z,q2_y) in enumerate(zip(quat_rot_z,quat_rot_y)):\n        df_tmp = df.copy()\n        quat_new = quaternion_multiply(q2_z,q1)\n        rot = R.from_quat(np.array([q2_y[1],q2_y[2],q2_y[3],q2_y[0]]))\n        acc_new = rot.apply(df_tmp[acc_cols]).astype(np.float32)\n        df_tmp[acc_cols] = acc_new\n        q2_y_reverse = quaternion_conjugate(q2_y)\n        quat_rot = quaternion_multiply(quat_new,q2_y_reverse)\n        df_tmp[['rot_w','rot_x','rot_y','rot_z']] = quat_rot\n        df_tmp['sequence_id'] = df_tmp['sequence_id']+f'__aug{i}'\n        df_new.iloc[i*len(df):(i+1)*len(df)] = df_tmp\n    return df_new\n\n\n############ GAF #############\n\ndef gaf_cosine_similarity(time_series_2d):\n    \"\"\"\n    Compute multivariate GAF using cosine similarity between vectors.\n    \n    Parameters\n    ----------\n    time_series_2d : array-like, shape (n_timesteps, n_features)\n        Single multivariate time series\n        \n    Returns\n    -------\n    gaf_matrix : array, shape (n_timesteps, n_timesteps)\n        GAF matrix with values in [-1, 1]\n    \"\"\"\n    vectors = time_series_2d\n    n_timesteps = len(vectors)\n    \n    # Compute cosine similarity\n    dot_products = np.dot(vectors, vectors.T)\n    norms = np.linalg.norm(vectors, axis=1)\n    norm_matrix = np.outer(norms, norms)\n    \n    # Avoid division by zero\n    mask = norm_matrix > 1e-10\n    gaf_matrix = np.zeros_like(dot_products)\n    gaf_matrix[mask] = np.clip(dot_products[mask] / norm_matrix[mask], -1, 1)\n    \n    return gaf_matrix\n\ndef gaf_angular_sum(time_series_2d, reference_vector=None):\n    \"\"\"\n    Compute multivariate GAF using angular sum approach (true GAF extension).\n    \n    Parameters\n    ----------\n    time_series_2d : array-like, shape (n_timesteps, n_features)\n        Single multivariate time series\n    reference_vector : array-like, optional\n        Reference vector for angle computation. If None, uses [1,0,0,...]\n        \n    Returns\n    -------\n    gaf_matrix : array, shape (n_timesteps, n_timesteps)\n        GAF matrix with values in [-1, 1]\n    \"\"\"\n    vectors = time_series_2d\n    n_timesteps, n_features = vectors.shape\n    \n    # Normalize vectors to unit sphere\n    norms = np.linalg.norm(vectors, axis=1)\n    norms_safe = np.where(norms == 0, 1, norms)\n    unit_vectors = vectors / norms_safe[:, np.newaxis]\n    \n    # Set reference vector\n    if reference_vector is None:\n        reference_vector = np.zeros(n_features)\n        reference_vector[0] = 1\n    else:\n        reference_vector = np.array(reference_vector) / np.linalg.norm(reference_vector)\n    \n    # Compute angles with reference vector\n    angles = []\n    for i in range(len(unit_vectors)):\n        dot_prod = np.dot(unit_vectors[i], reference_vector)\n        dot_prod = np.clip(dot_prod, -1, 1)\n        angle = np.arccos(dot_prod)\n        angles.append(angle)\n    \n    angles = np.array(angles)\n    \n    # Compute cos(θ_i + θ_j) for all pairs (original GAF approach)\n    gaf_matrix = np.cos(np.add.outer(angles, angles))\n    \n    return gaf_matrix\n\ndef gaf_gramian(time_series_2d, normalize_output=True):\n    \"\"\"\n    Compute multivariate GAF using Gramian matrix approach.\n    \n    Parameters\n    ----------\n    time_series_2d : array-like, shape (n_timesteps, n_features)\n        Single multivariate time series\n    normalize_output : bool, optional (default=True)\n        Whether to normalize output to [-1, 1] range\n        \n    Returns\n    -------\n    gaf_matrix : array, shape (n_timesteps, n_timesteps)\n        GAF matrix (values depend on normalize_output parameter)\n    \"\"\"\n    vectors = time_series_2d\n    \n    # Compute Gramian matrix (dot products)\n    gaf_matrix = np.dot(vectors, vectors.T)\n    \n    if normalize_output:\n        # Normalize to [-1, 1] range\n        min_val, max_val = np.min(gaf_matrix), np.max(gaf_matrix)\n        if max_val > min_val:\n            gaf_matrix = 2 * (gaf_matrix - min_val) / (max_val - min_val) - 1\n        else:\n            gaf_matrix = np.zeros_like(gaf_matrix)\n        return np.clip(gaf_matrix, -1, 1)\n    else:\n        return gaf_matrix\n\ndef compute_multivariate_gaf(X_3d, method='cosine', **kwargs):\n    \"\"\"\n    Compute GAF for multivariate time series data.\n    \n    Parameters\n    ----------\n    X_3d : array-like, shape (n_samples, n_timesteps, n_features)\n        Normalized multivariate time series data\n    method : str, optional (default='cosine')\n        GAF computation method:\n        - 'cosine': Cosine similarity between vectors\n        - 'angular_sum': Angular sum approach (true GAF extension)\n        - 'gramian': Gramian matrix approach\n    **kwargs : additional arguments passed to GAF functions\n        \n    Returns\n    -------\n    X_gaf : array, shape (n_samples, n_timesteps, n_timesteps)\n        GAF matrices for each sample\n    \"\"\"\n    \n    n_samples, n_timesteps, n_features = X_3d.shape\n    X_gaf = np.zeros((n_samples, n_timesteps, n_timesteps))\n    \n    # Select GAF function\n    if method == 'cosine':\n        gaf_func = gaf_cosine_similarity\n    elif method == 'angular_sum':\n        gaf_func = lambda x: gaf_angular_sum(x, **kwargs)\n    elif method == 'gramian':\n        gaf_func = lambda x: gaf_gramian(x, **kwargs)\n    else:\n        raise ValueError(f\"Unknown GAF method: {method}\")\n    \n    # Compute GAF for each sample\n    for i in range(n_samples):\n        X_gaf[i] = gaf_func(X_3d[i])\n        \n    return X_gaf\n\n############# fixing outliers #############\n\ndef Fix_outlier(df,imu_only):\n    outlier_idx = df.subject.isin(['SUBJ_019262','SUBJ_045235'])\n    df_outlier = df.loc[outlier_idx]\n    quat = df_outlier[['rot_x','rot_y','rot_z','rot_w']].fillna(\n        {'rot_x':-0.119916,'rot_y':-0.059953,'rot_z':-0.188298,'rot_w':0.360375})\n\n    d = 180\n    quat_rot_z = np.array([np.cos(np.deg2rad(d)/2),0,0,np.sin(np.deg2rad(d)/2)])\n    quat_rot = quaternion_multiply(quat[['rot_w','rot_x','rot_y','rot_z']].values,\n                                   quat_rot_z)\n    quat_rot = quat_rot*((quat_rot[:,0:1]>0)*2-1)\n    df.loc[outlier_idx,['rot_w','rot_x','rot_y','rot_z']]=quat_rot\n    df.loc[outlier_idx,'acc_x'] *= -1\n    df.loc[outlier_idx,'acc_y'] *= -1\n\n    if not imu_only:\n        tof_columns = [f'tof_{i}_v{j}' for i in range(1,6) for j in range(64)]\n        thm_columns = [f'thm_{i}' for i in range(1,6)]\n        tof = df_outlier[tof_columns].values.reshape(-1,5,8,8)\n        tof = np.rot90(tof, k=2, axes=(2, 3))\n        tof[:,[0,1,2,3,4]] = tof[:,[0,3,4,1,2]]\n        tof = tof.reshape(-1,5*64)\n\n        thm = df_outlier[thm_columns].values\n        thm[:,[0,1,2,3,4]] =  thm[:,[0,3,4,1,2]]\n        df.loc[outlier_idx,tof_columns] = tof\n        df.loc[outlier_idx,thm_columns] = thm\n        \n    return df\n\n\n########## miscellaneous ###############\n\n\ndef transform_to_match_distribution(x1, x2):\n    # Sort both arrays\n    x1_sorted = np.sort(x1)\n    x2_sorted = np.sort(x2)\n\n    # Create percentile positions\n    x1_percentiles = np.linspace(0, 100, len(x1_sorted))\n    x2_percentiles = np.linspace(0, 100, len(x2_sorted))\n\n    # Interpolation function to map x2 percentiles to x1 values\n    def transform(values):\n        # Get percentiles for each value in x2\n        # Using 'nearest' interpolation to avoid extrapolation issues\n        percentiles = np.interp(values, x2_sorted, x2_percentiles, \n                                left=0, right=100)\n        # Map these percentiles to x1 values\n        transformed_values = np.interp(percentiles, x1_percentiles, x1_sorted, \n                                     left=x1_sorted[0], right=x1_sorted[-1])\n        return transformed_values\n\n    return transform(x2)\n\n########### Competition metric ###############\n\n\"\"\"\nHierarchical macro F1 metric for the CMI 2025 Challenge.\n\nThis script defines a single entry point `score(solution, submission, row_id_column_name)`\nthat the Kaggle metrics orchestrator will call.\nIt performs validation on submission IDs and computes a combined binary & multiclass F1 score.\n\"\"\"\n\nclass ParticipantVisibleError(Exception):\n    \"\"\"Errors raised here will be shown directly to the competitor.\"\"\"\n    pass\n\n\nclass CompetitionMetric:\n    \"\"\"Hierarchical macro F1 for the CMI 2025 challenge.\"\"\"\n    def __init__(self, data_type = 'str'):\n        self.target_gestures = [\n            'Above ear - pull hair',\n            'Cheek - pinch skin',\n            'Eyebrow - pull hair',\n            'Eyelash - pull hair',\n            'Forehead - pull hairline',\n            'Forehead - scratch',\n            'Neck - pinch skin',\n            'Neck - scratch',\n        ]\n        self.non_target_gestures = [\n            'Write name on leg',\n            'Wave hello',\n            'Glasses on/off',\n            'Text on phone',\n            'Write name in air',\n            'Feel around in tray and pull out an object',\n            'Scratch knee/leg skin',\n            'Pull air toward your face',\n            'Drink from bottle/cup',\n            'Pinch knee/leg skin'\n        ]\n        self.all_classes = self.target_gestures + self.non_target_gestures\n        self.non_target = 'non_target'\n        if data_type != 'str':\n            sorted_gestures = sorted(self.all_classes)\n            self.target_gestures = [sorted_gestures.index(g) for g in self.target_gestures]\n            self.non_target_gestures = [sorted_gestures.index(g) for g in self.non_target_gestures]\n            self.all_classes = np.arange(len(self.all_classes))\n            self.non_target = -1\n        \n\n    def calculate_hierarchical_f1(\n        self,\n        sol: pd.DataFrame,\n        sub: pd.DataFrame\n    ) -> float:\n\n        # Validate gestures\n        invalid_types = {i for i in sub['gesture'].unique() if i not in self.all_classes}\n        if invalid_types:\n            raise ParticipantVisibleError(\n                f\"Invalid gesture values in submission: {invalid_types}\"\n            )\n\n        # Compute binary F1 (Target vs Non-Target)\n        y_true_bin = sol['gesture'].isin(self.target_gestures).values\n        y_pred_bin = sub['gesture'].isin(self.target_gestures).values\n        f1_binary = f1_score(\n            y_true_bin,\n            y_pred_bin,\n            pos_label=True,\n            zero_division=0,\n            average='binary'\n        )\n\n        # Build multi-class labels for gestures\n        y_true_mc = sol['gesture'].apply(lambda x: x if x in self.target_gestures else self.non_target)\n        y_pred_mc = sub['gesture'].apply(lambda x: x if x in self.target_gestures else self.non_target)\n\n        # Compute macro F1 over all gesture classes\n        f1_macro = f1_score(\n            y_true_mc,\n            y_pred_mc,\n            average='macro',\n            zero_division=0\n        )\n        print(f1_binary,f1_macro)\n        return 0.5 * f1_binary + 0.5 * f1_macro\n\n\ndef score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    row_id_column_name: str,\n    data_type: str = 'str'\n) -> float:\n    \"\"\"\n    Compute hierarchical macro F1 for the CMI 2025 challenge.\n\n    Expected input:\n      - solution and submission as pandas.DataFrame\n      - Column 'sequence_id': unique identifier for each sequence\n      - 'gesture': one of the eight target gestures or \"Non-Target\"\n\n    This metric averages:\n    1. Binary F1 on SequenceType (Target vs Non-Target)\n    2. Macro F1 on gesture (mapping non-targets to \"Non-Target\")\n\n    Raises ParticipantVisibleError for invalid submissions,\n    including invalid SequenceType or gesture values.\n\n\n    Examples\n    --------\n    >>> import pandas as pd\n    >>> row_id_column_name = \"id\"\n    >>> solution = pd.DataFrame({'id': range(4), 'gesture': ['Eyebrow - pull hair']*4})\n    >>> submission = pd.DataFrame({'id': range(4), 'gesture': ['Forehead - pull hairline']*4})\n    >>> score(solution, submission, row_id_column_name=row_id_column_name)\n    0.5\n    >>> submission = pd.DataFrame({'id': range(4), 'gesture': ['Text on phone']*4})\n    >>> score(solution, submission, row_id_column_name=row_id_column_name)\n    0.0\n    >>> score(solution, solution, row_id_column_name=row_id_column_name)\n    1.0\n    \"\"\"\n    # Validate required columns\n    for col in (row_id_column_name, 'gesture'):\n        if col not in solution.columns:\n            raise ParticipantVisibleError(f\"Solution file missing required column: '{col}'\")\n        if col not in submission.columns:\n            raise ParticipantVisibleError(f\"Submission file missing required column: '{col}'\")\n\n    metric = CompetitionMetric(data_type)\n    return metric.calculate_hierarchical_f1(solution, submission)\n\n# Competition scoring function based on a specific encoding. \n# Target gestures are indexed from 0 to 7.\ndef score_from_int(true,pred):\n    y0 = np.where(true <= 7,1,0)\n    p0 = np.where(pred <= 7,1,0)\n    y1 = np.where(true <= 7, true, 8)\n    p1 = np.where(pred <= 7, pred, 8)\n    return 0.5 * (f1_score(y0,p0) + f1_score(y1,p1,average='macro'))\n\n########### tf79_utils ##################\n\nimport os, numpy as np, pandas as pd\nKAGGLE = 'kaggle' in os.getcwd()\n\nif KAGGLE:\n    import tensorflow as tf\n    from tensorflow.keras.utils import Sequence, to_categorical, pad_sequences\n    from tensorflow.keras.models import Model, load_model\n    from tensorflow.keras import backend as K\n    from tensorflow.keras.layers import (\n        Input, Conv1D, BatchNormalization, Activation, add, MaxPooling1D, Dropout,\n        Bidirectional, LSTM, GlobalAveragePooling1D, Dense, Multiply, Reshape,\n        Lambda, Concatenate, GRU, GaussianNoise\n    )\n    from sklearn.preprocessing import StandardScaler\n    from scipy.spatial.transform import Rotation as R\n    from pathlib import Path\n    import joblib\n    \n    \n    #Tensor Manipulations\n    def time_sum(x):\n        return K.sum(x, axis=1)\n    \n    def squeeze_last_axis(x):\n        return tf.squeeze(x, axis=-1)\n    \n    def expand_last_axis(x):\n        return tf.expand_dims(x, axis=-1)\n    \n    def se_block(x, reduction=8):\n        ch = x.shape[-1]\n        se = GlobalAveragePooling1D()(x)\n        se = Dense(ch // reduction, activation='relu')(se)\n        se = Dense(ch, activation='sigmoid')(se)\n        se = Reshape((1, ch))(se)\n        return Multiply()([x, se])\n    \n    # Residual CNN Block with SE\n    def residual_se_cnn_block(x, filters, kernel_size, pool_size=2, drop=0.3, wd=1e-4):\n        shortcut = x\n        for _ in range(2):\n            x = Conv1D(filters, kernel_size, padding='same', use_bias=False,\n                       kernel_regularizer=l2(wd))(x)\n            x = BatchNormalization()(x)\n            x = Activation('relu')(x)\n        x = se_block(x)\n        if shortcut.shape[-1] != filters:\n            shortcut = Conv1D(filters, 1, padding='same', use_bias=False,\n                              kernel_regularizer=l2(wd))(shortcut)\n            shortcut = BatchNormalization()(shortcut)\n        x = add([x, shortcut])\n        x = Activation('relu')(x)\n        x = MaxPooling1D(pool_size)(x)\n        x = Dropout(drop)(x)\n        return x\n    \n    def attention_layer(inputs):\n        score = Dense(1, activation='tanh')(inputs)\n        score = Lambda(squeeze_last_axis)(score)\n        weights = Activation('softmax')(score)\n        weights = Lambda(expand_last_axis)(weights)\n        context = Multiply()([inputs, weights])\n        context = Lambda(time_sum)(context)\n        return context\n    \n    \n    # Normalizes and cleans the time series sequence. \n    \n    def preprocess_sequence(df_seq: pd.DataFrame, feature_cols: list[str], \n                            scaler: StandardScaler):\n        mat = df_seq[feature_cols].ffill().bfill().fillna(0).values\n        return scaler.transform(mat).astype('float32')\n    \n    \n    def 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    \n    def 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    \n    \n    def prepare_tf79_input(df_seq: pd.DataFrame):\n        df_seq['acc_mag'] = np.sqrt(df_seq['acc_x']**2 + df_seq['acc_y']**2 + df_seq['acc_z']**2)\n        df_seq['rot_angle'] = 2 * np.arccos(df_seq['rot_w'].clip(-1, 1))\n        df_seq['acc_mag_jerk'] = df_seq['acc_mag'].diff().fillna(0)\n        df_seq['rot_angle_vel'] = df_seq['rot_angle'].diff().fillna(0)\n    \n        acc_cols_for_gravity_removal = ['acc_x', 'acc_y', 'acc_z']\n        rot_cols_for_gravity_removal = ['rot_x', 'rot_y', 'rot_z', 'rot_w']\n    \n        if not all(col in df_seq.columns for col in acc_cols_for_gravity_removal + rot_cols_for_gravity_removal):\n            print(f\"Warning: Missing raw acc/rot columns for gravity removal in predict for sequence. Using raw acc as linear.\")\n            df_seq['linear_acc_x'] = df_seq.get('acc_x', 0)\n            df_seq['linear_acc_y'] = df_seq.get('acc_y', 0)\n            df_seq['linear_acc_z'] = df_seq.get('acc_z', 0)\n        else:\n            acc_data_seq = df_seq[acc_cols_for_gravity_removal]\n            rot_data_seq = df_seq[rot_cols_for_gravity_removal]\n            linear_accel_seq_arr = remove_gravity_from_acc(acc_data_seq, rot_data_seq)\n            \n            df_seq['linear_acc_x'] = linear_accel_seq_arr[:, 0]\n            df_seq['linear_acc_y'] = linear_accel_seq_arr[:, 1]\n            df_seq['linear_acc_z'] = linear_accel_seq_arr[:, 2]\n        \n        df_seq['linear_acc_mag'] = np.sqrt(df_seq['linear_acc_x']**2 + df_seq['linear_acc_y']**2 + df_seq['linear_acc_z']**2)\n        df_seq['linear_acc_mag_jerk'] = df_seq['linear_acc_mag'].diff().fillna(0)\n        \n        # Calculate angular velocity from quaternions in predict function\n        if all(col in df_seq.columns for col in rot_cols_for_gravity_removal):\n            angular_vel_seq_arr = calculate_angular_velocity_from_quat(df_seq[rot_cols_for_gravity_removal])\n            df_seq['angular_vel_x'] = angular_vel_seq_arr[:, 0]\n            df_seq['angular_vel_y'] = angular_vel_seq_arr[:, 1]\n            df_seq['angular_vel_z'] = angular_vel_seq_arr[:, 2]\n    \n            # Calculate angular jerk from angular velocity\n            df_seq['angular_jerk_x'] = df_seq['angular_vel_x'].diff().fillna(0)\n            df_seq['angular_jerk_y'] = df_seq['angular_vel_y'].diff().fillna(0)\n            df_seq['angular_jerk_z'] = df_seq['angular_vel_z'].diff().fillna(0)\n    \n            # Calculate angular snap from angular jerk\n            df_seq['angular_snap_x'] = df_seq['angular_jerk_x'].diff().fillna(0)\n            df_seq['angular_snap_y'] = df_seq['angular_jerk_y'].diff().fillna(0)\n            df_seq['angular_snap_z'] = df_seq['angular_jerk_z'].diff().fillna(0)\n    \n        else:\n            print(f\"Warning: Missing quaternion columns for angular velocity, jerk, and snap calculation in predict. Filling with 0.\")\n            df_seq['angular_vel_x'] = 0\n            df_seq['angular_vel_y'] = 0\n            df_seq['angular_vel_z'] = 0\n            df_seq['angular_jerk_x'] = 0\n            df_seq['angular_jerk_y'] = 0\n            df_seq['angular_jerk_z'] = 0\n            df_seq['angular_snap_x'] = 0\n            df_seq['angular_snap_y'] = 0\n            df_seq['angular_snap_z'] = 0\n    \n        for i in range(1, 6): \n            pixel_cols_tof = [f\"tof_{i}_v{p}\" for p in range(64)]\n            if not all(col in df_seq.columns for col in pixel_cols_tof):\n                print(f\"Warning: Missing some TOF pixel columns for tof_{i} in predict. Filling aggregates with 0.\")\n                df_seq[f'tof_{i}_mean'] = 0\n                df_seq[f'tof_{i}_std']  = 0\n                df_seq[f'tof_{i}_min']  = 0\n                df_seq[f'tof_{i}_max']  = 0\n                continue\n    \n            tof_sensor_data = df_seq[pixel_cols_tof].replace(-1, np.nan)\n            df_seq[f'tof_{i}_mean'] = tof_sensor_data.mean(axis=1)\n            df_seq[f'tof_{i}_std']  = tof_sensor_data.std(axis=1)\n            df_seq[f'tof_{i}_min']  = tof_sensor_data.min(axis=1)\n            df_seq[f'tof_{i}_max']  = tof_sensor_data.max(axis=1)\n            \n        if 'tof_range_across_sensors' in final_feature_cols:\n            tof_mean_cols_for_contrast = [f'tof_{i}_mean' for i in range(1, 6) if f'tof_{i}_mean' in df_seq.columns]\n            thm_cols_for_contrast = [f'thm_{i}' for i in range(1, 6) if f'thm_{i}' in df_seq.columns]\n    \n            if tof_mean_cols_for_contrast:\n                tof_values_for_contrast = df_seq[tof_mean_cols_for_contrast]\n                df_seq['tof_range_across_sensors'] = tof_values_for_contrast.max(axis=1) - tof_values_for_contrast.min(axis=1)\n                df_seq['tof_std_across_sensors'] = tof_values_for_contrast.std(axis=1)\n            else:\n                df_seq['tof_range_across_sensors'] = 0\n                df_seq['tof_std_across_sensors'] = 0\n    \n            if thm_cols_for_contrast:\n                thm_values_for_contrast = df_seq[thm_cols_for_contrast]\n                df_seq['thm_range_across_sensors'] = thm_values_for_contrast.max(axis=1) - thm_values_for_contrast.min(axis=1)\n                df_seq['thm_std_across_sensors'] = thm_values_for_contrast.std(axis=1)\n            else:\n                df_seq['thm_range_across_sensors'] = 0\n                df_seq['thm_std_across_sensors'] = 0\n            \n        df_seq_final_features = pd.DataFrame(index=df_seq.index)\n        for col_name in final_feature_cols:\n            if col_name in df_seq.columns:\n                df_seq_final_features[col_name] = df_seq[col_name]\n            else:\n                print(f\"CRITICAL ERROR IN PREDICT: Feature '{col_name}' expected by model (from final_feature_cols) was NOT generated in df_seq. Filling with 0. THIS IS LIKELY A BUG.\")\n                df_seq_final_features[col_name] = 0 \n                \n        mat_unscaled = df_seq_final_features.ffill().bfill().fillna(0).values.astype('float32')\n        \n        mat_scaled = scaler.transform(mat_unscaled)\n        \n        pad_input = pad_sequences([mat_scaled], maxlen=pad_len, padding='post', truncating='post', dtype='float32')\n        return pad_input\n    \n    \n    custom_objs = {\n        'time_sum': time_sum,\n        'squeeze_last_axis': squeeze_last_axis,\n        'expand_last_axis': expand_last_axis,\n        'se_block': se_block,\n        'residual_se_cnn_block': residual_se_cnn_block,\n        'attention_layer': attention_layer,\n    }\n    \n    PRETRAINED_DIR = Path(\"/kaggle/input/cmi-models-new/tf79-8223\")\n    print(\"tf79_utils: loading artefacts from\", PRETRAINED_DIR)\n    final_feature_cols = np.load(PRETRAINED_DIR / \"feature_cols.npy\", allow_pickle=True).tolist()\n    pad_len        = int(np.load(PRETRAINED_DIR / \"sequence_maxlen.npy\"))\n    scaler         = joblib.load(PRETRAINED_DIR / \"scaler.pkl\")\n    \n    # Re-calculate imu_dim_final based on the actual features that will be used\n    imu_features_in_final_cols = [c for c in final_feature_cols if any(c.startswith(prefix) for prefix in ['linear_acc_', 'acc_', 'rot_', 'angular_vel_', 'angular_jerk_', 'angular_snap_'])]\n    imu_dim_final = len(imu_features_in_final_cols)\n    tof_thm_aggregated_dim_final = len(final_feature_cols) - imu_dim_final\n    \n    \n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null}]}