{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":102335,"databundleVersionId":12518947,"sourceType":"competition"},{"sourceId":242954653,"sourceType":"kernelVersion"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os, json, joblib, numpy as np, pandas as pd\nfrom pathlib import Path\nimport warnings \nwarnings.filterwarnings(\"ignore\")\n\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.utils.class_weight import compute_class_weight\n\nfrom tensorflow.keras.utils import Sequence, to_categorical, pad_sequences\nfrom tensorflow.keras.models import Model, load_model\nfrom 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)\nfrom tensorflow.keras.regularizers import l2\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.callbacks import EarlyStopping\nfrom tensorflow.keras import backend as K\nimport tensorflow as tf\nimport polars as pl\nfrom sklearn.model_selection import StratifiedGroupKFold\nfrom scipy.spatial.transform import Rotation as R\n\nimport gc\nimport librosa\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport plotly.express as px\nimport plotly.graph_objects as go\nimport plotly.io as pio\nimport random\nimport time\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\n\nfrom copy import deepcopy\nfrom glob import glob\nfrom plotly.subplots import make_subplots\nfrom scipy.signal import butter, lfilter\nfrom torch.utils.data import DataLoader, Dataset, Subset\nfrom tqdm import tqdm\n\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\npio.renderers.default = 'iframe'","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-07-08T08:50:50.958993Z","iopub.execute_input":"2025-07-08T08:50:50.959453Z","iopub.status.idle":"2025-07-08T08:51:05.425126Z","shell.execute_reply.started":"2025-07-08T08:50:50.959435Z","shell.execute_reply":"2025-07-08T08:51:05.424070Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import random\ndef seed_everything(seed):\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    random.seed(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    tf.experimental.numpy.random.seed(seed)\n    os.environ['TF_CUDNN_DETERMINISTIC'] = '1'\n    os.environ['TF_DETERMINISTIC_OPS'] = '1'\nseed_everything(seed=42)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.861621Z","iopub.execute_input":"2025-07-08T07:54:53.861876Z","iopub.status.idle":"2025-07-08T07:54:53.867319Z","shell.execute_reply.started":"2025-07-08T07:54:53.861851Z","shell.execute_reply":"2025-07-08T07:54:53.866461Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# (Competition metric will only be imported when TRAINing)\nTRAIN = True                     # ← set to True when you want to train\nRAW_DIR = Path(\"/kaggle/input/cmi-detect-behavior-with-sensor-data\")\nPRETRAINED_DIR = Path(\"/kaggle/input/lb-0-78-quaternions-tf-bilstm-gru-attention\")  # used when TRAIN=False\nEXPORT_DIR = Path(\"./\")                                    # artefacts will be saved here\nBATCH_SIZE = 64\nPAD_PERCENTILE = 95\nLR_INIT = 5e-4\nWD = 3e-3\nMIXUP_ALPHA = 0.4\nEPOCHS = 160\nPATIENCE = 40\n\n\nprint(\"▶ imports ready · tensorflow\", tf.__version__)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.868939Z","iopub.execute_input":"2025-07-08T07:54:53.869179Z","iopub.status.idle":"2025-07-08T07:54:53.887860Z","shell.execute_reply.started":"2025-07-08T07:54:53.869162Z","shell.execute_reply":"2025-07-08T07:54:53.887267Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#Tensor Manipulations\ndef time_sum(x):\n    return K.sum(x, axis=1)\n\ndef squeeze_last_axis(x):\n    return tf.squeeze(x, axis=-1)\n\ndef expand_last_axis(x):\n    return tf.expand_dims(x, axis=-1)\n\ndef 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\ndef 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\ndef 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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.888534Z","iopub.execute_input":"2025-07-08T07:54:53.888740Z","iopub.status.idle":"2025-07-08T07:54:53.905765Z","shell.execute_reply.started":"2025-07-08T07:54:53.888725Z","shell.execute_reply":"2025-07-08T07:54:53.905087Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Normalizes and cleans the time series sequence. \n\ndef preprocess_sequence(df_seq: pd.DataFrame, feature_cols: list[str], scaler: StandardScaler):\n    mat = df_seq[feature_cols].ffill().bfill().fillna(0).values\n    return scaler.transform(mat).astype('float32')\n\n# MixUp the data argumentation in order to regularize the neural network. \n\nclass MixupGenerator(Sequence):\n    def __init__(self, X, y, batch_size, alpha=0.2):\n        self.X, self.y = X, y\n        self.batch = batch_size\n        self.alpha = alpha\n        self.indices = np.arange(len(X))\n    def __len__(self):\n        return int(np.ceil(len(self.X) / self.batch))\n    def __getitem__(self, i):\n        idx = self.indices[i*self.batch:(i+1)*self.batch]\n        Xb, yb = self.X[idx], self.y[idx]\n        lam = np.random.beta(self.alpha, self.alpha)\n        perm = np.random.permutation(len(Xb))\n        X_mix = lam * Xb + (1-lam) * Xb[perm]\n        y_mix = lam * yb + (1-lam) * yb[perm]\n        return X_mix, y_mix\n    def on_epoch_end(self):\n        np.random.shuffle(self.indices)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.906433Z","iopub.execute_input":"2025-07-08T07:54:53.906641Z","iopub.status.idle":"2025-07-08T07:54:53.925435Z","shell.execute_reply.started":"2025-07-08T07:54:53.906618Z","shell.execute_reply":"2025-07-08T07:54:53.924904Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"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\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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.926431Z","iopub.execute_input":"2025-07-08T07:54:53.926712Z","iopub.status.idle":"2025-07-08T07:54:53.950572Z","shell.execute_reply.started":"2025-07-08T07:54:53.926689Z","shell.execute_reply":"2025-07-08T07:54:53.949866Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def 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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.951326Z","iopub.execute_input":"2025-07-08T07:54:53.951625Z","iopub.status.idle":"2025-07-08T07:54:53.971453Z","shell.execute_reply.started":"2025-07-08T07:54:53.951600Z","shell.execute_reply":"2025-07-08T07:54:53.970880Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def build_two_branch_model(pad_len, imu_dim, tof_dim, n_classes, wd=1e-4):\n    inp = Input(shape=(pad_len, imu_dim+tof_dim))\n    imu = Lambda(lambda t: t[:, :, :imu_dim])(inp)\n    tof = Lambda(lambda t: t[:, :, imu_dim:])(inp)\n\n    # IMU deep branch\n    x1 = residual_se_cnn_block(imu, 64, 3, drop=0.1, wd=wd)\n    x1 = residual_se_cnn_block(x1, 128, 5, drop=0.1, wd=wd)\n\n    # TOF/Thermal lighter branch\n    x2 = Conv1D(64, 3, padding='same', use_bias=False, kernel_regularizer=l2(wd))(tof)\n    x2 = BatchNormalization()(x2); x2 = Activation('relu')(x2)\n    x2 = MaxPooling1D(2)(x2); x2 = Dropout(0.2)(x2)\n    x2 = Conv1D(128, 3, padding='same', use_bias=False, kernel_regularizer=l2(wd))(x2)\n    x2 = BatchNormalization()(x2); x2 = Activation('relu')(x2)\n    x2 = MaxPooling1D(2)(x2); x2 = Dropout(0.2)(x2)\n\n    merged = Concatenate()([x1, x2])\n\n    xa = Bidirectional(LSTM(128, return_sequences=True, kernel_regularizer=l2(wd)))(merged)\n    xb = Bidirectional(GRU(128, return_sequences=True, kernel_regularizer=l2(wd)))(merged)\n    xc = GaussianNoise(0.09)(merged)\n    xc = Dense(16, activation='elu')(xc)\n    \n    x = Concatenate()([xa, xb, xc])\n    x = Dropout(0.4)(x)\n    x = attention_layer(x)\n\n    for units, drop in [(256, 0.5), (128, 0.3)]:\n        x = Dense(units, use_bias=False, kernel_regularizer=l2(wd))(x)\n        x = BatchNormalization()(x); x = Activation('relu')(x)\n        x = Dropout(drop)(x)\n\n    out = Dense(n_classes, activation='softmax', kernel_regularizer=l2(wd))(x)\n    return Model(inp, out)\n\ntmp_model = build_two_branch_model(127,7,325,18)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:53.972144Z","iopub.execute_input":"2025-07-08T07:54:53.972398Z","iopub.status.idle":"2025-07-08T07:54:55.212925Z","shell.execute_reply.started":"2025-07-08T07:54:53.972381Z","shell.execute_reply":"2025-07-08T07:54:55.212153Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class SignalTransform:\n    def __init__(self, always_apply: bool = False, p: float = 0.5):\n        self.always_apply = always_apply\n        self.p = p\n\n    def __call__(self, y: np.ndarray):\n        if self.always_apply:\n            return self.apply(y)\n        else:\n            if np.random.rand() < self.p:\n                return self.apply(y)\n            else:\n                return y\n\n    def apply(self, y: np.ndarray):\n        raise NotImplementedError\n\n\nclass Compose:\n    def __init__(self, transforms: list):\n        self.transforms = transforms\n\n    def __call__(self, y: np.ndarray):\n        for trns in self.transforms:\n            y = trns(y)\n        return y\n\n\nclass OneOf:\n    def __init__(self, transforms: list):\n        self.transforms = transforms\n\n    def __call__(self, y: np.ndarray):\n        n_trns = len(self.transforms)\n        trns_idx = np.random.choice(n_trns)\n        trns = self.transforms[trns_idx]\n        return trns(y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:55.215475Z","iopub.execute_input":"2025-07-08T07:54:55.215869Z","iopub.status.idle":"2025-07-08T07:54:55.221870Z","shell.execute_reply.started":"2025-07-08T07:54:55.215850Z","shell.execute_reply":"2025-07-08T07:54:55.220986Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class GaussianNoise1(SignalTransform):\n    def __init__(\n        self, always_apply: bool = False, \n        p: float = 0.5, max_noise_amplitude: float = 0.20, **kwargs\n    ):\n        super().__init__(always_apply, p)\n        self.noise_amplitude = (0.0, max_noise_amplitude)\n\n    def apply(self, x: np.ndarray, **params):\n        noise_amplitude = np.random.uniform(*self.noise_amplitude)\n        noise = np.random.randn(*x.shape)  # shape (L, N)\n        augmented = (x + noise * noise_amplitude).astype(x.dtype)\n        return augmented\n        \n# ✅ Hàm tạo noise theo luật 1/f^β (β = 1 là pink noise)\ndef powerlaw_psd_gaussian(exponent, size):\n    \"\"\"Sinh tín hiệu nhiễu Gaussian với phổ công suất dạng 1/f^exponent\"\"\"\n    rng = np.random.default_rng()\n    freqs = np.fft.rfftfreq(size)\n    freqs[0] = freqs[1]  # tránh chia cho 0\n    amplitudes = 1 / np.power(freqs, exponent / 2.0)\n    phases = rng.uniform(0, 2 * np.pi, len(freqs))\n    spectrum = amplitudes * np.exp(1j * phases)\n    noise = np.fft.irfft(spectrum, n=size)\n    return noise\n\n# ✅ Phiên bản không dùng colorednoise\nclass PinkNoiseSNR(SignalTransform):\n    def __init__(self, always_apply=False, p=0.5, min_snr=5.0, max_snr=20.0, **kwargs):\n        super().__init__(always_apply, p)\n        self.min_snr = min_snr\n        self.max_snr = max_snr\n\n    def apply(self, y: np.ndarray, **params):\n        snr = np.random.uniform(self.min_snr, self.max_snr)\n        a_signal = np.sqrt((y ** 2).max(axis=0))  # độ lớn tín hiệu, từng kênh\n        a_noise = a_signal / (10 ** (snr / 20))   # biên độ nhiễu dựa trên SNR\n\n        # Tạo pink noise cho từng kênh\n        pink_noise = np.stack([\n            powerlaw_psd_gaussian(1, len(y)) for _ in range(y.shape[1])\n        ], axis=1)\n\n        a_pink = np.sqrt((pink_noise ** 2).max(axis=0))  # chuẩn hóa\n        pink_noise_normalized = pink_noise * (a_noise / a_pink)\n\n        augmented = (y + pink_noise_normalized).astype(y.dtype)\n        return augmented\n        \nclass TimeStretch(SignalTransform):\n    def __init__(self, max_rate=1.5, min_rate=0.5, always_apply=False, p=0.5):\n        super().__init__(always_apply, p)\n        self.max_rate = max_rate\n        self.min_rate = min_rate\n        self.always_apply = always_apply\n        self.p = p\n\n    def apply(self, x: np.ndarray):\n        \"\"\"\n        Stretch a 1D or 2D array in time using linear interpolation.\n        - x: np.ndarray of shape (L,) or (L, N)\n        - rate: float, e.g., 1.2 for 20% longer, 0.8 for 20% shorter\n        \"\"\"\n        rate = np.random.uniform(self.min_rate, self.max_rate)\n        L = x.shape[0]\n        L_new = int(L / rate)\n        orig_idx = np.linspace(0, L - 1, num=L)\n        new_idx = np.linspace(0, L - 1, num=L_new)\n\n        if x.ndim == 1:\n            return np.interp(new_idx, orig_idx, x)\n        elif x.ndim == 2:\n            return np.stack([\n                np.interp(new_idx, orig_idx, x[:, i]) for i in range(x.shape[1])\n            ], axis=1)\n        else:\n            raise ValueError(\"Only 1D or 2D arrays are supported.\")\n\nclass TimeShift(SignalTransform):\n    def __init__(self, always_apply=False, p=0.5, max_shift_pct=0.25, padding_mode=\"replace\"):\n        super().__init__(always_apply, p)\n        \n        assert 0 <= max_shift_pct <= 1.0, \"`max_shift_pct` must be between 0 and 1\"\n        assert padding_mode in [\"replace\", \"zero\"], \"`padding_mode` must be either 'replace' or 'zero'\"\n        \n        self.max_shift_pct = max_shift_pct\n        self.padding_mode = padding_mode\n\n    def apply(self, x: np.ndarray, **params):\n        assert x.ndim == 2, \"`x` must be a 2D array with shape (L, N)\"\n        \n        L = x.shape[0]\n        max_shift = int(L * self.max_shift_pct)\n        shift = np.random.randint(-max_shift, max_shift + 1)\n\n        # Roll along time axis (axis=0)\n        augmented = np.roll(x, shift, axis=0)\n\n        if self.padding_mode == \"zero\":\n            if shift > 0:\n                augmented[:shift, :] = 0\n            elif shift < 0:\n                augmented[shift:, :] = 0\n\n        return augmented\n\n\nclass ButterFilter(SignalTransform):\n    def __init__(self, always_apply=False, p=0.5, cutoff_freq=20, sampling_rate=200, order=4):\n        super().__init__(always_apply, p)\n        \n        self.cutoff_freq = cutoff_freq\n        self.sampling_rate = sampling_rate\n        self.order = order\n\n    def apply(self, x: np.ndarray, **params):\n        assert x.ndim == 2, \"`x` must be a 2D array with shape (L, N)\"\n        return self.butter_lowpass_filter(x)\n\n    def butter_lowpass_filter(self, data: np.ndarray):\n        nyquist = 0.5 * self.sampling_rate\n        normal_cutoff = self.cutoff_freq / nyquist\n        b, a = butter(self.order, normal_cutoff, btype='low', analog=False)\n        filtered_data = lfilter(b, a, data, axis=0)  # filter each channel independently\n        return filtered_data\n\nmy_transforms = Compose([\n    OneOf([\n        GaussianNoise1(p=0.3, max_noise_amplitude=0.02),\n        PinkNoiseSNR(p=0.3, min_snr=6.0, max_snr=20.0),\n        ButterFilter(p=0.3, cutoff_freq=15),\n    ]),\n    TimeStretch(p=0.2, min_rate=0.8, max_rate=1.2),\n    TimeShift(p=0.2, max_shift_pct=0.1, padding_mode=\"zero\"),\n])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T07:54:55.222626Z","iopub.execute_input":"2025-07-08T07:54:55.222879Z","iopub.status.idle":"2025-07-08T07:54:55.241721Z","shell.execute_reply.started":"2025-07-08T07:54:55.222857Z","shell.execute_reply":"2025-07-08T07:54:55.241038Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if TRAIN:\n    print(\"▶ TRAIN MODE – loading dataset …\")\n    df = pd.read_csv(RAW_DIR / \"train.csv\")\n\n    train_dem_df = pd.read_csv(RAW_DIR / \"train_demographics.csv\")\n    df_for_groups = pd.merge(df.copy(), train_dem_df, on='subject', how='left')\n\n    le = LabelEncoder()\n    df['gesture_int'] = le.fit_transform(df['gesture'])\n    np.save(EXPORT_DIR / \"gesture_classes.npy\", le.classes_)\n    gesture_classes = le.classes_\n\n    print(\"  Calculating base engineered IMU features (magnitude, angle)...\")\n    df['acc_mag'] = np.sqrt(df['acc_x']**2 + df['acc_y']**2 + df['acc_z']**2)\n    df['rot_angle'] = 2 * np.arccos(df['rot_w'].clip(-1, 1))\n    \n    print(\"  Calculating engineered IMU derivatives (jerk, angular velocity) for original acc_mag...\")\n    df['acc_mag_jerk'] = df.groupby('sequence_id')['acc_mag'].diff().fillna(0)\n    df['rot_angle_vel'] = df.groupby('sequence_id')['rot_angle'].diff().fillna(0)\n\n    print(\"  Removing gravity and calculating linear acceleration features...\")\n    \n    linear_accel_list = []\n    for _, group in df.groupby('sequence_id'):\n        acc_data_group = group[['acc_x', 'acc_y', 'acc_z']]\n        rot_data_group = group[['rot_x', 'rot_y', 'rot_z', 'rot_w']]\n        linear_accel_group = remove_gravity_from_acc(acc_data_group, rot_data_group)\n        linear_accel_list.append(pd.DataFrame(linear_accel_group, columns=['linear_acc_x', 'linear_acc_y', 'linear_acc_z'], index=group.index))\n    \n    df_linear_accel = pd.concat(linear_accel_list)\n    df = pd.concat([df, df_linear_accel], axis=1)\n\n    df['linear_acc_mag'] = np.sqrt(df['linear_acc_x']**2 + df['linear_acc_y']**2 + df['linear_acc_z']**2)\n    df['linear_acc_mag_jerk'] = df.groupby('sequence_id')['linear_acc_mag'].diff().fillna(0)\n\n    print(\"  Calculating angular velocity from quaternion derivatives...\")\n    angular_vel_list = []\n    for _, group in df.groupby('sequence_id'):\n        rot_data_group = group[['rot_x', 'rot_y', 'rot_z', 'rot_w']]\n        angular_vel_group = calculate_angular_velocity_from_quat(rot_data_group)\n        angular_vel_list.append(pd.DataFrame(angular_vel_group, columns=['angular_vel_x', 'angular_vel_y', 'angular_vel_z'], index=group.index))\n    \n    df_angular_vel = pd.concat(angular_vel_list)\n    df = pd.concat([df, df_angular_vel], axis=1)\n\n    print(\"  Calculating angular jerk from angular velocity...\")\n    df['angular_jerk_x'] = df.groupby('sequence_id')['angular_vel_x'].diff().fillna(0)\n    df['angular_jerk_y'] = df.groupby('sequence_id')['angular_vel_y'].diff().fillna(0)\n    df['angular_jerk_z'] = df.groupby('sequence_id')['angular_vel_z'].diff().fillna(0)\n\n    print(\"  Calculating angular snap from angular jerk...\")\n    df['angular_snap_x'] = df.groupby('sequence_id')['angular_jerk_x'].diff().fillna(0)\n    df['angular_snap_y'] = df.groupby('sequence_id')['angular_jerk_y'].diff().fillna(0)\n    df['angular_snap_z'] = df.groupby('sequence_id')['angular_jerk_z'].diff().fillna(0)\n\n    meta_cols = { } # This was an empty dict in your provided code, keeping it as is.\n\n    imu_cols_base = ['linear_acc_x', 'linear_acc_y', 'linear_acc_z']\n    imu_cols_base.extend([c for c in df.columns if c.startswith('rot_') and c not in ['rot_angle', 'rot_angle_vel']])\n    \n    imu_engineered_features = [\n        'acc_mag', 'rot_angle',\n        'acc_mag_jerk', 'rot_angle_vel',\n        'linear_acc_mag', 'linear_acc_mag_jerk',\n        'angular_vel_x', 'angular_vel_y', 'angular_vel_z',\n        'angular_jerk_x', 'angular_jerk_y', 'angular_jerk_z',\n        'angular_snap_x', 'angular_snap_y', 'angular_snap_z' # Added new angular snap features\n    ]\n    imu_cols = imu_cols_base + imu_engineered_features\n    imu_cols = list(dict.fromkeys(imu_cols))\n\n    thm_cols_original = [c for c in df.columns if c.startswith('thm_')]\n    \n    tof_aggregated_cols_template = []\n    for i in range(1, 6):\n        tof_aggregated_cols_template.extend([f'tof_{i}_mean', f'tof_{i}_std', f'tof_{i}_min', f'tof_{i}_max'])\n\n    final_feature_cols = imu_cols + thm_cols_original + tof_aggregated_cols_template\n    imu_dim_final = len(imu_cols)\n    tof_thm_aggregated_dim_final = len(thm_cols_original) + len(tof_aggregated_cols_template)\n    \n    print(f\"  IMU (incl. engineered & derivatives) {imu_dim_final} | THM + Aggregated TOF {tof_thm_aggregated_dim_final} | total {len(final_feature_cols)} features\")\n    np.save(EXPORT_DIR / \"feature_cols.npy\", np.array(final_feature_cols))\n\n    print(\"  Building sequences with aggregated TOF and preparing data for scaler...\")\n    seq_gp = df.groupby('sequence_id') \n    \n    all_steps_for_scaler_list = []\n    X_list_unscaled, y_list_int_for_stratify, lens, phase_list = [], [], [], []\n    \n    print(\"  Building sequences with aggregated TOF and preparing data for scaler...\")\n    seq_gp = df.groupby('sequence_id') \n    \n    for seq_id, seq_df_orig in seq_gp:\n        seq_df = seq_df_orig.copy()\n    \n        for i in range(1, 6):\n            pixel_cols_tof = [f\"tof_{i}_v{p}\" for p in range(64)]\n            tof_sensor_data = seq_df[pixel_cols_tof].replace(-1, np.nan)\n            seq_df[f'tof_{i}_mean'] = tof_sensor_data.mean(axis=1)\n            seq_df[f'tof_{i}_std']  = tof_sensor_data.std(axis=1)\n            seq_df[f'tof_{i}_min']  = tof_sensor_data.min(axis=1)\n            seq_df[f'tof_{i}_max']  = tof_sensor_data.max(axis=1)\n    \n        mat_with_phase = seq_df[final_feature_cols + ['phase']].ffill().bfill().fillna(0)\n        phase_column = mat_with_phase['phase'].values\n        mat_unscaled = mat_with_phase.drop(columns=['phase']).values.astype('float32')\n    \n        X_list_unscaled.append(mat_unscaled)\n        y_list_int_for_stratify.append(seq_df['gesture_int'].iloc[0])\n        lens.append(len(mat_unscaled))\n        phase_list.append(phase_column)\n    \n    print(\"  Fitting StandardScaler...\")\n    all_steps_concatenated = np.concatenate(X_list_unscaled, axis=0)\n    scaler = StandardScaler().fit(all_steps_concatenated)\n    joblib.dump(scaler, EXPORT_DIR / \"scaler.pkl\")\n    del all_steps_concatenated\n    \n    print(\"  Scaling sequences...\")\n    X_scaled_list = [scaler.transform(x_seq) for x_seq in X_list_unscaled]\n    del X_list_unscaled\n    \n    pad_len = int(np.percentile(lens, PAD_PERCENTILE))\n    np.save(EXPORT_DIR / \"sequence_maxlen.npy\", pad_len)\n    \n    X = pad_sequences(X_scaled_list, maxlen=pad_len, padding='post', truncating='post', dtype='float32')\n    y_int_for_stratify = np.array(y_list_int_for_stratify)\n    y = to_categorical(y_int_for_stratify, num_classes=len(le.classes_))\n    \n    print(\"  Splitting data and preparing for training...\")\n    X_tr, X_val, y_tr, y_val, phase_tr, _ = train_test_split(\n        X, y, phase_list, test_size=0.2, random_state=31195, stratify=y_int_for_stratify\n    )\n\n    #######################################################################################\n    # DATA AUGMENTATION FOR TRAINING SET\n    #######################################################################################\n    \n    print(f\"  Original training set size: {len(X_tr)} sequences\")\n    print(\"  Applying data augmentation to training set...\")\n    \n    X_augmented = []\n    y_augmented = []\n    augmented_count = 0\n    skipped_count = 0\n    \n    for i in range(len(X_tr)):\n        x_seq = X_tr[i]\n        y_seq = y_tr[i]\n        \n        # Tính độ dài thực tế của sequence (bỏ qua padding zeros)\n        # Tìm hàng cuối cùng không phải toàn zeros\n        non_zero_rows = np.any(x_seq != 0, axis=1)\n        if np.any(non_zero_rows):\n            actual_length = np.max(np.where(non_zero_rows)[0]) + 1\n        else:\n            actual_length = 0\n        \n        # Chỉ augment nếu sequence có độ dài >50\n        if actual_length <= 50:\n            skipped_count += 1\n            continue\n        \n        # Apply augmentation transforms cho toàn bộ sequence\n        try:\n            x_seq_aug = my_transforms(x_seq)\n            \n            # Xử lý trường hợp TimeStretch thay đổi độ dài\n            if x_seq_aug.shape[0] != pad_len:\n                if x_seq_aug.shape[0] < pad_len:\n                    # Pad thêm nếu ngắn hơn\n                    pad_amount = pad_len - x_seq_aug.shape[0]\n                    padding = np.zeros((pad_amount, x_seq_aug.shape[1]), dtype=np.float32)\n                    x_seq_aug = np.vstack([x_seq_aug, padding])\n                else:\n                    # Truncate nếu dài hơn\n                    x_seq_aug = x_seq_aug[:pad_len]\n            \n            X_augmented.append(x_seq_aug)\n            y_augmented.append(y_seq)\n            augmented_count += 1\n        except Exception as e:\n            print(f\"  Warning: Augmentation failed for sequence {i}: {e}\")\n            skipped_count += 1\n    \n    # Kết hợp dữ liệu gốc với dữ liệu augmented\n    if X_augmented:\n        X_tr = np.concatenate([X_tr, np.array(X_augmented)], axis=0)\n        y_tr = np.concatenate([y_tr, np.array(y_augmented)], axis=0)\n        \n        print(f\"  ✅ Successfully augmented {augmented_count} sequences\")\n        print(f\"  ⚠️  Skipped {skipped_count} sequences (≤50 steps)\")\n    else:\n        print(f\"  ⚠️  No sequences were augmented (all had ≤50 steps)\")\n    \n    print(f\"  Final training set size: {len(X_tr)} sequences (after augmentation)\")\n    \n    # Shuffle lại training data để trộn augmented data\n    indices = np.random.permutation(len(X_tr))\n    X_tr = X_tr[indices]\n    y_tr = y_tr[indices]\n    \n    print(\"  Training data shuffled after augmentation\")\n\n    #######################################################################################\n    cw_vals = compute_class_weight('balanced', classes=np.arange(len(le.classes_)), y=y_int_for_stratify)\n    class_weight = dict(enumerate(cw_vals))\n\n    model = build_two_branch_model(pad_len, imu_dim_final, tof_thm_aggregated_dim_final, len(le.classes_), wd=WD)\n    \n    steps = len(X_tr) // BATCH_SIZE\n    lr_sched = tf.keras.optimizers.schedules.CosineDecayRestarts(5e-4, first_decay_steps=15 * steps) \n    \n    model.compile(optimizer=Adam(lr_sched),\n                  loss=tf.keras.losses.CategoricalCrossentropy(label_smoothing=0.1),\n                  metrics=['accuracy'])\n\n    train_gen = MixupGenerator(X_tr, y_tr, batch_size=BATCH_SIZE, alpha=MIXUP_ALPHA)\n    cb = EarlyStopping(patience=PATIENCE, restore_best_weights=True, verbose=1, monitor='val_accuracy', mode='max')\n    \n    print(\"  Starting model training...\")\n    model.fit(train_gen, epochs=EPOCHS, validation_data=(X_val, y_val),\n              class_weight=class_weight, callbacks=[cb], verbose=1)\n\n    model.save(EXPORT_DIR / \"gesture_two_branch_mixup.h5\")\n    print(\"✔ Training done – artefacts saved in\", EXPORT_DIR)\n\n    from cmi_2025_metric_copy_for_import import CompetitionMetric\n    preds_val = model.predict(X_val).argmax(1)\n    true_val_int  = y_val.argmax(1)\n    \n    h_f1 = CompetitionMetric().calculate_hierarchical_f1(\n        pd.DataFrame({'gesture': le.classes_[true_val_int]}),\n        pd.DataFrame({'gesture': le.classes_[preds_val]}))\n    print(\"Hold‑out H‑F1 =\", round(h_f1, 4))\nelse:\n    print(\"▶ INFERENCE MODE – 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    gesture_classes = np.load(PRETRAINED_DIR / \"gesture_classes.npy\", allow_pickle=True)\n\n    # Re-calculate imu_dim_final based on the actual features that will be used\n    # Убедитесь, что 'angular_distance' учитывается здесь при инференсе\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_distance'])]\n    imu_dim_final = len(imu_features_in_final_cols)\n\n    tof_thm_aggregated_dim_final = len(final_feature_cols) - imu_dim_final\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    model = load_model(PRETRAINED_DIR / \"gesture_two_branch_mixup.h5\",\n                       compile=False, custom_objects=custom_objs)\n    print(\"  Model, scaler, feature_cols, pad_len loaded – ready for evaluation\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T08:00:54.643122Z","iopub.execute_input":"2025-07-08T08:00:54.643832Z","iopub.status.idle":"2025-07-08T08:07:08.569881Z","shell.execute_reply.started":"2025-07-08T08:00:54.643804Z","shell.execute_reply":"2025-07-08T08:07:08.569275Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"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\nmodels = []\n\nmodel = load_model(\"/kaggle/working/gesture_two_branch_mixup.h5\",\n                   compile=False, custom_objects=custom_objs)\nmodels.append(model)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T08:07:08.571339Z","iopub.execute_input":"2025-07-08T08:07:08.571656Z","iopub.status.idle":"2025-07-08T08:07:08.918708Z","shell.execute_reply.started":"2025-07-08T08:07:08.571635Z","shell.execute_reply":"2025-07-08T08:07:08.918114Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def predict(sequence: pl.DataFrame, demographics: pl.DataFrame) -> str:\n    df_seq = sequence.to_pandas()\n\n    # --- Feature engineering ---\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 = ['acc_x', 'acc_y', 'acc_z']\n    rot_cols = ['rot_x', 'rot_y', 'rot_z', 'rot_w']\n\n    if not all(col in df_seq.columns for col in acc_cols + rot_cols):\n        print(\"Warning: Missing acc/rot columns for gravity removal. 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]\n        rot_data_seq = df_seq[rot_cols]\n        linear_accel_seq_arr = remove_gravity_from_acc(acc_data_seq, rot_data_seq)\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    # --- Angular velocity ---\n    if all(col in df_seq.columns for col in rot_cols):\n        angular_vel_seq_arr = calculate_angular_velocity_from_quat(df_seq[rot_cols])\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    else:\n        print(\"Warning: Missing quaternion columns for angular velocity. 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\n    # --- Angular jerk and snap ---\n    for axis in ['x', 'y', 'z']:\n        vel_col = f'angular_vel_{axis}'\n        if vel_col in df_seq.columns:\n            df_seq[f'angular_jerk_{axis}'] = df_seq[vel_col].diff().fillna(0)\n            df_seq[f'angular_snap_{axis}'] = df_seq[f'angular_jerk_{axis}'].diff().fillna(0)\n        else:\n            df_seq[f'angular_jerk_{axis}'] = 0\n            df_seq[f'angular_snap_{axis}'] = 0\n\n    # --- Angular distance ---\n    if all(col in df_seq.columns for col in rot_cols):\n        angular_dist_seq_arr = calculate_angular_distance(df_seq[rot_cols])\n        df_seq['angular_distance'] = angular_dist_seq_arr\n    else:\n        print(\"Warning: Missing quaternion columns for angular distance. Filling with 0.\")\n        df_seq['angular_distance'] = 0\n\n    # --- TOF sensors ---\n    for i in range(1, 6):\n        pixel_cols = [f\"tof_{i}_v{j}\" for j in range(64)]\n        if not all(col in df_seq.columns for col in pixel_cols):\n            print(f\"Warning: Missing TOF pixel columns for tof_{i}. 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_data = df_seq[pixel_cols].replace(-1, np.nan)\n        df_seq[f'tof_{i}_mean'] = tof_data.mean(axis=1)\n        df_seq[f'tof_{i}_std'] = tof_data.std(axis=1)\n        df_seq[f'tof_{i}_min'] = tof_data.min(axis=1)\n        df_seq[f'tof_{i}_max'] = tof_data.max(axis=1)\n\n    # --- Cross-sensor contrasts ---\n    if 'tof_range_across_sensors' in final_feature_cols:\n        tof_mean_cols = [f'tof_{i}_mean' for i in range(1, 6) if f'tof_{i}_mean' in df_seq.columns]\n        thm_cols = [f'thm_{i}' for i in range(1, 6) if f'thm_{i}' in df_seq.columns]\n\n        if tof_mean_cols:\n            tof_vals = df_seq[tof_mean_cols]\n            df_seq['tof_range_across_sensors'] = tof_vals.max(axis=1) - tof_vals.min(axis=1)\n            df_seq['tof_std_across_sensors'] = tof_vals.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:\n            thm_vals = df_seq[thm_cols]\n            df_seq['thm_range_across_sensors'] = thm_vals.max(axis=1) - thm_vals.min(axis=1)\n            df_seq['thm_std_across_sensors'] = thm_vals.std(axis=1)\n        else:\n            df_seq['thm_range_across_sensors'] = 0\n            df_seq['thm_std_across_sensors'] = 0\n\n    # --- Final feature selection ---\n    df_seq_final = pd.DataFrame(index=df_seq.index)\n    for col in final_feature_cols:\n        if col in df_seq.columns:\n            df_seq_final[col] = df_seq[col]\n        else:\n            print(f\"CRITICAL ERROR IN PREDICT: Feature '{col}' expected by model (from final_feature_cols) was NOT generated. Filling with 0. THIS IS LIKELY A BUG.\")\n            df_seq_final[col] = 0\n\n    mat_unscaled = df_seq_final.ffill().bfill().fillna(0).values.astype('float32')\n    mat_scaled = scaler.transform(mat_unscaled)\n    pad_input = pad_sequences([mat_scaled], maxlen=pad_len, padding='post', truncating='post', dtype='float32')\n\n    # --- Model prediction ---\n    predictions = []\n    for model in models:\n        idx = int(model.predict(pad_input, verbose=0).argmax(1)[0])\n        predictions.append(idx)\n\n    final_idx = max(set(predictions), key=predictions.count)\n    return str(gesture_classes[final_idx])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T08:09:59.498710Z","iopub.execute_input":"2025-07-08T08:09:59.499343Z","iopub.status.idle":"2025-07-08T08:09:59.516135Z","shell.execute_reply.started":"2025-07-08T08:09:59.499313Z","shell.execute_reply":"2025-07-08T08:09:59.515156Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import kaggle_evaluation.cmi_inference_server\ninference_server = kaggle_evaluation.cmi_inference_server.CMIInferenceServer(predict)\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    inference_server.serve()\nelse:\n    inference_server.run_local_gateway(\n        data_paths=(\n            '/kaggle/input/cmi-detect-behavior-with-sensor-data/test.csv',\n            '/kaggle/input/cmi-detect-behavior-with-sensor-data/test_demographics.csv',\n        )\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-08T08:09:59.673173Z","iopub.execute_input":"2025-07-08T08:09:59.673389Z","iopub.status.idle":"2025-07-08T08:10:00.754089Z","shell.execute_reply.started":"2025-07-08T08:09:59.673373Z","shell.execute_reply":"2025-07-08T08:10:00.753325Z"}},"outputs":[],"execution_count":null}]}