{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nfrom scipy.spatial.transform import Rotation as R\nfrom scipy.stats import skew, kurtosis\nfrom scipy.signal import find_peaks\nfrom sklearn.preprocessing import LabelEncoder\nimport joblib\n\n# Configuration\nRAW_DIR = Path(\"/kaggle/input/cmi-detect-behavior-with-sensor-data\")  # Update this to your data directory\nOUTPUT_DIR = Path(\"./features\")\nOUTPUT_DIR.mkdir(exist_ok=True)\n\nprint(\"🔧 BFRB Feature Extraction Pipeline for Tree-based Models\")\n\n#=============================================================================\n# Core Feature Engineering Functions\n#=============================================================================\n\ndef remove_gravity_from_acc(acc_data, rot_data):\n    \"\"\"Remove gravity component from accelerometer data\"\"\"\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    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):\n    \"\"\"Calculate angular velocity from quaternion 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_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            delta_rot = rot_t.inv() * rot_t_plus_dt\n            angular_vel[i, :] = delta_rot.as_rotvec() / time_delta\n        except ValueError:\n            pass\n    \n    return angular_vel\n\ndef calculate_angular_distance(rot_data):\n    \"\"\"Calculate angular distance between consecutive quaternions\"\"\"\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\n            continue\n        try:\n            r1 = R.from_quat(q1)\n            r2 = R.from_quat(q2)\n            relative_rotation = r1.inv() * r2\n            angle = np.linalg.norm(relative_rotation.as_rotvec())\n            angular_dist[i] = angle\n        except ValueError:\n            angular_dist[i] = 0\n    \n    return angular_dist\n\n#=============================================================================\n# Statistical Feature Extraction Functions\n#=============================================================================\n\ndef calculate_statistical_features(series, prefix):\n    \"\"\"Calculate comprehensive statistical features for a time series\"\"\"\n    features = {}\n    \n    # Basic statistics\n    features[f'{prefix}_mean'] = series.mean()\n    features[f'{prefix}_std'] = series.std()\n    features[f'{prefix}_min'] = series.min()\n    features[f'{prefix}_max'] = series.max()\n    features[f'{prefix}_median'] = series.median()\n    features[f'{prefix}_range'] = series.max() - series.min()\n    features[f'{prefix}_iqr'] = series.quantile(0.75) - series.quantile(0.25)\n    \n    # Percentiles\n    for q in [10, 25, 75, 90]:\n        features[f'{prefix}_q{q}'] = series.quantile(q/100)\n    \n    # Shape statistics\n    features[f'{prefix}_skew'] = skew(series.dropna())\n    features[f'{prefix}_kurtosis'] = kurtosis(series.dropna())\n    \n    # Energy and power\n    features[f'{prefix}_energy'] = np.sum(series**2)\n    features[f'{prefix}_power'] = np.mean(series**2)\n    features[f'{prefix}_rms'] = np.sqrt(np.mean(series**2))\n    \n    # Zero crossing rate\n    zero_crossings = np.diff(np.signbit(series.fillna(0))).sum()\n    features[f'{prefix}_zero_crossings'] = zero_crossings\n    \n    # Peak detection\n    try:\n        peaks, _ = find_peaks(series.fillna(0))\n        features[f'{prefix}_num_peaks'] = len(peaks)\n        features[f'{prefix}_peak_density'] = len(peaks) / len(series) if len(series) > 0 else 0\n    except:\n        features[f'{prefix}_num_peaks'] = 0\n        features[f'{prefix}_peak_density'] = 0\n    \n    # Trend features\n    if len(series) > 1:\n        x = np.arange(len(series))\n        slope = np.polyfit(x, series.fillna(method='ffill').fillna(0), 1)[0]\n        features[f'{prefix}_trend_slope'] = slope\n    else:\n        features[f'{prefix}_trend_slope'] = 0\n    \n    # Variability features\n    features[f'{prefix}_cv'] = series.std() / abs(series.mean()) if series.mean() != 0 else 0\n    features[f'{prefix}_mad'] = (series - series.mean()).abs().mean()  # Mean absolute deviation\n    \n    return features\n\ndef calculate_derivative_features(series, prefix):\n    \"\"\"Calculate features based on derivatives (velocity, acceleration)\"\"\"\n    features = {}\n    \n    # First derivative (velocity)\n    velocity = series.diff().fillna(0)\n    features.update(calculate_statistical_features(velocity, f'{prefix}_vel'))\n    \n    # Second derivative (acceleration) \n    acceleration = velocity.diff().fillna(0)\n    features.update(calculate_statistical_features(acceleration, f'{prefix}_accel'))\n    \n    # Jerk (third derivative)\n    jerk = acceleration.diff().fillna(0)\n    features[f'{prefix}_jerk_mean'] = jerk.mean()\n    features[f'{prefix}_jerk_std'] = jerk.std()\n    features[f'{prefix}_jerk_max'] = jerk.max()\n    \n    return features\n\n# def calculate_frequency_features(series, prefix, sampling_rate=200):\n#     \"\"\"Calculate frequency domain features using FFT\"\"\"\n#     features = {}\n    \n#     try:\n#         # Remove DC component and apply window\n#         signal = series.fillna(0) - series.mean()\n        \n#         if len(signal) > 1:\n#             # FFT\n#             fft = np.fft.fft(signal)\n#             freqs = np.fft.fftfreq(len(signal), 1/sampling_rate)\n            \n#             # Power spectral density\n#             psd = np.abs(fft)**2\n            \n#             # Frequency domain statistics\n#             features[f'{prefix}_spectral_centroid'] = np.sum(freqs[:len(freqs)//2] * psd[:len(psd)//2]) / np.sum(psd[:len(psd)//2])\n#             features[f'{prefix}_spectral_rolloff'] = freqs[np.where(np.cumsum(psd[:len(psd)//2]) >= 0.85 * np.sum(psd[:len(psd)//2]))[0][0]]\n#             features[f'{prefix}_spectral_bandwidth'] = np.sqrt(np.sum(((freqs[:len(freqs)//2] - features[f'{prefix}_spectral_centroid'])**2) * psd[:len(psd)//2]) / np.sum(psd[:len(psd)//2]))\n            \n#             # Dominant frequency\n#             dominant_freq_idx = np.argmax(psd[:len(psd)//2])\n#             features[f'{prefix}_dominant_freq'] = freqs[dominant_freq_idx]\n#             features[f'{prefix}_dominant_freq_power'] = psd[dominant_freq_idx]\n            \n#             # Energy in frequency bands\n#             low_freq_mask = (freqs >= 0) & (freqs <= 5)\n#             mid_freq_mask = (freqs > 5) & (freqs <= 15)\n#             high_freq_mask = (freqs > 15) & (freqs <= 50)\n            \n#             total_energy = np.sum(psd[:len(psd)//2])\n#             features[f'{prefix}_low_freq_energy'] = np.sum(psd[low_freq_mask]) / total_energy if total_energy > 0 else 0\n#             features[f'{prefix}_mid_freq_energy'] = np.sum(psd[mid_freq_mask]) / total_energy if total_energy > 0 else 0\n#             features[f'{prefix}_high_freq_energy'] = np.sum(psd[high_freq_mask]) / total_energy if total_energy > 0 else 0\n        \n#         else:\n#             # Default values for short sequences\n#             for feature_name in ['spectral_centroid', 'spectral_rolloff', 'spectral_bandwidth', \n#                                'dominant_freq', 'dominant_freq_power', 'low_freq_energy', \n#                                'mid_freq_energy', 'high_freq_energy']:\n#                 features[f'{prefix}_{feature_name}'] = 0\n                \n#     except Exception as e:\n#         print(f\"Frequency feature calculation failed for {prefix}: {e}\")\n#         for feature_name in ['spectral_centroid', 'spectral_rolloff', 'spectral_bandwidth', \n#                            'dominant_freq', 'dominant_freq_power', 'low_freq_energy', \n#                            'mid_freq_energy', 'high_freq_energy']:\n#             features[f'{prefix}_{feature_name}'] = 0\n    \n#     return features\n\n\ndef calculate_frequency_features(series, prefix, sampling_rate=200):\n    \"\"\"Calculate frequency domain features using FFT\"\"\"\n    features = {}\n    \n    try:\n        # Remove NaN values first - CRITICAL FIX\n        series_clean = series.dropna()\n        if len(series_clean) < 2:  # Need at least 2 points for FFT\n            raise ValueError(\"Not enough data points after removing NaN\")\n        \n        # Remove DC component and apply window\n        signal = series_clean - series_clean.mean()\n        \n        # FFT\n        fft = np.fft.fft(signal)\n        freqs = np.fft.fftfreq(len(signal), 1/sampling_rate)\n        \n        # Power spectral density\n        psd = np.abs(fft)**2\n        \n        # Only use positive frequencies (first half)\n        positive_freq_mask = freqs >= 0\n        positive_freqs = freqs[positive_freq_mask]\n        positive_psd = psd[positive_freq_mask]\n        \n        if len(positive_freqs) == 0 or np.sum(positive_psd) == 0:\n            raise ValueError(\"No valid frequency components\")\n        \n        # Frequency domain statistics\n        features[f'{prefix}_spectral_centroid'] = np.sum(positive_freqs * positive_psd) / np.sum(positive_psd)\n        \n        # Spectral rolloff (85%)\n        cumulative_energy = np.cumsum(positive_psd)\n        rolloff_index = np.where(cumulative_energy >= 0.85 * np.sum(positive_psd))[0]\n        features[f'{prefix}_spectral_rolloff'] = positive_freqs[rolloff_index[0]] if len(rolloff_index) > 0 else 0\n        \n        # Spectral bandwidth\n        features[f'{prefix}_spectral_bandwidth'] = np.sqrt(\n            np.sum(((positive_freqs - features[f'{prefix}_spectral_centroid'])**2) * positive_psd) / np.sum(positive_psd)\n        )\n        \n        # Dominant frequency\n        dominant_freq_idx = np.argmax(positive_psd)\n        features[f'{prefix}_dominant_freq'] = positive_freqs[dominant_freq_idx]\n        features[f'{prefix}_dominant_freq_power'] = positive_psd[dominant_freq_idx]\n        \n        # Energy in frequency bands\n        low_freq_mask = (positive_freqs >= 0) & (positive_freqs <= 5)\n        mid_freq_mask = (positive_freqs > 5) & (positive_freqs <= 15)\n        high_freq_mask = (positive_freqs > 15) & (positive_freqs <= 50)\n        \n        total_energy = np.sum(positive_psd)\n        features[f'{prefix}_low_freq_energy'] = np.sum(positive_psd[low_freq_mask]) / total_energy if total_energy > 0 else 0\n        features[f'{prefix}_mid_freq_energy'] = np.sum(positive_psd[mid_freq_mask]) / total_energy if total_energy > 0 else 0\n        features[f'{prefix}_high_freq_energy'] = np.sum(positive_psd[high_freq_mask]) / total_energy if total_energy > 0 else 0\n        \n    except Exception as e:\n        # Set all frequency features to 0 when calculation fails\n        for feature_name in ['spectral_centroid', 'spectral_rolloff', 'spectral_bandwidth', \n                           'dominant_freq', 'dominant_freq_power', 'low_freq_energy', \n                           'mid_freq_energy', 'high_freq_energy']:\n            features[f'{prefix}_{feature_name}'] = 0\n    \n    return features\n\n\ndef calculate_cross_signal_features(signal1, signal2, prefix1, prefix2):\n    \"\"\"Calculate features between two signals (correlations, etc.)\"\"\"\n    features = {}\n    \n    # Cross-correlation\n    if len(signal1) > 1 and len(signal2) > 1:\n        correlation = np.corrcoef(signal1.fillna(0), signal2.fillna(0))[0, 1]\n        features[f'{prefix1}_{prefix2}_correlation'] = correlation if not np.isnan(correlation) else 0\n        \n        # Cross-correlation lag features\n        cross_corr = np.correlate(signal1.fillna(0), signal2.fillna(0), mode='full')\n        max_corr_idx = np.argmax(np.abs(cross_corr))\n        lag = max_corr_idx - len(signal1) + 1\n        features[f'{prefix1}_{prefix2}_max_cross_corr'] = cross_corr[max_corr_idx]\n        features[f'{prefix1}_{prefix2}_lag'] = lag\n        \n        # Coherence (simplified as correlation of power)\n        power1 = signal1**2\n        power2 = signal2**2\n        power_corr = np.corrcoef(power1.fillna(0), power2.fillna(0))[0, 1]\n        features[f'{prefix1}_{prefix2}_coherence'] = power_corr if not np.isnan(power_corr) else 0\n    else:\n        features[f'{prefix1}_{prefix2}_correlation'] = 0\n        features[f'{prefix1}_{prefix2}_max_cross_corr'] = 0\n        features[f'{prefix1}_{prefix2}_lag'] = 0\n        features[f'{prefix1}_{prefix2}_coherence'] = 0\n    \n    return features\n\n#=============================================================================\n# Comprehensive Feature Extraction for Single Sequence\n#=============================================================================\n\ndef extract_sequence_features(df_seq):\n    \"\"\"Extract comprehensive features from a single sequence\"\"\"\n    features = {}\n    \n    # Basic sequence metadata\n    features['sequence_length'] = len(df_seq)\n    features['sequence_duration'] = features['sequence_length'] / 200  # Assuming 200Hz\n    \n    # Apply core feature engineering first\n    df_seq = df_seq.copy()\n    \n    # Remove gravity from accelerometer data\n    linear_accel = remove_gravity_from_acc(df_seq, df_seq)\n    df_seq['linear_acc_x'] = linear_accel[:, 0]\n    df_seq['linear_acc_y'] = linear_accel[:, 1] \n    df_seq['linear_acc_z'] = linear_accel[:, 2]\n    df_seq['linear_acc_mag'] = np.sqrt(df_seq['linear_acc_x']**2 + \n                                      df_seq['linear_acc_y']**2 + \n                                      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\n    angular_vel = calculate_angular_velocity_from_quat(df_seq)\n    df_seq['angular_vel_x'] = angular_vel[:, 0]\n    df_seq['angular_vel_y'] = angular_vel[:, 1]\n    df_seq['angular_vel_z'] = angular_vel[:, 2]\n    df_seq['angular_vel_mag'] = np.sqrt(angular_vel[:, 0]**2 + angular_vel[:, 1]**2 + angular_vel[:, 2]**2)\n    df_seq['angular_distance'] = calculate_angular_distance(df_seq)\n    \n    # Original accelerometer magnitude\n    df_seq['acc_mag'] = np.sqrt(df_seq['acc_x']**2 + df_seq['acc_y']**2 + df_seq['acc_z']**2)\n    \n    # ===== IMU FEATURES =====\n    imu_signals = {\n        'acc_x': df_seq['acc_x'],\n        'acc_y': df_seq['acc_y'], \n        'acc_z': df_seq['acc_z'],\n        'acc_mag': df_seq['acc_mag'],\n        'linear_acc_x': df_seq['linear_acc_x'],\n        'linear_acc_y': df_seq['linear_acc_y'],\n        'linear_acc_z': df_seq['linear_acc_z'],\n        'linear_acc_mag': df_seq['linear_acc_mag'],\n        'angular_vel_x': df_seq['angular_vel_x'],\n        'angular_vel_y': df_seq['angular_vel_y'],\n        'angular_vel_z': df_seq['angular_vel_z'],\n        'angular_vel_mag': df_seq['angular_vel_mag'],\n        'angular_distance': df_seq['angular_distance'],\n    }\n    \n    # Extract statistical, derivative, and frequency features for each IMU signal\n    for signal_name, signal_data in imu_signals.items():\n        features.update(calculate_statistical_features(signal_data, signal_name))\n        features.update(calculate_derivative_features(signal_data, signal_name))\n        features.update(calculate_frequency_features(signal_data, signal_name))\n    \n    # ===== THERMAL FEATURES =====\n    for i in range(1, 6):\n        thm_col = f'thm_{i}'\n        if thm_col in df_seq.columns:\n            features.update(calculate_statistical_features(df_seq[thm_col], thm_col))\n            features.update(calculate_derivative_features(df_seq[thm_col], thm_col))\n            features.update(calculate_frequency_features(df_seq[thm_col], thm_col))\n    \n    # ===== TIME-OF-FLIGHT FEATURES =====\n    for i in range(1, 6):\n        # Individual pixel features (aggregated)\n        pixel_cols = [f\"tof_{i}_v{p}\" for p in range(64)]\n        if all(col in df_seq.columns for col in pixel_cols):\n            tof_data = df_seq[pixel_cols].replace(-1, np.nan)\n            \n            # Spatial aggregations over pixels at each timestep\n            tof_mean_series = tof_data.mean(axis=1)\n            tof_std_series = tof_data.std(axis=1)\n            tof_min_series = tof_data.min(axis=1)\n            tof_max_series = tof_data.max(axis=1)\n            tof_range_series = tof_max_series - tof_min_series\n            \n            # Temporal features of spatial aggregations\n            for agg_name, agg_series in [\n                ('mean', tof_mean_series), ('std', tof_std_series),\n                ('min', tof_min_series), ('max', tof_max_series), ('range', tof_range_series)\n            ]:\n                prefix = f'tof_{i}_{agg_name}'\n                features.update(calculate_statistical_features(agg_series, prefix))\n                features.update(calculate_derivative_features(agg_series, prefix))\n            \n            # Overall ToF sensor statistics\n            all_tof_values = tof_data.values.flatten()\n            all_tof_values = all_tof_values[~np.isnan(all_tof_values)]\n            if len(all_tof_values) > 0:\n                features[f'tof_{i}_overall_mean'] = np.mean(all_tof_values)\n                features[f'tof_{i}_overall_std'] = np.std(all_tof_values)\n                features[f'tof_{i}_overall_min'] = np.min(all_tof_values)\n                features[f'tof_{i}_overall_max'] = np.max(all_tof_values)\n                features[f'tof_{i}_coverage_ratio'] = len(all_tof_values) / (64 * len(df_seq))\n            else:\n                features[f'tof_{i}_overall_mean'] = 0\n                features[f'tof_{i}_overall_std'] = 0\n                features[f'tof_{i}_overall_min'] = 0\n                features[f'tof_{i}_overall_max'] = 0\n                features[f'tof_{i}_coverage_ratio'] = 0\n    \n    # ===== CROSS-SIGNAL FEATURES =====\n    # Accelerometer cross-correlations\n    features.update(calculate_cross_signal_features(df_seq['acc_x'], df_seq['acc_y'], 'acc_x', 'acc_y'))\n    features.update(calculate_cross_signal_features(df_seq['acc_x'], df_seq['acc_z'], 'acc_x', 'acc_z'))\n    features.update(calculate_cross_signal_features(df_seq['acc_y'], df_seq['acc_z'], 'acc_y', 'acc_z'))\n    \n    # Linear acceleration cross-correlations\n    features.update(calculate_cross_signal_features(df_seq['linear_acc_x'], df_seq['linear_acc_y'], 'lin_acc_x', 'lin_acc_y'))\n    features.update(calculate_cross_signal_features(df_seq['linear_acc_x'], df_seq['linear_acc_z'], 'lin_acc_x', 'lin_acc_z'))\n    features.update(calculate_cross_signal_features(df_seq['linear_acc_y'], df_seq['linear_acc_z'], 'lin_acc_y', 'lin_acc_z'))\n    \n    # Angular velocity cross-correlations\n    features.update(calculate_cross_signal_features(df_seq['angular_vel_x'], df_seq['angular_vel_y'], 'ang_vel_x', 'ang_vel_y'))\n    features.update(calculate_cross_signal_features(df_seq['angular_vel_x'], df_seq['angular_vel_z'], 'ang_vel_x', 'ang_vel_z'))\n    features.update(calculate_cross_signal_features(df_seq['angular_vel_y'], df_seq['angular_vel_z'], 'ang_vel_y', 'ang_vel_z'))\n    \n    # ===== BEHAVIOR-BASED FEATURES =====\n    if 'behavior' in df_seq.columns:\n        behavior_counts = df_seq['behavior'].value_counts()\n        for behavior in ['Transition', 'Pause', 'Gesture']:\n            features[f'behavior_{behavior}_count'] = behavior_counts.get(behavior, 0)\n            features[f'behavior_{behavior}_ratio'] = behavior_counts.get(behavior, 0) / len(df_seq)\n        \n        # Behavior transitions\n        behavior_changes = (df_seq['behavior'] != df_seq['behavior'].shift()).sum()\n        features['behavior_transitions'] = behavior_changes\n    \n    return features\n\n#=============================================================================\n# Main Feature Extraction Pipeline\n#=============================================================================\n\ndef extract_features_from_dataset(csv_path, output_path, demographics_path=None):\n    \"\"\"Extract features from entire dataset and save as DataFrame\"\"\"\n    print(f\"Loading data from {csv_path}\")\n    \n    # Load data in chunks to manage memory\n    chunk_size = 10000\n    all_features = []\n    \n    # Get unique sequence IDs first\n    sequence_ids = pd.read_csv(csv_path, usecols=['sequence_id'])['sequence_id'].unique()\n    total_sequences = len(sequence_ids)\n    \n    print(f\"Found {total_sequences} unique sequences\")\n    \n    # Process sequences in batches\n    batch_size = 100\n    for batch_start in range(0, total_sequences, batch_size):\n        batch_end = min(batch_start + batch_size, total_sequences)\n        batch_seq_ids = sequence_ids[batch_start:batch_end]\n        \n        print(f\"Processing batch {batch_start//batch_size + 1}/{(total_sequences-1)//batch_size + 1} \"\n              f\"(sequences {batch_start+1}-{batch_end})\")\n        \n        # Load data for this batch\n        batch_data = pd.read_csv(csv_path)\n        batch_data = batch_data[batch_data['sequence_id'].isin(batch_seq_ids)]\n        \n        batch_features = []\n        for seq_id in batch_seq_ids:\n            try:\n                seq_data = batch_data[batch_data['sequence_id'] == seq_id].copy()\n                \n                if len(seq_data) == 0:\n                    continue\n                \n                # Extract features\n                seq_features = extract_sequence_features(seq_data)\n                seq_features['sequence_id'] = seq_id\n                \n                # Add metadata\n                if 'subject' in seq_data.columns:\n                    seq_features['subject'] = seq_data['subject'].iloc[0]\n                if 'gesture' in seq_data.columns:\n                    seq_features['gesture'] = seq_data['gesture'].iloc[0]\n                if 'orientation' in seq_data.columns:\n                    seq_features['orientation'] = seq_data['orientation'].iloc[0]\n                if 'sequence_type' in seq_data.columns:\n                    seq_features['sequence_type'] = seq_data['sequence_type'].iloc[0]\n                \n                batch_features.append(seq_features)\n                \n            except Exception as e:\n                print(f\"Error processing sequence {seq_id}: {e}\")\n                continue\n        \n        all_features.extend(batch_features)\n        \n        # Free memory\n        del batch_data\n    \n    # Convert to DataFrame\n    feature_df = pd.DataFrame(all_features)\n    \n    # Add demographic features if available\n    if demographics_path and Path(demographics_path).exists():\n        print(\"Adding demographic features...\")\n        demographics = pd.read_csv(demographics_path)\n        feature_df = feature_df.merge(demographics, on='subject', how='left')\n    \n    # Save features\n    print(f\"Saving {len(feature_df)} feature vectors with {len(feature_df.columns)} features to {output_path}\")\n    feature_df.to_csv(output_path, index=False)\n    \n    # Save feature names for later use\n    feature_names_path = output_path.parent / f\"{output_path.stem}_feature_names.txt\"\n    feature_columns = [col for col in feature_df.columns \n                      if col not in ['sequence_id', 'subject', 'gesture', 'orientation', 'sequence_type']]\n    \n    with open(feature_names_path, 'w') as f:\n        for feature in feature_columns:\n            f.write(f\"{feature}\\n\")\n    \n    print(f\"Feature names saved to {feature_names_path}\")\n    print(f\"Feature extraction completed! Shape: {feature_df.shape}\")\n    \n    return feature_df\n\ndef create_tree_model_ready_dataset():\n    \"\"\"Create train/test datasets ready for tree-based models\"\"\"\n    print(\"🌳 Creating Tree Model Ready Dataset\")\n    \n    # Extract training features\n    train_features = extract_features_from_dataset(\n        csv_path=RAW_DIR / \"train.csv\",\n        output_path=OUTPUT_DIR / \"train_features.csv\",\n        demographics_path=RAW_DIR / \"train_demographics.csv\"\n    )\n    \n    # Extract test features (if available)\n    test_path = RAW_DIR / \"test.csv\"\n    if test_path.exists():\n        test_features = extract_features_from_dataset(\n            csv_path=test_path,\n            output_path=OUTPUT_DIR / \"test_features.csv\",\n            demographics_path=RAW_DIR / \"test_demographics.csv\"\n        )\n    \n    # Create summary report\n    with open(OUTPUT_DIR / \"feature_extraction_report.txt\", 'w') as f:\n        f.write(\"BFRB Feature Extraction Report\\n\")\n        f.write(\"=\"*50 + \"\\n\\n\")\n        \n        f.write(f\"Training Features: {train_features.shape}\\n\")\n        if 'test_features' in locals():\n            f.write(f\"Test Features: {test_features.shape}\\n\")\n        \n        f.write(f\"\\nFeature Categories:\\n\")\n        f.write(\"- Statistical features (mean, std, min, max, percentiles, skew, kurtosis)\\n\")\n        f.write(\"- Derivative features (velocity, acceleration, jerk)\\n\")\n        f.write(\"- Frequency domain features (spectral centroid, rolloff, energy bands)\\n\")\n        f.write(\"- Cross-signal features (correlations, coherence)\\n\")\n        f.write(\"- Engineered features (linear acceleration, angular velocity)\\n\")\n        f.write(\"- Sensor-specific features (ToF spatial aggregations, thermal)\\n\")\n        f.write(\"- Behavioral features (transition counts, phase ratios)\\n\")\n        \n        if 'gesture' in train_features.columns:\n            f.write(f\"\\nGesture Distribution:\\n\")\n            gesture_counts = train_features['gesture'].value_counts()\n            for gesture, count in gesture_counts.items():\n                f.write(f\"  {gesture}: {count}\\n\")\n    \n    print(f\"\\n✅ Feature extraction completed!\")\n    print(f\"📁 Files saved in {OUTPUT_DIR}/\")\n    print(f\"📊 Training features: train_features.csv ({train_features.shape})\")\n    if 'test_features' in locals():\n        print(f\"📊 Test features: test_features.csv ({test_features.shape})\")\n    print(f\"📋 Feature names: train_features_feature_names.txt\")\n    print(f\"📄 Report: feature_extraction_report.txt\")\n    \n    return train_features","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-21T13:34:22.337789Z","iopub.execute_input":"2025-08-21T13:34:22.338654Z","iopub.status.idle":"2025-08-21T13:34:23.557219Z","shell.execute_reply.started":"2025-08-21T13:34:22.338621Z","shell.execute_reply":"2025-08-21T13:34:23.556351Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# Run the feature extraction pipeline\nfeature_df = create_tree_model_ready_dataset()\n\nprint(\"\\n🚀 Ready for tree-based models!\")\nprint(\"You can now use these features with:\")\nprint(\"- XGBoost: xgb.XGBClassifier()\")\nprint(\"- LightGBM: lgb.LGBMClassifier()\")  \nprint(\"- CatBoost: cb.CatBoostClassifier()\")\nprint(\"- Random Forest: sklearn.ensemble.RandomForestClassifier()\")\nprint(\"\\nExample usage:\")\nprint(\"\"\"\nimport pandas as pd\nimport xgboost as xgb\nfrom sklearn.model_selection import train_test_split\n\n# Load features\ndf = pd.read_csv('./features/train_features.csv')\n\n# Prepare data\nfeature_cols = [col for col in df.columns \n               if col not in ['sequence_id', 'subject', 'gesture', 'orientation', 'sequence_type']]\nX = df[feature_cols]\ny = df['gesture']\n\n# Train model\nX_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, stratify=y)\nmodel = xgb.XGBClassifier()\nmodel.fit(X_train, y_train)\n\"\"\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-08-21T13:34:40.675008Z","iopub.execute_input":"2025-08-21T13:34:40.675914Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}