{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81933,"databundleVersionId":9643020,"sourceType":"competition"},{"sourceId":10235049,"sourceType":"datasetVersion","datasetId":6328700}],"dockerImageVersionId":30822,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":72.905702,"end_time":"2024-12-18T10:26:31.543756","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-12-18T10:25:18.638054","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"b002ee23-ed4c-423d-adcd-c58dfcdbcc28","cell_type":"markdown","source":"**PB 0.435, LB 0.460**\n\nI didn't select it for final score for it's mundane PB performance, what a pityXD\n\nUpvote if it's helpful for u","metadata":{}},{"id":"fc79ff03","cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nimport gc\nfrom pathlib import Path\nimport warnings\nwarnings.filterwarnings('ignore')\nfrom scipy import stats, signal\nfrom scipy.fft import fft, fftfreq\nfrom typing import Dict, List, Tuple, Optional, Union\nimport logging\nfrom tqdm import tqdm\nfrom datetime import datetime, time\nimport re\nimport lightgbm as lgb\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import StratifiedKFold, ParameterGrid\nfrom sklearn.metrics import cohen_kappa_score","metadata":{"execution":{"iopub.execute_input":"2024-12-18T10:25:21.304760Z","iopub.status.busy":"2024-12-18T10:25:21.304292Z","iopub.status.idle":"2024-12-18T10:25:26.672750Z","shell.execute_reply":"2024-12-18T10:25:26.671899Z"},"papermill":{"duration":5.37737,"end_time":"2024-12-18T10:25:26.675116","exception":false,"start_time":"2024-12-18T10:25:21.297746","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"id":"d6454066-869f-49fb-91de-c0df86401301","cell_type":"markdown","source":"# TimeSeries Processor","metadata":{}},{"id":"e11b92f9-2952-49cc-9285-3410023516fd","cell_type":"code","source":"# TimeSeries Processor, including:\n\n# Time series data loading and validation\n# Data Preprocessing\n# Time alignment and resampling\n# Missing value handling\n# Outlier handling\n# Batch processing support\n# FFT feature extraction# Band energy analysis\n# Spectral entropy and moments\n# Spectrum slope\n# Activity feature extraction\n# Activity level classification\n# Campaign continuity analysis\n# Time activity mode\n# Social pattern characteristics\n# Social time period analysis\n# Weekday/weekend comparison\n# Social regularity indicator\n# Auxiliary functions\n# Activity duration analysis\n# The longest continuous value of the sequence is calculated\n\n# 时间序列数据加载和验证\n# 数据预处理\n# 时间对齐和重采样\n# 缺少值处理\n# 异常值处理\n# 批处理支持\n# FFT特征提取\n# 频带能量分析\n# 谱熵和谱矩\n# 频谱斜率\n# 活动特征提取\n# 活动水平分类\n# 活动持续性分析\n# 时段活动模式\n# 社交模式特征\n# 社交时间段分析\n# 工作日/周末对比\n# 社交规律性指标\n# 辅助函数\n# 活动持续时间分析\n# 序列最长连续值计算","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"de734b69","cell_type":"code","source":"train_basic_info = pd.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/train.csv')\ntest_basic_info = pd.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/test.csv')\n\n# Train time series data has already been processed to current datasets\n# You can regenerate it just by uncommenting the code in the main method below\ntrain_time_features = pd.read_csv('/kaggle/input/train-ts-features2/features2.csv')\n\ndata_dict = pd.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/data_dictionary.csv')\n\ncsv_output_dest = '/kaggle/working/submission.csv'","metadata":{"execution":{"iopub.execute_input":"2024-12-18T10:25:26.693938Z","iopub.status.busy":"2024-12-18T10:25:26.693400Z","iopub.status.idle":"2024-12-18T10:25:26.884879Z","shell.execute_reply":"2024-12-18T10:25:26.883534Z"},"papermill":{"duration":0.199309,"end_time":"2024-12-18T10:25:26.887218","exception":false,"start_time":"2024-12-18T10:25:26.687909","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"id":"b095270e","cell_type":"code","source":"class TimeSeriesProcessor:\n\n    def __init__(self, config: Optional[Dict] = None):\n        \"\"\"\n        Initialize TimeSeriesProcessor\n        \n        Args:\n            config: Configuration dictionary with processing parameters\n                - min_length: Minimum required length for time series\n                - max_gap: Maximum allowed gap in data points\n                - resample_rate: Rate for resampling time series\n                - validation_rules: Rules for data validation\n        \"\"\"\n        self.config = config or {\n            'min_length': 24 * 60,  # 1 day minimum\n            'max_gap': 30,          # 30 minute maximum gap\n            'resample_rate': '5T',  # 5 minute resampling\n            'validation_rules': {\n                'required_columns': ['X', 'Y', 'Z', 'enmo', 'time_of_day'],\n                'value_ranges': {\n                    'enmo': (0, 100),\n                    'anglez': (-180, 180)\n                }\n            }\n        }\n\n        # Initialize cache for intermediate results\n        self.cache = {}\n\n    def load_parquet_data(self, base_dir: str, batch_size: int = 50) -> Dict[str, pd.DataFrame]:\n        \"\"\"\n        Load parquet files from directory structure with batching\n        \n        Args:\n            base_dir: Base directory containing ID folders\n            batch_size: Number of files to load at once\n            \n        Returns:\n            Dictionary mapping IDs to DataFrames\n        \"\"\"\n        parquet_data = {}\n        base_path = Path(base_dir)\n\n        try:\n            # Get list of all ID directories\n            id_dirs = [d for d in base_path.iterdir() if d.is_dir() and d.name.startswith('id=')]\n\n            # Process in batches\n            for i in tqdm(range(0, len(id_dirs), batch_size), desc=\"Loading parquet files\"):\n                batch_dirs = id_dirs[i:i + batch_size]\n\n                for id_dir in batch_dirs:\n                    # Extract ID from directory name\n                    subject_id = id_dir.name.split('=')[1]\n\n                    # Find and read parquet file\n                    parquet_file = id_dir / 'part-0.parquet'\n                    if parquet_file.exists():\n                        df = pd.read_parquet(parquet_file)\n\n                        # Validate data\n                        if self._validate_timeseries(df):\n                            parquet_data[subject_id] = df\n\n\n                # Clear memory periodically\n                gc.collect()\n\n            return parquet_data\n\n        except Exception as e:\n            return {}\n\n    def _validate_timeseries(self, df: pd.DataFrame) -> bool:\n        \"\"\"Validate time series data structure and content\"\"\"\n        try:\n            # Check required columns\n            required_cols = self.config['validation_rules']['required_columns']\n            if not all(col in df.columns for col in required_cols):\n                return False\n\n            # Check data length\n            if len(df) < self.config['min_length']:\n                return False\n\n            # Check value ranges\n            for col, (min_val, max_val) in self.config['validation_rules']['value_ranges'].items():\n                if col in df.columns:\n                    if df[col].min() < min_val or df[col].max() > max_val:\n                        return False\n\n            return True\n\n        except Exception as e:\n            return False\n\n    def preprocess_timeseries(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"\n        Preprocess time series data with validation and error handling\n        \n        Includes:\n        - Time alignment and resampling\n        - Missing value handling\n        - Outlier detection\n        - Signal filtering\n        \"\"\"\n        try:\n            df = df.copy()\n\n            # Convert time fields\n            df['timestamp'] = pd.to_datetime(df['time_of_day'])\n            df['time_obj'] = df['time_of_day'].apply(self._convert_time)\n\n            # Resample to regular intervals\n            df = df.set_index('timestamp').resample(self.config['resample_rate']).agg({\n                'X': 'mean',\n                'Y': 'mean',\n                'Z': 'mean',\n                'enmo': 'mean',\n                'anglez': 'mean',\n                'time_obj': 'first',\n                'non-wear_flag': 'max'\n            })\n\n            # Handle missing values\n            df = self._handle_missing_values(df)\n\n            # Remove outliers\n            df = self._handle_outliers(df)\n\n            return df\n\n        except Exception as e:\n            raise\n\n    def _convert_time(self, nanoseconds: int) -> Optional[time]:\n        \"\"\"Convert nanosecond timestamp to time of day\"\"\"\n        try:\n            ns = int(nanoseconds)\n            total_seconds = ns // 1_000_000_000\n            total_seconds = total_seconds % 86400\n\n            hours = total_seconds // 3600\n            minutes = (total_seconds % 3600) // 60\n            seconds = total_seconds % 60\n\n            return time(hour=int(hours), minute=int(minutes), second=int(seconds))\n\n        except Exception as e:\n            self.logger.error(f\"Error converting time: {str(e)}\")\n            return None\n\n    def _handle_missing_values(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"Handle missing values in time series data\"\"\"\n        try:\n            # Forward fill small gaps\n            df = df.fillna(method='ffill', limit=self.config['max_gap'])\n\n            # Fill remaining with zeros or interpolate\n            numeric_cols = df.select_dtypes(include=[np.number]).columns\n            df[numeric_cols] = df[numeric_cols].fillna(0)\n\n            return df\n\n        except Exception as e:\n            return df\n\n    def _handle_outliers(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"Remove statistical outliers from time series\"\"\"\n        try:\n            for col in ['X', 'Y', 'Z', 'enmo', 'anglez']:\n                if col in df.columns:\n                    # Calculate z-scores\n                    z_scores = np.abs(stats.zscore(df[col], nan_policy='omit'))\n                    # Replace outliers with column median\n                    df.loc[z_scores > 3, col] = df[col].median()\n\n            return df\n\n        except Exception as e:\n            self.logger.error(f\"Error handling outliers: {str(e)}\")\n            return df\n\n\nclass AutoEncoder(nn.Module):\n    def __init__(self, input_dim, encoding_dim):\n        super(AutoEncoder, self).__init__()\n        self.encoder = nn.Sequential(\n            nn.Linear(input_dim, encoding_dim*3),\n            nn.ReLU(),\n            nn.Linear(encoding_dim*3, encoding_dim*2),\n            nn.ReLU(), \n            nn.Linear(encoding_dim*2, encoding_dim),\n            nn.ReLU()\n        )\n        self.decoder = nn.Sequential(\n            nn.Linear(encoding_dim, input_dim*2),\n            nn.ReLU(),\n            nn.Linear(input_dim*2, input_dim*3),\n            nn.ReLU(),\n            nn.Linear(input_dim*3, input_dim),\n            nn.Sigmoid()\n        )\n        \n    def forward(self, x):\n        encoded = self.encoder(x)\n        decoded = self.decoder(encoded)\n        return decoded\n\nclass FeatureExtractor:\n    \"\"\"Extract features from processed time series data\"\"\"\n\n    def _extract_autoencoder_features(self, df: pd.DataFrame, encoding_dim=60, epochs=100, batch_size=32) -> Dict[str, float]:\n        \"\"\"Extract autoencoder features from time series data\"\"\"\n        features = {}\n        try:\n            # 准备数据\n            scaler = StandardScaler()\n            df_scaled = scaler.fit_transform(df.select_dtypes(include=[np.number]))\n            data_tensor = torch.FloatTensor(df_scaled)\n            \n            # 初始化autoencoder\n            input_dim = data_tensor.shape[1] \n            autoencoder = AutoEncoder(input_dim, encoding_dim)\n        \n            # 训练\n            criterion = nn.MSELoss()\n            optimizer = optim.Adam(autoencoder.parameters())\n            \n            for epoch in range(epochs):\n                for i in range(0, len(data_tensor), batch_size):\n                    batch = data_tensor[i:i + batch_size]\n                    optimizer.zero_grad()\n                    reconstructed = autoencoder(batch)\n                    loss = criterion(reconstructed, batch)\n                    loss.backward()\n                    optimizer.step()\n                    \n            # Extract coded feature\n            with torch.no_grad():\n                encoded_data = autoencoder.encoder(data_tensor).numpy()\n                \n            # Convert to a feature dictionary\n            for i in range(encoded_data.shape[1]):\n                features[f'ae_feature_{i}'] = encoded_data[:,i].mean()\n                \n            return features\n            \n        except Exception as e:\n            self.logger.error(f\"Error extracting autoencoder features: {str(e)}\")\n            return {}\n            \n    def __init__(self, config: Optional[Dict] = None):\n        \"\"\"\n        Initialize FeatureExtractor\n        \n        Args:\n            config: Configuration dictionary with feature extraction parameters\n                - feature_sets: List of feature sets to extract\n                - activity_thresholds: Activity level thresholds\n                - window_sizes: Sizes for different analysis windows\n        \"\"\"\n        self.config = config or {\n            'feature_sets': ['statistical', 'rhythm', 'activity', 'fft', 'social', 'circadian'],\n            'activity_thresholds': {\n                'sedentary': 0.02,\n                'light': 0.1,\n                'moderate': 0.4\n            },\n            'window_sizes': {\n                'short': 300,    # 5 minutes\n                'medium': 1800,  # 30 minutes\n                'long': 3600     # 1 hour\n            }\n        }\n\n\n    def extract_all_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract all configured feature sets from time series data\"\"\"\n        features = {}\n\n        try:\n            # Extract each feature set based on configuration\n            if 'statistical' in self.config['feature_sets']:\n                features.update(self._extract_statistical_features(df))\n\n            if 'rhythm' in self.config['feature_sets']:\n                features.update(self._extract_rhythm_features(df))\n\n            if 'activity' in self.config['feature_sets']:\n                features.update(self._extract_activity_features(df))\n\n            if 'fft' in self.config['feature_sets']:\n                features.update(self._extract_fft_features(df['enmo']))\n\n            if 'social' in self.config['feature_sets']:\n                features.update(self._extract_social_features(df))\n                \n            if 'circadian' in self.config['feature_sets']:\n                features.update(self._extract_circadian_features(df))\n\n            if 'autoencoder' in self.config['feature_sets']:\n                features.update(self._extract_autoencoder_features(df))\n            \n            return features\n\n        except Exception as e:\n            self.logger.error(f\"Error extracting features: {str(e)}\")\n            return {}\n\n    def _extract_statistical_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract statistical features from time series\"\"\"\n        features = {}\n\n        try:\n            # Core statistics for each measurement\n            for col in ['X', 'Y', 'Z', 'enmo', 'anglez']:\n                if col in df.columns:\n                    data = df[col].dropna()\n\n                    # Basic statistics\n                    features.update({\n                        f'{col}_mean': data.mean(),\n                        f'{col}_std': data.std(),\n                        f'{col}_skew': data.skew(),\n                        f'{col}_kurtosis': data.kurtosis(),\n                        f'{col}_rms': np.sqrt(np.mean(np.square(data))),\n                        f'{col}_range': data.max() - data.min()\n                    })\n\n                    # Percentile features\n                    percentiles = [1, 5, 25, 50, 75, 95, 99]\n                    features.update({\n                        f'{col}_p{p}': np.percentile(data, p)\n                        for p in percentiles\n                    })\n\n                    # Additional metrics\n                    features.update({\n                        f'{col}_zero_crossings': ((data[:-1] * data[1:]) < 0).sum(),\n                        f'{col}_above_mean': (data > data.mean()).mean(),\n                        f'{col}_peak_to_peak': data.max() - data.min()\n                    })\n\n            # Inter-axis correlations for accelerometer data\n            if all(col in df.columns for col in ['X', 'Y', 'Z']):\n                corr_matrix = df[['X', 'Y', 'Z']].corr()\n                features.update({\n                    'corr_xy': corr_matrix.loc['X', 'Y'],\n                    'corr_xz': corr_matrix.loc['X', 'Z'],\n                    'corr_yz': corr_matrix.loc['Y', 'Z']\n                })\n\n            return features\n\n        except Exception as e:\n            return {}\n\n    def _extract_rhythm_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract circadian rhythm features using rolling statistics\"\"\"\n        features = {}\n\n        try:\n            # Initialize hourly activity patterns\n            df = df.copy()\n            activity = df['enmo'].fillna(0)\n            hour_groups = df.groupby(df.index.hour)\n\n            # Calculate hourly statistics\n            hourly_mean = hour_groups['enmo'].mean()\n            hourly_std = hour_groups['enmo'].std()\n\n            # Basic rhythm metrics\n            features.update({\n                'rhythm_strength': hourly_mean.std(),\n                'rhythm_stability': 1 - (hourly_std.mean() / (hourly_mean.mean() + 1e-6)),\n                'peak_hour': hourly_mean.idxmax(),\n                'trough_hour': hourly_mean.idxmin(),\n                'peak_trough_ratio': (hourly_mean.max() + 1e-6) / (hourly_mean.min() + 1e-6)\n            })\n\n            # Calculate activity onsets and offsets\n            smoothed = self._smooth_signal(hourly_mean.values)\n            threshold = smoothed.mean()\n\n            # Find transitions\n            above_threshold = smoothed > threshold\n            transitions = np.where(np.diff(above_threshold))[0]\n\n            if len(transitions) >= 2:\n                # Find major activity period\n                durations = []\n                for i in range(0, len(transitions)-1, 2):\n                    if i+1 < len(transitions):\n                        duration = transitions[i+1] - transitions[i]\n                        durations.append((duration, transitions[i], transitions[i+1]))\n\n                if durations:\n                    max_duration, onset, offset = max(durations, key=lambda x: x[0])\n                    features.update({\n                        'activity_onset_hour': onset,\n                        'activity_offset_hour': offset,\n                        'active_duration': max_duration\n                    })\n\n            # Calculate rhythm complexity metrics\n            rolling_24h = activity.rolling(window='24H', center=True)\n            features.update({\n                'rhythm_complexity': len(transitions),\n                'day_to_day_stability': 1 - (rolling_24h.std().mean() /\n                                             (rolling_24h.mean().mean() + 1e-6))\n            })\n\n            return features\n\n        except Exception as e:\n            self.logger.error(f\"Error extracting rhythm features: {str(e)}\")\n            return {}\n\n    def _smooth_signal(self, signal: np.ndarray, window_size: int = 3) -> np.ndarray:\n        \"\"\"Apply centered moving average smoothing\"\"\"\n        if len(signal) < window_size:\n            return signal\n\n        weights = np.ones(window_size) / window_size\n        return np.convolve(signal, weights, mode='same')\n\n    # Continuing FeatureExtractor class...\n    \n    def _extract_fft_features(self, series: Union[pd.Series, np.ndarray], prefix: str = '') -> Dict[str, float]:\n        \"\"\"Extract frequency domain features using FFT\"\"\"\n        features = {}\n    \n        try:\n            # Convert to numpy array if needed\n            data = series.to_numpy() if isinstance(series, pd.Series) else series\n            data = data[~np.isnan(data)]  # Remove NaN values\n    \n            if len(data) < 2:  # Need at least 2 points for FFT\n                return {}\n    \n            # Compute FFT\n            yf = fft(data - np.mean(data))  # Remove DC component\n            xf = fftfreq(len(data), 1/self.config.get('sampling_rate', 0.2))  # 5s intervals by default\n    \n            # Get positive frequencies only\n            pos_mask = xf > 0\n            freqs = xf[pos_mask]\n            power = np.abs(yf[pos_mask])\n    \n            # Normalize power spectrum\n            power_norm = power / np.sum(power)\n    \n            # Frequency bands of interest (in Hz)\n            bands = {\n                'ultradian': (0.0001, 0.001),  # ~2.8-28 hours\n                'circadian': (0.001, 0.0167),   # ~1-24 hours\n                'activity': (0.0167, 0.1)       # ~10min-1hour\n            }\n    \n            # Calculate power in each band\n            for band_name, (low, high) in bands.items():\n                band_mask = (freqs >= low) & (freqs <= high)\n                if np.any(band_mask):\n                    band_power = power_norm[band_mask]\n                    features.update({\n                        f'{prefix}power_{band_name}': np.sum(band_power),\n                        f'{prefix}peak_freq_{band_name}': freqs[band_mask][np.argmax(band_power)]\n                    })\n                else:\n                    features.update({\n                        f'{prefix}power_{band_name}': 0,\n                        f'{prefix}peak_freq_{band_name}': 0\n                    })\n    \n            # Calculate spectral features separately\n            spectral_centroid = np.sum(freqs * power_norm) / np.sum(power_norm)\n            spectral_spread = np.sqrt(np.sum(((freqs - spectral_centroid)**2) * power_norm))\n    \n            # Add spectral features\n            features.update({\n                f'{prefix}spectral_entropy': -np.sum(power_norm * np.log2(power_norm + 1e-10)),\n                f'{prefix}spectral_centroid': spectral_centroid,\n                f'{prefix}spectral_spread': spectral_spread,\n                f'{prefix}spectral_slope': np.polyfit(np.log10(freqs + 1e-10),\n                                                      np.log10(power_norm + 1e-10), 1)[0]\n            })\n    \n            return features\n    \n        except Exception as e:\n            self.logger.error(f\"Error extracting FFT features: {str(e)}\")\n            self.logger.debug(\"FFT extraction error details:\", exc_info=True)\n            # Return empty dict on error to allow processing to continue\n            return {}\n    def _extract_activity_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract activity pattern features\"\"\"\n        features = {}\n\n        try:\n            thresholds = self.config['activity_thresholds']\n            activity = df['enmo'].fillna(0)\n\n            # Activity level classification\n            features.update({\n                'sedentary_time': (activity <= thresholds['sedentary']).mean(),\n                'light_activity_time': ((activity > thresholds['sedentary']) &\n                                        (activity <= thresholds['light'])).mean(),\n                'moderate_activity_time': ((activity > thresholds['light']) &\n                                           (activity <= thresholds['moderate'])).mean(),\n                'vigorous_activity_time': (activity > thresholds['moderate']).mean()\n            })\n\n            # Activity bout analysis\n            for intensity in ['moderate', 'vigorous']:\n                threshold = thresholds[intensity.replace('vigorous', 'moderate')]\n                bouts = self._find_activity_bouts(activity > threshold)\n\n                if bouts:\n                    features.update({\n                        f'{intensity}_bout_count': len(bouts),\n                        f'{intensity}_bout_mean_duration': np.mean(bouts),\n                        f'{intensity}_bout_max_duration': np.max(bouts),\n                        f'{intensity}_bout_total_time': np.sum(bouts)\n                    })\n                else:\n                    features.update({\n                        f'{intensity}_bout_count': 0,\n                        f'{intensity}_bout_mean_duration': 0,\n                        f'{intensity}_bout_max_duration': 0,\n                        f'{intensity}_bout_total_time': 0\n                    })\n\n            # Time-based activity patterns\n            hour = df.index.hour\n            for period, (start, end) in {\n                'morning': (6, 12),\n                'afternoon': (12, 18),\n                'evening': (18, 24),\n                'night': (0, 6)\n            }.items():\n                mask = (hour >= start) & (hour < end)\n                period_activity = activity[mask]\n\n                if len(period_activity) > 0:\n                    features.update({\n                        f'{period}_activity_mean': period_activity.mean(),\n                        f'{period}_activity_peak': period_activity.max(),\n                        f'{period}_active_time': (period_activity > thresholds['sedentary']).mean()\n                    })\n\n            # Activity patterns and transitions\n            features.update({\n                'activity_autocorr_lag1': activity.autocorr(lag=1),\n                'activity_trend': np.polyfit(np.arange(len(activity)), activity.values, 1)[0],\n                'activity_changes': np.sum(np.abs(np.diff(activity)) > thresholds['sedentary']),\n                'longest_inactive_period': self._find_longest_streak(activity <= thresholds['sedentary']),\n                'longest_active_period': self._find_longest_streak(activity > thresholds['sedentary'])\n            })\n\n            return features\n\n        except Exception as e:\n            return {}\n\n    def _find_activity_bouts(self, activity_mask: np.ndarray,\n                             min_duration: int = 10) -> List[int]:\n        \"\"\"Find continuous activity bouts exceeding minimum duration\"\"\"\n        bouts = []\n        current_bout = 0\n\n        for is_active in activity_mask:\n            if is_active:\n                current_bout += 1\n            elif current_bout >= min_duration:\n                bouts.append(current_bout)\n                current_bout = 0\n            else:\n                current_bout = 0\n\n        if current_bout >= min_duration:\n            bouts.append(current_bout)\n\n        return bouts\n\n    def _find_longest_streak(self, condition: np.ndarray) -> int:\n        \"\"\"Find longest streak of True values in boolean array\"\"\"\n        if not len(condition):\n            return 0\n\n        streaks = np.diff(np.where(np.concatenate(([False], condition, [False])))[0])\n        return np.max(streaks) if len(streaks) > 0 else 0\n\n    def _extract_social_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract features related to social activity patterns\"\"\"\n        features = {}\n\n        try:\n            # Define time periods for social activity analysis\n            hour = df.index.hour\n            weekday = df.index.weekday\n\n            social_periods = {\n                'morning_routine': (hour >= 6) & (hour < 9),\n                'work_school': (hour >= 9) & (hour < 17),\n                'evening_social': (hour >= 17) & (hour < 23),\n                'late_night': (hour >= 23) | (hour < 6)\n            }\n\n            activity = df['enmo'].fillna(0)\n\n            # Activity patterns during social periods\n            for period, mask in social_periods.items():\n                period_activity = activity[mask]\n                if len(period_activity) > 0:\n                    features.update({\n                        f'{period}_activity': period_activity.mean(),\n                        f'{period}_variability': period_activity.std(),\n                        f'{period}_active_ratio': (period_activity > self.config['activity_thresholds']['sedentary']).mean()\n                    })\n\n            # Weekend vs weekday patterns\n            is_weekend = (weekday >= 5)\n            features.update({\n                'weekend_activity': activity[is_weekend].mean(),\n                'weekday_activity': activity[~is_weekend].mean(),\n                'weekend_weekday_ratio': (activity[is_weekend].mean() /\n                                          (activity[~is_weekend].mean() + 1e-6))\n            })\n\n            # Social timing features\n            evening_mask = social_periods['evening_social']\n            evening_activity = activity[evening_mask]\n\n            if len(evening_activity) > 0:\n                features.update({\n                    'social_hour_regularity': evening_activity.std() / (evening_activity.mean() + 1e-6),\n                    'social_hour_peak_time': hour[evening_mask][evening_activity.argmax()],\n                    'late_activity_ratio': activity[hour >= 22].mean() / (activity[hour < 22].mean() + 1e-6)\n                })\n\n            return features\n\n        except Exception as e:\n            self.logger.error(f\"Error extracting social features: {str(e)}\")\n            return {}\n\n\n    def _extract_circadian_features(self, df: pd.DataFrame) -> Dict[str, float]:\n        \"\"\"Extract detailed biorhythm features / 提取详细的生物节律特征\"\"\"\n        features = {}\n    \n        try:\n            # data preparation\n            activity = df['enmo'].fillna(0)\n            time_of_day = pd.to_datetime(df.index).time\n    \n            # 24 hour activity mode / 24小时活动模式\n            hour_activity = pd.DataFrame({\n                'activity': activity,\n                'hour': pd.to_datetime(df.index).hour\n            }).groupby('hour')['activity'].agg(['mean', 'std'])\n    \n            # Calculate the circadian rhythm strength / 计算昼夜节律强度\n            cosinor_features = self._compute_cosinor(activity)\n            features.update(cosinor_features)\n    \n            # Calculate the activity transformation characteristics / 计算活动转换特征\n            rest_onset, rest_offset = self._detect_rest_periods(activity)\n            features.update({\n                'rest_onset_hour': rest_onset,\n                'rest_offset_hour': rest_offset,\n                'rest_duration': (rest_offset - rest_onset) % 24\n            })\n    \n            # Calculate the feature of the active fragment / 计算活动片段特征\n            active_bouts = self._compute_bout_features(activity)\n            features.update(active_bouts)\n    \n            # Calculate the diurnal stability index / 计算昼夜稳定性指标\n            is_features = self._compute_IS_IV(activity)\n            features.update(is_features)\n    \n            return features\n    \n        except Exception as e:\n            return {}\n    \n    def _compute_cosinor(self, activity: pd.Series) -> Dict[str, float]:\n        \"\"\"Analyze 24-hour rhythm by Cosinor / 使用Cosinor方法分析24小时节律\"\"\"\n        try:\n            t = np.arange(len(activity)) * (24/(24*60/5))  # 5min sampling\n            w = 2 * np.pi / 24\n    \n            X = np.column_stack((\n                np.cos(w * t),\n                np.sin(w * t),\n                np.ones_like(t)\n            ))\n    \n            Y = activity.values\n            beta = np.linalg.lstsq(X, Y, rcond=None)[0]\n    \n            amplitude = np.sqrt(beta[0]**2 + beta[1]**2)\n            acrophase = np.arctan2(-beta[1], beta[0]) * (24 / (2*np.pi))\n            mesor = beta[2]\n    \n            return {\n                'rhythm_amplitude': amplitude,\n                'rhythm_acrophase': acrophase % 24,\n                'rhythm_mesor': mesor\n            }\n    \n        except Exception as e:\n            self.logger.error(f\"Cosinor分析错误: {str(e)}\")\n            return {}\n    \n    def _detect_rest_periods(self, activity: pd.Series) -> Tuple[float, float]:\n        \"\"\"Detect the rest time window / 检测休息时间窗口\"\"\"\n        try:\n            # Averaging by sliding windows / 使用滑动窗口平均\n            window_size = 12  # 1hour\n            smoothed = activity.rolling(window=window_size, center=True).mean()\n    \n            # Calculate your hourly activity level / 计算每小时活动水平\n            hourly = pd.DataFrame({\n                'activity': smoothed,\n                'hour': pd.to_datetime(activity.index).hour\n            }).groupby('hour')['activity'].mean()\n    \n            # Find the continuous period of lowest activity / 找出活动最低的连续时间段\n            hours = np.arange(24)\n            min_activity_sum = float('inf')\n            rest_onset = 0\n    \n            for i in range(24):\n                window = (hours[i:i+8]) % 24  # 8hours window\n                activity_sum = hourly.iloc[window].sum()\n    \n                if activity_sum < min_activity_sum:\n                    min_activity_sum = activity_sum\n                    rest_onset = i\n    \n            rest_offset = (rest_onset + 8) % 24\n    \n            return rest_onset, rest_offset\n    \n        except Exception as e:\n            return 0, 8\n    \n    def _compute_bout_features(self, activity: pd.Series) -> Dict[str, float]:\n        \"\"\"Calculate the feature of the active fragment / 计算活动片段特征\"\"\"\n        try:\n            threshold = activity.mean() + activity.std()\n            is_active = activity > threshold\n    \n            # 寻找活动片段\n            active_bouts = []\n            current_bout = 0\n    \n            for active in is_active:\n                if active:\n                    current_bout += 1\n                elif current_bout > 0:\n                    active_bouts.append(current_bout)\n                    current_bout = 0\n    \n            if current_bout > 0:\n                active_bouts.append(current_bout)\n    \n            if not active_bouts:\n                return {\n                    'bout_count': 0,\n                    'mean_bout_duration': 0,\n                    'max_bout_duration': 0\n                }\n    \n            return {\n                'bout_count': len(active_bouts),\n                'mean_bout_duration': np.mean(active_bouts) * 5 / 60,  # 转换为小时\n                'max_bout_duration': np.max(active_bouts) * 5 / 60\n            }\n    \n        except Exception as e:\n            self.logger.error(f\"活动片段分析错误: {str(e)}\")\n            return {}\n    \n    def _compute_IS_IV(self, activity: pd.Series) -> Dict[str, float]:\n        \"\"\"Calculate circadian stability (IS) and diurnal variability (IV)\n        计算昼夜节律稳定性(IS)和日间可变性(IV)\"\"\"\n        try:\n            # resample to 1 hour / 重采样到1小时\n            hourly = activity.resample('1H').mean()\n    \n            # Calc IS\n            daily_profile = hourly.groupby(hourly.index.hour).mean()\n            n = len(hourly)\n            p = 24  # period length (hour) / 周期长度(小时)\n    \n            h = daily_profile.values\n            x = hourly.values\n    \n            rss = np.sum((x - np.mean(x))**2)\n            rs = n * np.sum((h - np.mean(x))**2) / p\n    \n            is_value = rs / rss if rss > 0 else 0\n    \n            # Calc IV\n            hourly_diff = np.diff(hourly.values)\n            iv_value = np.sum(hourly_diff**2) / (n * np.var(x)) if np.var(x) > 0 else 0\n    \n            return {\n                'rhythm_stability': is_value,\n                'rhythm_variability': iv_value,\n                'rhythm_strength': is_value / (iv_value + 0.1)  # Comprehensive index / 综合指标\n            }\n    \n        except Exception as e:\n            return {}\n\n\n\nclass FeatureValidator:\n    \"\"\"Validate and clean extracted features\"\"\"\n\n    def __init__(self, config: Optional[Dict] = None):\n        \"\"\"\n        Initialize FeatureValidator\n        \n        Args:\n            config: Configuration dictionary with validation rules\n        \"\"\"\n        self.config = config or {\n            'missing_threshold': 0.3,  # Maximum allowed missing ratio\n            'correlation_threshold': 0.95,  # Threshold for correlation filtering\n            'min_variance_ratio': 1e-8,  # Minimum variance ratio to keep feature\n            'value_ranges': {\n                'activity_features': (0, 1),\n                'correlation_features': (-1, 1),\n                'ratio_features': (0, None),\n                'time_features': (0, 24)\n            }\n        }\n\n    def validate_features(self, features: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"\n        Validate and clean feature DataFrame\n        \n        Args:\n            features: DataFrame of extracted features\n            \n        Returns:\n            Cleaned and validated DataFrame\n        \"\"\"\n        try:\n            df = features.copy()\n\n            # Remove features with too many missing values\n            missing_ratio = df.isnull().mean()\n            df = df.drop(columns=missing_ratio[missing_ratio > self.config['missing_threshold']].index)\n\n            # Remove constant or near-constant features\n            variance = df.var()\n            df = df.drop(columns=variance[variance/variance.max() < self.config['min_variance_ratio']].index)\n\n            # Handle infinite values\n            df = df.replace([np.inf, -np.inf], np.nan)\n\n            # Fill remaining missing values\n            df = self._fill_missing_values(df)\n\n            # Value range validation\n            df = self._validate_value_ranges(df)\n\n            # Remove highly correlated features\n            df = self._remove_correlated_features(df)\n\n            return df\n\n        except Exception as e:\n            return features\n\n    def _fill_missing_values(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"Fill missing values using appropriate strategies\"\"\"\n        try:\n            # Group similar features\n            time_cols = [col for col in df.columns if any(x in col.lower() for x in ['time', 'hour', 'duration'])]\n            ratio_cols = [col for col in df.columns if any(x in col.lower() for x in ['ratio', 'percentage', 'prop'])]\n\n            # Fill missing values by group\n            for cols in [time_cols, ratio_cols]:\n                if cols:\n                    df[cols] = df[cols].fillna(df[cols].median())\n\n            # Fill remaining with median\n            df = df.fillna(df.median())\n\n            return df\n\n        except Exception as e:\n            return df\n\n    def _validate_value_ranges(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"Validate and clip values to acceptable ranges\"\"\"\n        try:\n            for feature_type, (min_val, max_val) in self.config['value_ranges'].items():\n                cols = [col for col in df.columns if feature_type.split('_')[0] in col.lower()]\n                if cols:\n                    if min_val is not None:\n                        df[cols] = df[cols].clip(lower=min_val)\n                    if max_val is not None:\n                        df[cols] = df[cols].clip(upper=max_val)\n\n            return df\n\n        except Exception as e:\n            self.logger.error(f\"Error validating value ranges: {str(e)}\")\n            return df\n\n    def _remove_correlated_features(self, df: pd.DataFrame) -> pd.DataFrame:\n        \"\"\"Remove highly correlated features\"\"\"\n        try:\n            # Calculate correlation matrix\n            corr_matrix = df.corr().abs()\n\n            # Create upper triangle mask\n            upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))\n\n            # Find features to drop\n            to_drop = [column for column in upper.columns\n                       if any(upper[column] > self.config['correlation_threshold'])]\n\n            if to_drop:\n                df = df.drop(columns=to_drop)\n\n            return df\n\n        except Exception as e:\n            self.logger.error(f\"Error removing correlated features: {str(e)}\")\n            return df\n\nclass BatchProcessor:\n    \"\"\"Process large datasets in batches\"\"\"\n\n    def __init__(self,\n                 time_processor: TimeSeriesProcessor,\n                 feature_extractor: FeatureExtractor,\n                 feature_validator: FeatureValidator,\n                 batch_size: int = 50):\n        \"\"\"\n        Initialize BatchProcessor\n        \n        Args:\n            time_processor: TimeSeriesProcessor instance\n            feature_extractor: FeatureExtractor instance\n            feature_validator: FeatureValidator instance\n            batch_size: Number of samples to process in each batch\n        \"\"\"\n        self.time_processor = time_processor\n        self.feature_extractor = feature_extractor\n        self.feature_validator = feature_validator\n        self.batch_size = batch_size\n\n    def process_dataset(self,\n                        parquet_dir: str,\n                        csv_path: Optional[str] = None) -> pd.DataFrame:\n        \"\"\"\n        Process entire dataset in batches\n        \n        Args:\n            parquet_dir: Directory containing parquet files\n            csv_path: Optional path to CSV with additional data\n            \n        Returns:\n            DataFrame with extracted and validated features\n        \"\"\"\n        try:\n            # Load CSV data if provided\n            csv_data = None\n            if csv_path:\n                csv_data = pd.read_csv(csv_path)\n    \n            # Get list of all parquet files\n            parquet_files = self._get_parquet_files(parquet_dir)\n    \n            # Process in batches\n            all_features = []\n            for i in tqdm(range(0, len(parquet_files), self.batch_size),\n                          desc=\"Processing batches\"):\n                # Get current batch\n                batch_files = parquet_files[i:i + self.batch_size]\n    \n                # Process batch\n                batch_features = self._process_batch(batch_files, csv_data)\n                all_features.append(batch_features)\n    \n                # Clear memory\n                gc.collect()\n    \n            # Combine all features\n            features_df = pd.concat(all_features, ignore_index=True)\n    \n            # 使用left join合并CSV数据，保留所有基本信息\n            if csv_data is not None:\n                features_df = pd.merge(\n                    csv_data,           # 左表：基本信息\n                    features_df,        # 右表：特征数据\n                    on='id',            # 关联键\n                    how='left'          # 使用left join保留所有基本信息\n                )\n    \n            return features_df\n    \n        except Exception as e:\n            self.logger.error(f\"Error processing dataset: {str(e)}\")\n            raise\n\n    def _get_parquet_files(self, parquet_dir: str) -> List[Path]:\n        \"\"\"Get list of all parquet files\"\"\"\n        parquet_path = Path(parquet_dir)\n        return [f for f in parquet_path.glob('**/part-0.parquet')]\n\n    def _process_batch(self,\n                       batch_files: List[Path],\n                       csv_data: Optional[pd.DataFrame]) -> pd.DataFrame:\n        \"\"\"Process a batch of files\"\"\"\n        batch_features = []\n    \n        for file_path in batch_files:\n            try:\n                # Get ID from path\n                subject_id = file_path.parent.name.split('=')[1]\n    \n                # Load and process time series\n                df = pd.read_parquet(file_path)\n                df = self.time_processor.preprocess_timeseries(df)\n    \n                # Extract features\n                features = self.feature_extractor.extract_all_features(df)\n                features['id'] = subject_id\n    \n                batch_features.append(features)\n    \n            except Exception as e:\n                continue\n    \n        # 转换为DataFrame\n        return pd.DataFrame(batch_features)\n\n\n# Configuration dictionaries\ntime_processor_config = {\n    'min_length': 24 * 60,  # 1 day minimum\n    'max_gap': 30,          # 30 minute maximum gap\n    'resample_rate': '5T',  # 5 minute resampling\n    'validation_rules': {\n        'required_columns': ['X', 'Y', 'Z', 'enmo', 'time_of_day'],\n        'value_ranges': {\n            'enmo': (0, 100),\n            'anglez': (-180, 180)\n        }\n    }\n}\n\nfeature_extractor_config = {\n    'feature_sets': ['statistical', 'rhythm', 'activity', 'fft', 'social', 'circadian', 'autoencoder'],\n    'sampling_rate': 0.2,  # 5 second intervals\n    'activity_thresholds': {\n        'sedentary': 0.02,\n        'light': 0.1,\n        'moderate': 0.4\n    },\n    'window_sizes': {\n        'short': 300,    # 5 minutes\n        'medium': 1800,  # 30 minutes\n        'long': 3600     # 1 hour\n    }\n}\n\nfeature_validator_config = {\n    'missing_threshold': 0.3,     # Maximum allowed missing ratio\n    'correlation_threshold': 0.95, # Threshold for correlation filtering\n    'min_variance_ratio': 1e-8,   # Minimum variance ratio to keep feature\n    'value_ranges': {\n        'activity_features': (0, 1),\n        'correlation_features': (-1, 1),\n        'ratio_features': (0, None),\n        'time_features': (0, 24)\n    }\n}\n\ndef setup_pipeline(batch_size: int = 50) -> BatchProcessor:\n    \"\"\"Setup the complete processing pipeline\"\"\"\n    # Initialize components\n    time_processor = TimeSeriesProcessor(config=time_processor_config)\n    feature_extractor = FeatureExtractor(config=feature_extractor_config)\n    feature_validator = FeatureValidator(config=feature_validator_config)\n\n    # Create batch processor\n    processor = BatchProcessor(\n        time_processor=time_processor,\n        feature_extractor=feature_extractor,\n        feature_validator=feature_validator,\n        batch_size=batch_size\n    )\n\n    return processor\n\ndef process_dataset(\n        parquet_dir: str,\n        csv_path: Optional[str] = None,\n        output_dir: str = \"output\",\n        batch_size: int = 50\n) -> pd.DataFrame:\n    \"\"\"\n    Process complete dataset with error handling and logging\n    \n    Args:\n        parquet_dir: Directory containing parquet files\n        csv_path: Optional path to CSV with additional data\n        output_dir: Directory to save results\n        batch_size: Number of samples to process in each batch\n    \n    Returns:\n        DataFrame with extracted features\n    \"\"\"\n    try:\n        # Create output directory\n        output_path = Path(output_dir)\n        output_path.mkdir(parents=True, exist_ok=True)\n\n        # Initialize pipeline\n        processor = setup_pipeline(batch_size=batch_size)\n\n        # Process dataset\n        features_df = processor.process_dataset(\n            parquet_dir=parquet_dir,\n            csv_path=csv_path\n        )\n\n        # Save results\n        output_file = output_path / f\"features.csv\"\n        features_df.to_csv(output_file, index=False)\n\n        for group in feature_extractor_config['feature_sets']:\n            cols = [col for col in features_df.columns if group in col]\n\n        return features_df\n\n    except Exception as e:\n        raise\n\n# Example usage\nif __name__ == \"__main__\":\n    # Train time series data has already been processed to current datasets in kaggle '/kaggle/input/train-ts-features2/features2.csv'\n    # train_features = process_dataset(\n    #     parquet_dir=\"/kaggle/input/child-mind-institute-problematic-internet-use/series_train.parquet/\",\n    #     output_dir=\"/kaggle/working/train/\",\n    #     batch_size=50\n    # )\n    \n    # Process test data\n    test_features = process_dataset(\n        parquet_dir=\"/kaggle/input/child-mind-institute-problematic-internet-use/series_test.parquet/\",\n        output_dir=\"/kaggle/working/test/\",\n        batch_size=50\n    )\n","metadata":{"execution":{"iopub.execute_input":"2024-12-18T10:25:26.898606Z","iopub.status.busy":"2024-12-18T10:25:26.898220Z","iopub.status.idle":"2024-12-18T10:25:34.483503Z","shell.execute_reply":"2024-12-18T10:25:34.482391Z"},"papermill":{"duration":7.594197,"end_time":"2024-12-18T10:25:34.485734","exception":false,"start_time":"2024-12-18T10:25:26.891537","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"id":"a9bb1cf3","cell_type":"markdown","source":"# BasicInfo Processing","metadata":{"papermill":{"duration":0.004135,"end_time":"2024-12-18T10:25:34.494414","exception":false,"start_time":"2024-12-18T10:25:34.490279","status":"completed"},"tags":[]}},{"id":"f69c229c","cell_type":"code","source":"def missing(df):\n    \"\"\"\n    Calculate the missing value and proportion of each column / 计算每一列的缺失值及占比\n    \"\"\"\n    missing_number = df.isnull().sum().sort_values(ascending=False)  # 每一列的缺失值求和后降序排序                  \n    missing_percent = (df.isnull().sum() / df.isnull().count()).sort_values(ascending=False)  # 每一列缺失值占比\n    missing_values = pd.concat([missing_number, missing_percent], axis=1,\n                               keys=['Missing_Number', 'Missing_Percent'])  # 合并为一个DataFrame\n    return missing_values\n\n\ntest_basic_info.head()\n\n# 分离有标签和无标签数据\nlabeled_data = train_basic_info[~train_basic_info['sii'].isna()]\nunlabeled_data = train_basic_info[train_basic_info['sii'].isna()]\n\n# 基本信息-特征工程\ndef engineer_features(df, data_dict):\n    # 1. 季节特征处理\n    def process_season_features(df):\n        season_cols = [col for col in df.columns if 'Season' in col]\n        season_map = {'Spring': 1, 'Summer': 2, 'Fall': 3, 'Winter': 4}\n        for col in season_cols:\n            if col in df.columns:\n                df[col] = df[col].map(season_map).fillna(np.nan)\n        return df\n\n    # 2. 体能测试特征\n    def process_fitness_features(df):\n        # 计算总运动时间（分钟）\n        if 'Fitness_Endurance-Time_Mins' in df.columns and 'Fitness_Endurance-Time_Sec' in df.columns:\n            df['Total_Exercise_Time'] = df['Fitness_Endurance-Time_Mins'] + df['Fitness_Endurance-Time_Sec'] / 60\n\n        # 计算各项体能测试的平均分\n        zone_cols = [col for col in df.columns if '_Zone' in col]\n        df['Average_Fitness_Zone'] = df[zone_cols].mean(axis=1)\n\n        # 创建综合体能得分\n        if all(x in df.columns for x in ['FGC-FGC_CU', 'FGC-FGC_PU', 'FGC-FGC_TL']):\n            df['Total_Strength_Score'] = df[['FGC-FGC_CU', 'FGC-FGC_PU', 'FGC-FGC_TL']].mean(axis=1)\n\n        return df\n\n    # 3. 身体测量特征\n    def process_physical_features(df):\n        # BMI分类\n        if 'Physical-BMI' in df.columns:\n            df['BMI_Category'] = pd.cut(df['Physical-BMI'],\n                                        bins=[0, 18.5, 24.9, 29.9, float('inf')],\n                                        labels=[-1, 0, 1, 2])\n            # labels=['Underweight', 'Normal', 'Overweight', 'Obese'])\n\n        # 血压分类\n        if all(x in df.columns for x in ['Physical-Systolic_BP', 'Physical-Diastolic_BP']):\n            df['BP_Category'] = df.apply(lambda x:\n                                         1 if x['Physical-Systolic_BP'] >= 130 or x['Physical-Diastolic_BP'] >= 80\n                                         else 0, axis=1)\n            # 'High' if x['Physical-Systolic_BP'] >= 130 or x['Physical-Diastolic_BP'] >= 80\n            #  else 'Normal', axis=1)\n\n        return df\n\n    # 4. BIA特征\n    def process_bia_features(df):\n        # 计算体脂率分类\n        if 'BIA-BIA_Fat' in df.columns:\n            df['Fat_Category'] = pd.qcut(df['BIA-BIA_Fat'], q=4,\n                                         labels=[-1, 0, 1, 2])\n            # labels=['Low', 'Normal', 'High', 'Very High'])\n\n        # 计算肌肉比例\n        if all(x in df.columns for x in ['BIA-BIA_SMM', 'BIA-BIA_FFM']):\n            df['Muscle_Ratio'] = df['BIA-BIA_SMM'] / df['BIA-BIA_FFM']\n\n        return df\n\n    # 5. 活动水平特征\n    def process_activity_features(df):\n        activity_cols = [col for col in df.columns if 'PAQ' in col and 'Total' in col]\n        df['Average_Activity_Score'] = df[activity_cols].mean(axis=1)\n        return df\n\n    # 6. 缺失值处理\n    def handle_missing_values(df, data_dict):\n        for column in df.columns:\n            col_type = data_dict[data_dict['Field'] == column]['Type'].values\n            if len(col_type) == 0:\n                continue\n\n            if col_type[0] == 'categorical int':\n                df[column].fillna(df[column].mode()[0], inplace=True)\n            elif col_type[0] in ['float', 'int']:\n                df[column].fillna(df[column].median(), inplace=True)\n            elif col_type[0] == 'str':\n                df[column].fillna(df[column].mode()[0], inplace=True)\n        return df\n\n\n    df = process_season_features(df)\n    df = process_fitness_features(df)\n    df = process_physical_features(df)\n    df = process_bia_features(df)\n    df = process_activity_features(df)\n    # df = handle_missing_values(df, data_dict)\n\n    return df\n\n# Apply feature engineering\nlabeled_train_basic_info = engineer_features(labeled_data, data_dict)\n# unlabeled_train_basic_info = engineer_features(unlabeled_data, data_dict)\ntest_basic_info = engineer_features(test_basic_info, data_dict)\n\n# 使用set集合运算获取不重合的列名\nnon_overlapping_ts = set(labeled_train_basic_info.columns) - set(test_basic_info.columns)  # ts中独有的列\nnon_overlapping_y = set(test_basic_info.columns) - set(labeled_train_basic_info.columns)  # y中独有的列\nnon_overlapping_ts -= {'sii'}\nnon_overlapping_ts\n\n# 删除labeled_train_basic_info中test_basic_info中不存在的列\nlabeled_train_basic_info.drop(columns=non_overlapping_ts, inplace=True)\n\n\n# 合并有运动数据的train&test\n\n# 有手环、训练集合并\ntrain_merged = pd.merge(\n    labeled_train_basic_info,  # 基本信息作为左表\n    train_time_features,  # 时序特征作为右表\n    on='id',  # 通过用户ID关联\n    how='left'  # 保留所有基本信息记录\n)\n\n# 测试集合并\ntest_time_features = pd.read_csv('/kaggle/working/test/features.csv')\n\ntest_merged = pd.merge(\n    test_basic_info,  # 基本信息作为左表\n    test_time_features,  # 时序特征作为右表\n    on='id',  # 通过用户ID关联\n    how='left'  # 保留所有基本信息记录\n)\n\nmsing = missing(train_merged)\n\nMISSING_VALUE_THRESHOLD = 0.85  # 缺失值阈值\nCORRELATION_THRESHOLD = 0.95  # 相关性阈值\nVARIANCE_THRESHOLD = 0.0  # 方差阈值\n\n# 1. Removes columns with too many missing values / 删除缺失值过多的列\nmissing_value_cols = msing[msing['Missing_Percent'] > MISSING_VALUE_THRESHOLD].index\ntrain_merged.drop(columns=missing_value_cols, inplace=True)\n\n# Skip the nonexistent column in test_merged / 跳过test_merged中不存在的列\nmissing_value_cols = [col for col in missing_value_cols if col in test_merged.columns]\ntest_merged.drop(columns=missing_value_cols, inplace=True)\n\n# 2. Remove 0 variance features / 删除方差为0的特征\nvariance = train_merged.var(numeric_only=True)\nzero_var_cols = variance[variance <= VARIANCE_THRESHOLD].index\ntrain_merged.drop(columns=zero_var_cols, inplace=True)\n\nzero_var_cols = [col for col in zero_var_cols if col in test_merged.columns]\ntest_merged.drop(columns=zero_var_cols, inplace=True)\n\n# 3. Delete highly related features / 删除高度相关特征\nnumeric_df = train_merged.select_dtypes(include=['int64', 'float64'])\n# Calculate correlation matrix for numeric columns\ncorr_matrix = numeric_df.corr().abs()\n# Create upper triangular matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))\n# Identify columns to drop\nto_drop = [column for column in upper.columns if any(upper[column] > CORRELATION_THRESHOLD)]\n\n# Return dataframe with original columns, dropping only correlated numeric columns\ntrain_merged.drop(columns=to_drop, inplace=True)\nto_drop = [col for col in to_drop if col in test_merged.columns]\ntest_merged.drop(columns=to_drop, inplace=True)\n\n\n# 4. Ensure that the feature column order is the same as train_merged\n# 确保特征列顺序与train_merged一致\ntest_merged = test_merged[train_merged.drop(columns='sii').columns]\n\n# delete ID\ntrain_merged.drop('id', axis=1, inplace=True)\ntrain_merged.to_csv('/kaggle/working/training.csv')\ntest_merged.to_csv('/kaggle/working/test.csv')","metadata":{"execution":{"iopub.execute_input":"2024-12-18T10:25:34.505196Z","iopub.status.busy":"2024-12-18T10:25:34.504785Z","iopub.status.idle":"2024-12-18T10:25:35.465059Z","shell.execute_reply":"2024-12-18T10:25:35.464086Z"},"papermill":{"duration":0.96866,"end_time":"2024-12-18T10:25:35.467364","exception":false,"start_time":"2024-12-18T10:25:34.498704","status":"completed"},"tags":[]},"outputs":[],"execution_count":null},{"id":"ac73ceaa","cell_type":"markdown","source":"# LightGBM","metadata":{"papermill":{"duration":0.004053,"end_time":"2024-12-18T10:25:35.475844","exception":false,"start_time":"2024-12-18T10:25:35.471791","status":"completed"},"tags":[]}},{"id":"d87778d8","cell_type":"code","source":"test_id = test_merged['id']\ntest_merged.drop('id', axis=1, inplace=True)\n\nX_train = train_merged.drop('sii', axis=1)\ny_train = train_merged['sii']\n\nX_train = X_train.rename(columns = lambda x:re.sub('[^A-Za-z0-9_]+', '', x))\ntest_merged = test_merged.rename(columns = lambda x:re.sub('[^A-Za-z0-9_]+', '', x))\n\ndef convert_to_integer_classes(predictions, rounding_strategy='round', min_class=0, max_class=4):\n    \"\"\"\n    Convert float predictions to integer classes with different rounding strategies\n    \n    Args:\n    predictions (array-like): Continuous predictions\n    rounding_strategy (str): 'round' (nearest), 'floor' (down), or 'ceil' (up)\n    min_class (int): Minimum possible class\n    max_class (int): Maximum possible class\n    \n    Returns:\n    numpy array of integer classes\n    \"\"\"\n    if rounding_strategy == 'round':\n        int_preds = np.round(predictions).astype(int)\n    elif rounding_strategy == 'floor':\n        int_preds = np.floor(predictions).astype(int)\n    elif rounding_strategy == 'ceil':\n        int_preds = np.ceil(predictions).astype(int)\n    else:\n        raise ValueError(\"Invalid rounding strategy. Choose 'round', 'floor', or 'ceil'.\")\n\n    # Clip predictions to ensure they're within the specified range\n    int_preds = np.clip(int_preds, min_class, max_class)\n\n    return int_preds\n\na = y_train.mean()\nb = y_train.var(ddof=0)\n\ny_min = y_train.min()\ny_max = y_train.max()\n\ndef quadratic_weighted_kappa(preds, data):\n    y_true = data.get_label()\n    y_pred = preds.clip(y_min, y_max).round()\n    qwk = cohen_kappa_score(y_true, y_pred, weights=\"quadratic\")\n    return 'QWK', qwk, True\n\n\ndef qwk_obj(preds, dtrain):\n    labels = dtrain.get_label()\n    preds = preds.clip(y_min, y_max)\n    f = 1/2 * np.sum((preds - labels)**2)\n    g = 1/2 * np.sum((preds - a)**2 + b)\n    df = preds - labels\n    dg = preds - a\n    grad = (df/g - f*dg/g**2)*len(labels)\n    hess = np.ones(len(labels))\n    return grad, hess\n\n# 修改网格搜索参数，确保所有参数都是列表\nparam_grid = {\n    'rounding_strategy': ['round'],\n    'learning_rate': [0.0092],\n    'max_depth': [4],\n    'objective': [qwk_obj],  # Custom objective functions / 将自定义目标函数包装在列表中\n    # 如果有其他参数，也同样包装在列表中\n}\n\n# 创建参数网格\ngrid_params = list(ParameterGrid(param_grid))\nresults = []\n\n# The initial score of the model is the key parameter.\n# I found that the mean value of the target is a bad choice for this data.\ninit_score = 2.0\n\nfor params in grid_params:\n    # Develop trade-in strategies / 舍入策略选择\n    rounding_strategy = params.pop('rounding_strategy')\n\n    # CV\n    skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)\n    folds = [(idx_train, idx_valid) for idx_train, idx_valid in skf.split(X_train, y_train)]\n\n    # Calc model performance by CV / 使用交叉验证计算模型性能\n    cv_results = lgb.cv(\n        params=params,\n        train_set=lgb.Dataset(X_train, y_train, init_score=[init_score]*len(X_train)),\n        num_boost_round=2000,\n        folds=folds,\n        feval=quadratic_weighted_kappa,\n        callbacks=[\n            lgb.early_stopping(stopping_rounds=100, verbose=True),\n            lgb.log_evaluation(100)\n        ],\n        return_cvbooster=True\n    )\n\n    # Get the average performance of cross-validation\n    performance = np.mean(cv_results['valid QWK-mean'])\n\n    # results\n    results.append({\n        'params': params,\n        'rounding_strategy': rounding_strategy,\n        'performance': performance\n    })\n    \n# find the best param\nbest_result = max(results, key=lambda x: x['performance'])\nprint(\"Best Parameters:\", best_result)\n    \nbest_params = best_result['params']\nbest_rounding_strategy = best_result['rounding_strategy']\n\n# Retrain with the best parameters / 使用最佳参数重新训练模型并生成最终提交结果\nfinal_models = lgb.cv(\n    params=best_params,\n    train_set=lgb.Dataset(X_train, y_train, init_score=[init_score]*len(X_train)),\n    num_boost_round=10000,\n    folds=folds,\n    feval=quadratic_weighted_kappa,\n    callbacks=[\n        lgb.early_stopping(stopping_rounds=100, verbose=True),\n        lgb.log_evaluation(100)\n    ],\n    return_cvbooster=True\n)[\"cvbooster\"].boosters\n\n# Generate final prediction\ntest_preds = np.zeros(len(test_merged))\nfor model in final_models:\n    test_preds += model.predict(test_merged) + init_score\ntest_preds /= len(final_models)\n\n# Transform the prediction using the best rounding strategy / 使用最佳舍入策略转换预测\nres = pd.DataFrame(convert_to_integer_classes(test_preds, rounding_strategy=best_rounding_strategy))\noutput_res = pd.concat([test_id, res], axis=1)\noutput_res.columns = ['id', 'sii']\n\noutput_res.to_csv(csv_output_dest, header=True, index=False)","metadata":{"execution":{"iopub.execute_input":"2024-12-18T10:25:35.487981Z","iopub.status.busy":"2024-12-18T10:25:35.487597Z","iopub.status.idle":"2024-12-18T10:26:30.089512Z","shell.execute_reply":"2024-12-18T10:26:30.088396Z"},"papermill":{"duration":54.61185,"end_time":"2024-12-18T10:26:30.092027","exception":false,"start_time":"2024-12-18T10:25:35.480177","status":"completed"},"tags":[]},"outputs":[],"execution_count":null}]}