import numpy as np
import os
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from tqdm import tqdm
import pandas as pd
from scipy import stats
from scipy.signal import savgol_filter, medfilt
from scipy.ndimage import uniform_filter1d, gaussian_filter1d
from scipy.interpolate import UnivariateSpline, interp1d
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
import warnings
import json
import time
from datetime import datetime
from itertools import product
warnings.filterwarnings('ignore')

def convert_index_to_grid_position(index):
    """Convert linear index to 13x13 grid position (row, col)."""
    row = int(index) // 13
    col = int(index) % 13
    return row, col

def create_grid_outlier_visualization(outlier_indices, param_name):
    """Create 13x13 grid showing outlier locations."""
    grid = np.zeros((13, 13))
    
    for idx in outlier_indices:
        if 0 <= idx < 169:
            row, col = convert_index_to_grid_position(idx)
            grid[row, col] = 1
    
    return grid

def ensure_json_serializable(obj):
    """Convert NumPy types to native Python types for JSON serialization."""
    if isinstance(obj, np.integer):
        return int(obj)
    elif isinstance(obj, np.floating):
        return float(obj)
    elif isinstance(obj, np.ndarray):
        return obj.tolist()
    elif isinstance(obj, list):
        return [ensure_json_serializable(item) for item in obj]
    elif isinstance(obj, tuple):
        return tuple(ensure_json_serializable(item) for item in obj)
    elif isinstance(obj, dict):
        return {key: ensure_json_serializable(value) for key, value in obj.items()}
    else:
        return obj

class AdvancedOutlierDetector:
    """
    Enhanced Advanced Outlier Detection with State Parameter Hyperparameter Tuning
    
    This class now includes the ability to automatically tune state parameter smoothing
    to optimize concordance correlation coefficient while maintaining high Pearson correlation.
    
    UPDATED: Location parameter now uses identical configuration to Alpha parameter.
    """
    
    def __init__(self, save_reports=True):
        # Base configuration - LOCATION NOW MATCHES ALPHA EXACTLY
        self.method_config = {
            'alpha': {
                'methods': ['extreme_threshold', 'robust_zscore'],
                'extreme_threshold': 20.0,
                'smoothing_method': 'lowess',
                'robust_zscore_threshold': 2.5,
                'lowess_frac': 0.3,
                'lowess_iterations': 3,
                'gaussian_sigma': 1.5,
                'spline_smoothing': 0.1,
                'max_outlier_percentage': 15,
                'final_smoothing': True,
                'final_smoothing_method': 'exponential_moving_average',
                'ema_alpha': 0.3
            },
            'beta': {
                'methods': ['extreme_threshold', 'iqr_method'],
                'extreme_threshold': 20.0,
                'smoothing_method': 'gaussian_kernel',
                'iqr_multiplier': 1.5,
                'gaussian_sigma': 2.0,
                'bilateral_sigma_spatial': 2.0,
                'bilateral_sigma_intensity': 0.5,
                'max_outlier_percentage': 18,
                'final_smoothing': True,
                'final_smoothing_method': 'bilateral_filter'
            },
            'location': {
                # CHANGED: Now identical to alpha configuration
                'methods': ['extreme_threshold', 'robust_zscore'],
                'extreme_threshold': 20.0,
                'smoothing_method': 'lowess',
                'robust_zscore_threshold': 2.5,
                'lowess_frac': 0.3,
                'lowess_iterations': 3,
                'gaussian_sigma': 1.5,
                'spline_smoothing': 0.1,
                'max_outlier_percentage': 15,
                'final_smoothing': True,
                'final_smoothing_method': 'exponential_moving_average',
                'ema_alpha': 0.3
            },
            'state': {
                'methods': ['extreme_threshold', 'robust_zscore'],
                'extreme_threshold': 15.0,
                'robust_zscore_threshold': 2.0,
                'smoothing_method': 'bilateral_filter',
                'bilateral_sigma_spatial': 1.5,
                'bilateral_sigma_intensity': 0.3,
                'max_outlier_percentage': 10,
                'final_smoothing': True,
                'final_smoothing_method': 'exponential_moving_average',
                'ema_alpha': 0.15
            }
        }
        
        self.save_reports = save_reports
        self.global_stats = {
            'files_processed': 0,
            'total_outliers_by_param': {'alpha': 0, 'beta': 0, 'location': 0, 'state': 0},
            'total_outliers_by_method': {},
            'smoothing_analysis': {'alpha': [], 'beta': [], 'location': [], 'state': []},
            'processing_log': []
        }

        # Hyperparameter tuning configurations for state parameter
        self.state_param_grid = {
            'gaussian_sigma_moderate': [1.8, 2.0, 2.2, 2.5, 2.8, 3.0],
            'spline_smoothing_factor': [0.25, 0.35, 0.4, 0.5, 0.6, 0.7],
            'spline_degree': [2, 3, 4],
            'ema_passes': [1, 2, 3, 4],
            'ema_alpha': [0.3, 0.4, 0.5, 0.6, 0.7],
            'final_gaussian_sigma': [1.2, 1.5, 1.8, 2.0, 2.3],
            'final_passes': [1, 2],
            'lof_contamination': [0.08, 0.10, 0.11, 0.12, 0.15],
            'lof_neighbors': [5, 7, 9, 11],
            'max_outlier_percentage': [18, 20, 22, 25, 28]
        }

        self.best_state_config = None
        self.tuning_results = []
    
    def concordance_correlation_coefficient(self, y_true, y_pred):
        """
        Calculate concordance correlation coefficient - the key metric we want to optimize.
        
        CCC measures both correlation and agreement, making it crucial for ML prediction quality.
        Unlike Pearson correlation, CCC penalizes systematic bias and scale differences.
        """
        if len(y_true.shape) == 1 or y_true.shape[0] == 1:
            if y_true.shape[0] == 1:
                y_true = y_true.flatten()
                y_pred = y_pred.flatten()
                
            if np.std(y_true) == 0 or np.std(y_pred) == 0:
                return 0
                
            mean_true = np.mean(y_true)
            mean_pred = np.mean(y_pred)
            
            covariance = np.mean((y_true - mean_true) * (y_pred - mean_pred))
            var_true = np.var(y_true, ddof=1)
            var_pred = np.var(y_pred, ddof=1)
            mean_diff_squared = (mean_true - mean_pred) ** 2
            
            numerator = 2 * covariance
            denominator = var_true + var_pred + mean_diff_squared
            
            if denominator == 0:
                return 0
            
            ccc = numerator / denominator
            return ccc
        else:
            mean_true = np.mean(y_true)
            mean_pred = np.mean(y_pred)
            covariance = np.cov(y_true, y_pred)[0, 1]
            var_true = np.var(y_true)
            var_pred = np.var(y_pred)
            mean_diff_squared = (mean_true - mean_pred) ** 2
            
            numerator = 2 * covariance
            denominator = var_true + var_pred + mean_diff_squared
            
            if denominator == 0:
                return 0
            
            ccc = numerator / denominator
            return ccc

    def tune_state_parameters(self, base_dir, num_test_files=8, max_configs=40):
        """
        Hyperparameter tuning function specifically for state parameter smoothing.
        
        This function systematically tests different parameter combinations to find
        the configuration that maximizes concordance correlation coefficient.
        
        Args:
            base_dir: Directory containing ideal and realistic subdirectories
            num_test_files: Number of files to use for evaluation (more = more robust)
            max_configs: Maximum parameter combinations to test (more = more thorough)
        
        Returns:
            best_config: Dictionary with optimal state parameter configuration
            best_score: The concordance score achieved by the best configuration
        """
        print("🎯 HYPERPARAMETER TUNING FOR STATE PARAMETER")
        print("="*60)
        print("Goal: Optimize concordance correlation coefficient")
        print("Focus: State parameter smoothing pipeline only")
        print("="*60)
        
        # Get test files for evaluation
        ideal_dir = os.path.join(base_dir, "ideal_inference")
        all_files = [f for f in os.listdir(ideal_dir) if f.endswith('.txt')]
        
        if len(all_files) < num_test_files:
            num_test_files = len(all_files)
            
        # Randomly select test files for robust evaluation
        test_files = np.random.choice(all_files, num_test_files, replace=False)
        print(f"Selected {len(test_files)} test files for evaluation")
        
        # Generate parameter combinations - prioritize most impactful parameters
        priority_params = ['gaussian_sigma_moderate', 'spline_smoothing_factor', 'ema_alpha', 'final_gaussian_sigma']
        
        # Create focused parameter grid for faster but effective tuning
        focused_grid = {}
        for param in priority_params:
            focused_grid[param] = self.state_param_grid[param]
        
        # Add some secondary parameters with reduced ranges
        focused_grid['ema_passes'] = [2, 3]
        focused_grid['spline_degree'] = [3]
        focused_grid['final_passes'] = [1, 2]
        focused_grid['lof_contamination'] = [0.10, 0.11, 0.12]
        
        # Generate all parameter combinations
        param_names = list(focused_grid.keys())
        param_values = list(focused_grid.values())
        all_combinations = list(product(*param_values))
        
        # Limit to max_configs for reasonable computation time
        if len(all_combinations) > max_configs:
            selected_indices = np.random.choice(len(all_combinations), max_configs, replace=False)
            selected_combinations = [all_combinations[i] for i in selected_indices]
        else:
            selected_combinations = all_combinations
        
        print(f"Testing {len(selected_combinations)} parameter combinations...")
        
        best_score = -1
        best_config = None
        
        # Test each parameter combination
        for i, param_combo in enumerate(tqdm(selected_combinations, desc="Tuning Parameters")):
            # Create configuration dictionary
            test_config = dict(zip(param_names, param_combo))
            
            # Fill in missing parameters with defaults from current config
            full_test_config = self.method_config['state'].copy()
            full_test_config.update(test_config)
            
            # Evaluate this configuration
            score = self.evaluate_state_config(full_test_config, test_files, base_dir)
            
            # Track results
            self.tuning_results.append({
                'config': test_config,
                'full_config': full_test_config,
                'score': score
            })
            
            # Update best configuration if this one is better
            if score > best_score:
                best_score = score
                best_config = full_test_config.copy()
                print(f"\n✅ NEW BEST CONFIG! Score: {score:.4f}")
                print(f"   Key params: σ={test_config.get('gaussian_sigma_moderate', 'default')}, "
                      f"spline={test_config.get('spline_smoothing_factor', 'default')}, "
                      f"ema_α={test_config.get('ema_alpha', 'default')}")
        
        # Store best configuration
        self.best_state_config = best_config
        
        # Print tuning summary
        print(f"\n🏆 TUNING COMPLETE!")
        print(f"Best concordance score: {best_score:.4f}")
        print(f"Configurations tested: {len(selected_combinations)}")
        
        # Show top 3 configurations
        sorted_results = sorted(self.tuning_results, key=lambda x: x['score'], reverse=True)
        print(f"\n📊 TOP 3 CONFIGURATIONS:")
        for i, result in enumerate(sorted_results[:3]):
            config = result['config']
            score = result['score']
            print(f"#{i+1}: Score {score:.4f} - "
                  f"σ={config.get('gaussian_sigma_moderate', 'def')}, "
                  f"spline={config.get('spline_smoothing_factor', 'def')}, "
                  f"ema_α={config.get('ema_alpha', 'def')}")
        
        return best_config, best_score

    def evaluate_state_config(self, test_config, test_files, base_dir):
        """
        Evaluate a specific state parameter configuration on test files.
        
        This function processes state parameters with the given configuration
        and measures the resulting concordance correlation coefficient.
        """
        ideal_dir = os.path.join(base_dir, "ideal_inference")
        realistic_dir = os.path.join(base_dir, "realistic_inference")
        
        all_concordances = []
        
        # Temporarily replace state configuration
        original_config = self.method_config['state'].copy()
        self.method_config['state'] = test_config
        
        try:
            for filename in test_files:
                try:
                    # Load file pair
                    ideal_path = os.path.join(ideal_dir, filename)
                    realistic_path = os.path.join(realistic_dir, filename)
                    
                    with open(ideal_path, 'r') as f:
                        ideal_data = np.array([float(line.strip()) for line in f.readlines()], dtype=np.float64)
                    
                    with open(realistic_path, 'r') as f:
                        realistic_data = np.array([float(line.strip()) for line in f.readlines()], dtype=np.float64)
                    
                    # Ensure same length
                    min_length = min(len(ideal_data), len(realistic_data))
                    ideal_data = ideal_data[:min_length]
                    realistic_data = realistic_data[:min_length]
                    
                    # Extract state parameter indices (every 4th starting from 3)
                    total_dims = len(ideal_data)
                    state_indices = list(range(3, total_dims, 4))
                    
                    # Get state parameter values
                    ideal_state = ideal_data[state_indices]
                    realistic_state = realistic_data[state_indices]
                    
                    # Process state parameters with test configuration
                    processed_ideal_state, _, _, _ = self.process_parameter_data(ideal_state, 'state')
                    processed_realistic_state, _, _, _ = self.process_parameter_data(realistic_state, 'state')
                    
                    # Calculate concordance correlation coefficient
                    concordance = self.concordance_correlation_coefficient(processed_ideal_state, processed_realistic_state)
                    
                    if not np.isnan(concordance) and not np.isinf(concordance):
                        all_concordances.append(concordance)
                        
                except Exception as e:
                    continue  # Skip files that cause errors
            
            # Return average concordance across all test files
            if all_concordances:
                return np.mean(all_concordances)
            else:
                return 0
                
        finally:
            # Always restore original configuration
            self.method_config['state'] = original_config

    def apply_optimized_state_config(self):
        """
        Apply the best state configuration found during tuning.
        
        Call this function after running tune_state_parameters() to use
        the optimized configuration for all subsequent processing.
        """
        if self.best_state_config is not None:
            self.method_config['state'] = self.best_state_config.copy()
            print("✅ Applied optimized state configuration")
            print(f"Key optimized parameters:")
            print(f"  - Gaussian sigma: {self.best_state_config['gaussian_sigma_moderate']}")
            print(f"  - Spline factor: {self.best_state_config['spline_smoothing_factor']}")
            print(f"  - EMA alpha: {self.best_state_config['ema_alpha']}")
            print(f"  - Final sigma: {self.best_state_config['final_gaussian_sigma']}")
        else:
            print("⚠️ No optimized configuration found. Run tune_state_parameters() first.")

    def save_tuning_results(self, output_dir):
        """Save hyperparameter tuning results for analysis."""
        if not os.path.exists(output_dir):
            os.makedirs(output_dir)
        
        # Save detailed results
        results_file = os.path.join(output_dir, f"state_tuning_results_{int(time.time())}.json")
        
        serializable_results = []
        for result in self.tuning_results:
            serializable_result = {
                'config': {k: float(v) if isinstance(v, np.number) else v for k, v in result['config'].items()},
                'score': float(result['score'])
            }
            serializable_results.append(serializable_result)
        
        tuning_summary = {
            'best_config': {k: float(v) if isinstance(v, np.number) else v for k, v in self.best_state_config.items()} if self.best_state_config else None,
            'best_score': float(max(r['score'] for r in self.tuning_results)) if self.tuning_results else 0,
            'all_results': serializable_results,
            'total_configs_tested': len(self.tuning_results)
        }
        
        with open(results_file, 'w') as f:
            json.dump(tuning_summary, f, indent=2)
        
        print(f"💾 Tuning results saved to: {results_file}")

    # All your existing methods remain exactly the same
    def detect_extreme_threshold_outliers(self, data, config):
        """Detect extreme outliers using absolute threshold."""
        threshold = config['extreme_threshold']
        outlier_mask = np.abs(data) > threshold
        outlier_indices = np.where(outlier_mask)[0].astype(int)
        
        report = {
            'method': 'extreme_threshold',
            'threshold': threshold,
            'outliers_found': len(outlier_indices),
            'outlier_indices': [int(x) for x in outlier_indices],
            'outlier_values': [float(data[int(idx)]) for idx in outlier_indices] if len(outlier_indices) > 0 else [],
            'max_positive_value': float(np.max(data)),
            'min_negative_value': float(np.min(data)),
            'data_range': float(np.max(data) - np.min(data)),
            'percentage_removed': float(len(outlier_indices) / len(data) * 100)
        }
        
        return outlier_indices, report
    
    def detect_robust_zscore_outliers(self, data, config):
        """Robust Z-score using Median Absolute Deviation (MAD)."""
        median = np.median(data)
        mad = np.median(np.abs(data - median))
        
        if mad == 0:
            mad = np.mean(np.abs(data - median))
            if mad == 0:
                return np.array([], dtype=int), {'method': 'robust_zscore', 'outliers_found': 0, 'error': 'No variation in data'}
        
        modified_z_scores = 0.6745 * (data - median) / mad
        threshold = config['robust_zscore_threshold']
        outlier_indices = np.where(np.abs(modified_z_scores) > threshold)[0].astype(int)
        
        report = {
            'method': 'robust_zscore',
            'threshold': threshold,
            'median': float(median),
            'mad': float(mad),
            'max_modified_zscore': float(np.max(np.abs(modified_z_scores))),
            'outliers_found': len(outlier_indices),
            'outlier_indices': [int(x) for x in outlier_indices],
            'outlier_zscores': [float(modified_z_scores[int(idx)]) for idx in outlier_indices] if len(outlier_indices) > 0 else []
        }
        
        return outlier_indices, report
    
    def detect_iqr_outliers(self, data, config):
        """IQR method for outlier detection."""
        Q1 = np.percentile(data, 25)
        Q3 = np.percentile(data, 75)
        IQR = Q3 - Q1
        
        multiplier = config['iqr_multiplier']
        lower_bound = Q1 - multiplier * IQR
        upper_bound = Q3 + multiplier * IQR
        
        outlier_mask = (data < lower_bound) | (data > upper_bound)
        outlier_indices = np.where(outlier_mask)[0].astype(int)
        
        report = {
            'method': 'iqr_method',
            'Q1': float(Q1),
            'Q3': float(Q3),
            'IQR': float(IQR),
            'multiplier': multiplier,
            'lower_bound': float(lower_bound),
            'upper_bound': float(upper_bound),
            'outliers_found': len(outlier_indices),
            'outlier_indices': [int(x) for x in outlier_indices],
            'outlier_values': [float(data[int(idx)]) for idx in outlier_indices] if len(outlier_indices) > 0 else []
        }
        
        return outlier_indices, report
    
    def detect_local_outlier_factor_outliers(self, data, config):
        """Local Outlier Factor (LOF) for outlier detection."""
        if len(data) < config['lof_neighbors'] + 1:
            return np.array([], dtype=int), {'method': 'local_outlier_factor', 'outliers_found': 0, 'error': 'Insufficient data for LOF'}
        
        data_reshaped = data.reshape(-1, 1)
        n_neighbors = min(config['lof_neighbors'], len(data) - 1)
        contamination = config['lof_contamination']
        
        lof = LocalOutlierFactor(n_neighbors=n_neighbors, contamination=contamination)
        outlier_labels = lof.fit_predict(data_reshaped)
        outlier_indices = np.where(outlier_labels == -1)[0].astype(int)
        lof_scores = -lof.negative_outlier_factor_
        
        report = {
            'method': 'local_outlier_factor',
            'n_neighbors': n_neighbors,
            'contamination': contamination,
            'outliers_found': len(outlier_indices),
            'outlier_indices': [int(x) for x in outlier_indices],
            'lof_scores': [float(lof_scores[int(idx)]) for idx in outlier_indices] if len(outlier_indices) > 0 else [],
            'mean_lof_score': float(np.mean(lof_scores)),
            'max_lof_score': float(np.max(lof_scores))
        }
        
        return outlier_indices, report
    
    def apply_minimal_interpolation_smoothing(self, data, config):
        """Minimal interpolation for location parameter."""
        smoothed = data.copy()
        threshold = config['extreme_threshold']
        extreme_indices = np.where(np.abs(data) > threshold)[0]
        
        for idx in extreme_indices:
            left_good = None
            right_good = None
            
            for i in range(idx-1, max(0, idx-10), -1):
                if np.abs(data[i]) <= threshold:
                    left_good = data[i]
                    break
            
            for i in range(idx+1, min(len(data), idx+11)):
                if np.abs(data[i]) <= threshold:
                    right_good = data[i]
                    break
            
            if left_good is not None and right_good is not None:
                smoothed[idx] = (left_good + right_good) / 2
            elif left_good is not None:
                smoothed[idx] = left_good
            elif right_good is not None:
                smoothed[idx] = right_good
            else:
                smoothed[idx] = np.sign(data[idx]) * threshold
        
        return smoothed
    
    def apply_lowess_smoothing(self, data, config):
        """LOWESS smoothing with enhanced parameters."""
        try:
            from statsmodels.nonparametric.smoothers_lowess import lowess
            
            x = np.arange(len(data))
            frac = config['lowess_frac']
            it = config.get('lowess_iterations', 3)
            
            smoothed = lowess(data, x, frac=frac, it=it, return_sorted=False)
            return smoothed
            
        except ImportError:
            print("Warning: statsmodels not available, using moving average fallback")
            window = max(int(len(data) * config['lowess_frac']), 3)
            return uniform_filter1d(data, size=window, mode='nearest')
    
    def apply_gaussian_kernel_smoothing(self, data, config):
        """Gaussian kernel smoothing."""
        sigma = config['gaussian_sigma']
        return gaussian_filter1d(data, sigma=sigma, mode='nearest')
    
    def apply_bilateral_filter(self, data, config):
        """Enhanced bilateral filter for edge-preserving smoothing."""
        sigma_spatial = config['bilateral_sigma_spatial']
        sigma_intensity = config['bilateral_sigma_intensity']
        
        smoothed = data.copy()
        
        for i in range(len(data)):
            window_size = int(3 * sigma_spatial)
            start = max(0, i - window_size)
            end = min(len(data), i + window_size + 1)
            
            spatial_weights = np.exp(-0.5 * ((np.arange(start, end) - i) / sigma_spatial) ** 2)
            intensity_weights = np.exp(-0.5 * ((data[start:end] - data[i]) / sigma_intensity) ** 2)
            
            weights = spatial_weights * intensity_weights
            weights /= np.sum(weights)
            
            smoothed[i] = np.sum(weights * data[start:end])
        
        return smoothed
    
    def apply_balanced_smooth_pipeline(self, data, config):
        """Balanced Smooth Pipeline for State Parameter Processing - now tunable."""
        stage1_data = self.apply_moderate_gaussian_smoothing(data, config)
        stage2_data = self.apply_balanced_spline_fitting(stage1_data, config)
        stage3_data = self.apply_light_ema(stage2_data, config)
        return stage3_data
    
    def apply_moderate_gaussian_smoothing(self, data, config):
        """Moderate Gaussian smoothing - tunable parameter."""
        sigma = config.get('gaussian_sigma_moderate', 2.2)
        return gaussian_filter1d(data, sigma=sigma, mode='nearest')
    
    def apply_balanced_spline_fitting(self, data, config):
        """Balanced spline fitting - tunable parameters."""
        smoothing_factor = config.get('spline_smoothing_factor', 0.4)
        degree = config.get('spline_degree', 3)
        
        x = np.arange(len(data))
        
        try:
            spline = UnivariateSpline(x, data, s=smoothing_factor * len(data), k=degree)
            return spline(x)
        except:
            return gaussian_filter1d(data, sigma=2.0, mode='nearest')
    
    def apply_light_ema(self, data, config):
        """Light exponential moving average - tunable parameters."""
        passes = config.get('ema_passes', 2)
        alpha = config.get('ema_alpha', 0.5)
        
        smoothed = data.copy()
        
        for pass_num in range(passes):
            temp_smoothed = np.zeros_like(smoothed)
            temp_smoothed[0] = smoothed[0]
            
            for i in range(1, len(smoothed)):
                temp_smoothed[i] = alpha * smoothed[i] + (1 - alpha) * temp_smoothed[i-1]
            
            if pass_num == passes - 1:
                temp_smoothed[-1] = smoothed[-1]
                for i in range(len(smoothed) - 2, -1, -1):
                    temp_smoothed[i] = 0.7 * temp_smoothed[i] + 0.3 * temp_smoothed[i+1]
            
            smoothed = temp_smoothed
        
        return smoothed
    
    def apply_gentle_gaussian_final(self, data, config):
        """Gentle final Gaussian smoothing - tunable parameters."""
        sigma = config.get('final_gaussian_sigma', 1.8)
        passes = config.get('final_passes', 1)
        
        smoothed = data.copy()
        
        for _ in range(passes):
            smoothed = gaussian_filter1d(smoothed, sigma=sigma, mode='nearest')
        
        return smoothed
    
    def apply_exponential_moving_average(self, data, config):
        """Exponential Moving Average smoothing."""
        alpha = config['ema_alpha']
        smoothed = np.zeros_like(data)
        smoothed[0] = data[0]
        
        for i in range(1, len(data)):
            smoothed[i] = alpha * data[i] + (1 - alpha) * smoothed[i-1]
        
        return smoothed
    
    def detect_outliers_multi_method(self, data, param_name):
        """Apply multiple outlier detection methods in sequence."""
        config = self.method_config[param_name]
        all_outlier_indices = set()
        method_reports = []
        
        extreme_outliers, extreme_report = self.detect_extreme_threshold_outliers(data, config)
        all_outlier_indices.update(extreme_outliers)
        method_reports.append(extreme_report)
        
        cleaned_data = self.replace_extreme_outliers(data, extreme_outliers, config)
        
        method_name = 'extreme_threshold'
        if method_name not in self.global_stats['total_outliers_by_method']:
            self.global_stats['total_outliers_by_method'][method_name] = 0
        self.global_stats['total_outliers_by_method'][method_name] += len(extreme_outliers)
        
        secondary_methods = config['methods'][1:]
        
        for method_name in secondary_methods:
            try:
                if method_name == 'robust_zscore':
                    outliers, report = self.detect_robust_zscore_outliers(cleaned_data, config)
                elif method_name == 'iqr_method':
                    outliers, report = self.detect_iqr_outliers(cleaned_data, config)
                elif method_name == 'local_outlier_factor':
                    outliers, report = self.detect_local_outlier_factor_outliers(cleaned_data, config)
                else:
                    continue
                
                all_outlier_indices.update(outliers)
                method_reports.append(report)
                
                if method_name not in self.global_stats['total_outliers_by_method']:
                    self.global_stats['total_outliers_by_method'][method_name] = 0
                self.global_stats['total_outliers_by_method'][method_name] += len(outliers)
                
            except Exception as e:
                error_report = {
                    'method': method_name,
                    'error': str(e),
                    'outliers_found': 0
                }
                method_reports.append(error_report)
        
        final_outlier_indices = np.array(list(all_outlier_indices), dtype=int)
        outlier_percentage = (len(final_outlier_indices) / len(data)) * 100
        
        if outlier_percentage > config['max_outlier_percentage']:
            final_outliers = list(extreme_outliers)
            remaining_slots = int(len(data) * config['max_outlier_percentage'] / 100) - len(final_outliers)
            
            if remaining_slots > 0:
                other_outliers = [idx for idx in final_outlier_indices if idx not in extreme_outliers]
                final_outliers.extend(other_outliers[:remaining_slots])
            
            final_outlier_indices = np.array(final_outliers, dtype=int)
        
        self.global_stats['total_outliers_by_param'][param_name] += len(final_outlier_indices)
        
        return final_outlier_indices, method_reports, cleaned_data
    
    def replace_extreme_outliers(self, data, outlier_indices, config):
        """Replace extreme outliers with interpolated values."""
        if len(outlier_indices) == 0:
            return data.copy()
        
        cleaned_data = data.copy()
        threshold = config['extreme_threshold']
        
        for idx in outlier_indices:
            idx = int(idx)
            left_good = None
            right_good = None
            
            for i in range(idx-1, max(0, idx-5), -1):
                if np.abs(data[i]) <= threshold:
                    left_good = data[i]
                    break
            
            for i in range(idx+1, min(len(data), idx+6)):
                if np.abs(data[i]) <= threshold:
                    right_good = data[i]
                    break
            
            if left_good is not None and right_good is not None:
                cleaned_data[idx] = (left_good + right_good) / 2
            elif left_good is not None:
                cleaned_data[idx] = left_good
            elif right_good is not None:
                cleaned_data[idx] = right_good
            else:
                cleaned_data[idx] = np.sign(data[idx]) * threshold
        
        return cleaned_data
    
    def apply_smoothing_method(self, data, config):
        """Apply the configured smoothing method."""
        method = config['smoothing_method']
        
        try:
            if method == 'minimal_interpolation':
                return self.apply_minimal_interpolation_smoothing(data, config)
            elif method == 'lowess':
                return self.apply_lowess_smoothing(data, config)
            elif method == 'gaussian_kernel':
                return self.apply_gaussian_kernel_smoothing(data, config)
            elif method == 'bilateral_filter':
                return self.apply_bilateral_filter(data, config)
            elif method == 'balanced_smooth_pipeline':
                return self.apply_balanced_smooth_pipeline(data, config)
            else:
                print(f"Unknown smoothing method '{method}', using Gaussian fallback")
                return gaussian_filter1d(data, sigma=2.0, mode='nearest')
                
        except Exception as e:
            print(f"Warning: {method} smoothing failed: {e}. Using Gaussian fallback.")
            return gaussian_filter1d(data, sigma=2.0, mode='nearest')
    
    def apply_final_smoothing(self, data, config):
        """Apply final smoothing step."""
        if not config.get('final_smoothing', False):
            return data
        
        method = config.get('final_smoothing_method', 'exponential_moving_average')
        
        try:
            if method == 'exponential_moving_average':
                return self.apply_exponential_moving_average(data, config)
            elif method == 'bilateral_filter':
                return self.apply_bilateral_filter(data, config)
            elif method == 'lowess':
                return self.apply_lowess_smoothing(data, config)
            elif method == 'gaussian_kernel':
                return self.apply_gaussian_kernel_smoothing(data, config)
            elif method == 'gentle_gaussian':
                return self.apply_gentle_gaussian_final(data, config)
            else:
                return data
                
        except Exception as e:
            print(f"Warning: Final smoothing method {method} failed: {e}")
            return data
    
    def process_parameter_data(self, data, param_name):
        """Complete processing pipeline for a single parameter."""
        config = self.method_config[param_name]
        
        outlier_indices, method_reports, cleaned_data = self.detect_outliers_multi_method(data, param_name)
        smoothed_data = self.apply_smoothing_method(cleaned_data, config)
        
        processed_data = cleaned_data.copy()
        if len(outlier_indices) > 0:
            processed_data[outlier_indices] = smoothed_data[outlier_indices]
        
        final_data = self.apply_final_smoothing(processed_data, config)
        
        original_std = np.std(data)
        final_std = np.std(final_data)
        noise_reduction = (1 - final_std / original_std) if original_std > 0 else 0
        
        smoothing_report = {
            'original_std': float(original_std),
            'final_std': float(final_std),
            'noise_reduction': float(noise_reduction),
            'smoothing_method': config['smoothing_method'],
            'final_smoothing_method': config.get('final_smoothing_method', 'none'),
            'outliers_replaced': len(outlier_indices),
            'outlier_percentage': float(len(outlier_indices) / len(data) * 100)
        }
        
        return final_data, outlier_indices, method_reports, smoothing_report

# Keep all your existing processing functions exactly the same
def process_file_pair_advanced(ideal_file_path, realistic_file_path, detector):
    """Process a pair of ideal and realistic files using enhanced methods."""
    
    def read_file(file_path):
        with open(file_path, 'r') as f:
            lines = f.readlines()
        
        data = []
        for line in lines:
            try:
                value = float(line.strip())
                if np.isnan(value) or np.isinf(value):
                    value = 0.0
                data.append(value)
            except ValueError:
                data.append(0.0)
        return np.array(data, dtype=np.float64)
    
    ideal_data = read_file(ideal_file_path)
    realistic_data = read_file(realistic_file_path)
    
    min_length = min(len(ideal_data), len(realistic_data))
    ideal_data = ideal_data[:min_length]
    realistic_data = realistic_data[:min_length]
    
    total_dims = len(ideal_data)
    
    parameter_indices = {
        'alpha': list(range(0, total_dims, 4)),
        'beta': list(range(1, total_dims, 4)),
        'location': list(range(2, total_dims, 4)),
        'state': list(range(3, total_dims, 4))
    }
    
    original_ideal = ideal_data.copy()
    original_realistic = realistic_data.copy()
    processed_ideal = ideal_data.copy()
    processed_realistic = realistic_data.copy()
    
    file_report = {
        'ideal_file_path': ideal_file_path,
        'realistic_file_path': realistic_file_path,
        'total_data_points': int(total_dims),
        'parameters': {}
    }
    
    for param_name, indices in parameter_indices.items():
        ideal_param_values = ideal_data[indices]
        realistic_param_values = realistic_data[indices]
        
        proc_ideal, ideal_outliers, ideal_method_reports, ideal_smoothing_report = detector.process_parameter_data(
            ideal_param_values, param_name
        )
        
        proc_realistic, realistic_outliers, realistic_method_reports, realistic_smoothing_report = detector.process_parameter_data(
            realistic_param_values, param_name
        )
        
        combined_outliers = np.unique(np.concatenate([ideal_outliers, realistic_outliers])).astype(int)
        
        for i, idx in enumerate(indices):
            processed_ideal[idx] = proc_ideal[i]
            processed_realistic[idx] = proc_realistic[i]
        
        param_report = {
            'original_count': int(len(ideal_param_values)),
            'combined_outlier_count': int(len(combined_outliers)),
            'ideal_outlier_count': int(len(ideal_outliers)),
            'realistic_outlier_count': int(len(realistic_outliers)),
            'outlier_percentage': float(len(combined_outliers) / len(ideal_param_values) * 100),
            'outlier_indices_in_param': [int(x) for x in combined_outliers],
            'outlier_global_indices': [int(indices[int(i)]) for i in combined_outliers],
            'original_ideal_range': [float(np.min(ideal_param_values)), float(np.max(ideal_param_values))],
            'original_realistic_range': [float(np.min(realistic_param_values)), float(np.max(realistic_param_values))],
            'processed_ideal_range': [float(np.min(proc_ideal)), float(np.max(proc_ideal))],
            'processed_realistic_range': [float(np.min(proc_realistic)), float(np.max(proc_realistic))],
            'ideal_method_reports': ideal_method_reports,
            'realistic_method_reports': realistic_method_reports,
            'ideal_smoothing_report': ideal_smoothing_report,
            'realistic_smoothing_report': realistic_smoothing_report,
            'methods_used': detector.method_config[param_name]['methods'],
            'smoothing_method': detector.method_config[param_name]['smoothing_method']
        }
        
        file_report['parameters'][param_name] = param_report
    
    result = {
        'original_ideal': original_ideal,
        'original_realistic': original_realistic,
        'processed_ideal': processed_ideal,
        'processed_realistic': processed_realistic,
        'file_report': file_report
    }
    
    return result

def create_advanced_outlier_plot(file_name, original_ideal, original_realistic, 
                                processed_ideal, processed_realistic, file_report):
    """Create comprehensive plots showing enhanced outlier detection results."""
    parameter_indices = {
        'Alpha': list(range(0, len(original_ideal), 4)),
        'Beta': list(range(1, len(original_ideal), 4)),
        'Location': list(range(2, len(original_ideal), 4)),
        'State': list(range(3, len(original_ideal), 4))
    }
    
    figs = []
    
    fig1, axes1 = plt.subplots(2, 2, figsize=(20, 16))
    fig1.suptitle(f'ALPHA-STYLE Location + OPTIMIZED State Processing: {file_name}\n'
                 f'(Location: Alpha-style LOWESS | State: Hyperparameter Tuned)', 
                 fontsize=16, fontweight='bold')
    
    axes1 = axes1.flatten()
    
    for i, (param_name, indices) in enumerate(parameter_indices.items()):
        ax = axes1[i]
        param_key = param_name.lower()
        
        orig_ideal = np.array([original_ideal[idx] for idx in indices])
        orig_realistic = np.array([original_realistic[idx] for idx in indices])
        proc_ideal = np.array([processed_ideal[idx] for idx in indices])
        proc_realistic = np.array([processed_realistic[idx] for idx in indices])
        
        param_report = file_report['parameters'][param_key]
        
        x_indices = np.arange(len(orig_ideal))
        
        ax.plot(x_indices, orig_ideal, 'o-', color='lightblue', alpha=0.6, 
               markersize=2, linewidth=0.8, label='Original Ideal')
        ax.plot(x_indices, orig_realistic, 'o-', color='lightcoral', alpha=0.6, 
               markersize=2, linewidth=0.8, label='Original Realistic')
        
        ax.plot(x_indices, proc_ideal, '-', color='darkblue', alpha=0.9, 
               linewidth=3, label='Processed Ideal')
        ax.plot(x_indices, proc_realistic, '-', color='darkred', alpha=0.9, 
               linewidth=3, label='Processed Realistic')
        
        outlier_indices = param_report['outlier_indices_in_param']
        if len(outlier_indices) > 0:
            outlier_ideal_values = [orig_ideal[idx] for idx in outlier_indices]
            outlier_realistic_values = [orig_realistic[idx] for idx in outlier_indices]
            
            ax.scatter(outlier_indices, outlier_ideal_values, 
                      color='orange', s=80, marker='X', linewidth=3, 
                      label=f'Outliers: {len(outlier_indices)}', zorder=10)
            ax.scatter(outlier_indices, outlier_realistic_values, 
                      color='orange', s=80, marker='X', linewidth=3, zorder=10)
        
        methods_used = param_report['methods_used']
        smoothing_method = param_report['smoothing_method']
        noise_reduction_ideal = param_report['ideal_smoothing_report']['noise_reduction']
        noise_reduction_realistic = param_report['realistic_smoothing_report']['noise_reduction']
        
        if param_name.lower() == 'state':
            title = f'{param_name} - HYPERPARAMETER OPTIMIZED 🎯\n'
            title += f'Tuned for Maximum Concordance | Methods: {", ".join(methods_used)}'
        elif param_name.lower() == 'location':
            title = f'{param_name} - ALPHA-STYLE PROCESSING 🔄\n'
            title += f'Identical to Alpha | Methods: {", ".join(methods_used)} | LOWESS Smoothing'
        else:
            title = f'{param_name} - Methods: {", ".join(methods_used)}\n'
            title += f'Smoothing: {smoothing_method}, Noise Reduction: I:{noise_reduction_ideal:.1%} R:{noise_reduction_realistic:.1%}'
        
        ax.set_title(title, fontweight='bold', fontsize=10)
        ax.set_xlabel('Index')
        ax.set_ylabel('Value')
        ax.grid(True, alpha=0.3)
        ax.legend(fontsize=8)
        
        all_processed_values = np.concatenate([proc_ideal, proc_realistic])
        y_min = np.percentile(all_processed_values, 2)
        y_max = np.percentile(all_processed_values, 98)
        y_range = y_max - y_min
        margin = y_range * 0.1
        ax.set_ylim(y_min - margin, y_max + margin)
        
        if param_name.lower() == 'state':
            box_color = 'lightgreen'
            status = "OPTIMIZED 🎯"
        elif param_name.lower() == 'location':
            box_color = 'lightcyan'
            status = "ALPHA-STYLE 🔄"
        elif param_name.lower() == 'alpha':
            box_color = 'lightyellow'
            status = "REFERENCE ✓"
        else:
            box_color = 'lightyellow'
            status = "STANDARD ✓"
        
        methods_text = f"{status}\nMETHODS USED:\n"
        for method_report in param_report['ideal_method_reports']:
            method_name = method_report['method']
            outliers_found = method_report['outliers_found']
            methods_text += f"• {method_name}: {outliers_found} outliers\n"
        
        methods_text += f"\nSMOOTHING:\n"
        methods_text += f"• Primary: {smoothing_method}\n"
        
        final_method = param_report['ideal_smoothing_report'].get('final_smoothing_method', 'none')
        if final_method != 'none':
            methods_text += f"• Final: {final_method}\n"
        
        if param_name.lower() == 'location':
            methods_text += f"• IDENTICAL TO ALPHA\n"
        
        methods_text += f"\nRESULTS:\n"
        methods_text += f"• Outliers: {param_report['combined_outlier_count']}\n"
        methods_text += f"• Noise reduction: {noise_reduction_ideal:.1%}"
        
        ax.text(0.02, 0.98, methods_text, transform=ax.transAxes, fontsize=8,
               verticalalignment='top', bbox=dict(boxstyle='round,pad=0.4', 
               facecolor=box_color, alpha=0.8))
    
    plt.tight_layout()
    plt.subplots_adjust(top=0.88, bottom=0.08)
    figs.append(fig1)
    
    fig2, axes2 = plt.subplots(2, 2, figsize=(20, 16))
    fig2.suptitle(f'13x13 Grid Outlier Locations (ALPHA-STYLE Location): {file_name}', 
                 fontsize=16, fontweight='bold')
    
    axes2 = axes2.flatten()
    
    for i, (param_name, indices) in enumerate(parameter_indices.items()):
        ax = axes2[i]
        param_key = param_name.lower()
        param_report = file_report['parameters'][param_key]
        outlier_indices = param_report['outlier_indices_in_param']
        
        grid = create_grid_outlier_visualization(outlier_indices, param_name)
        
        im = ax.imshow(grid, cmap='Reds', aspect='equal', interpolation='nearest', 
                      vmin=0, vmax=1, origin='lower')
        
        ax.set_xticks(np.arange(-0.5, 13, 1), minor=True)
        ax.set_yticks(np.arange(-0.5, 13, 1), minor=True)
        ax.grid(which='minor', color='gray', linestyle='-', linewidth=0.5, alpha=0.7)
        
        ax.set_xticks(np.arange(0, 13, 2))
        ax.set_yticks(np.arange(0, 13, 2))
        ax.set_xlabel('X Grid Position')
        ax.set_ylabel('Y Grid Position')
        
        outlier_count = len(outlier_indices)
        methods_str = "+".join(param_report['methods_used'])
        
        if param_name.lower() == 'state':
            title = f'{param_name} Grid - OPTIMIZED 🎯 ({outlier_count} outliers)\nMethods: {methods_str} (Hyperparameter Tuned)'
        elif param_name.lower() == 'location':
            title = f'{param_name} Grid - ALPHA-STYLE 🔄 ({outlier_count} outliers)\nMethods: {methods_str} (Identical to Alpha)'
        elif param_name.lower() == 'alpha':
            title = f'{param_name} Grid - REFERENCE ({outlier_count} outliers)\nMethods: {methods_str} (LOWESS + Robust Z-score)'
        else:
            title = f'{param_name} Grid ({outlier_count} outliers)\nMethods: {methods_str}'
        
        ax.set_title(title, fontweight='bold', fontsize=11)
        
        cbar = plt.colorbar(im, ax=ax, shrink=0.8)
        cbar.set_label('Outlier Detected', rotation=270, labelpad=20)
    
    plt.tight_layout()
    plt.subplots_adjust(top=0.88, bottom=0.08)
    figs.append(fig2)
    
    return figs

def process_all_files_enhanced_optimized_with_tuning(base_dir, output_suffix="_method_15_inference", 
                                                   batch_size=50, pdf_limit=50, run_tuning=True):
    """
    Enhanced processing with state parameter hyperparameter tuning and alpha-style location processing.
    
    This function adds hyperparameter tuning capability to optimize state parameter
    smoothing for maximum concordance correlation coefficient, while making location
    processing identical to alpha processing.
    """
    print("🚀 ENHANCED PROCESSING: ALPHA-STYLE LOCATION + OPTIMIZED STATE")
    print("="*80)
    print("🔄 LOCATION: Now identical to Alpha (LOWESS + Robust Z-score + EMA)")
    print("🎯 STATE: Hyperparameter tuned for maximum concordance correlation")
    print("✅ ALPHA/BETA: Optimized balanced approach (unchanged)")
    print("📊 PROCESSING: PDF reports for first 50 files only")
    print("="*80)
    
    detector = AdvancedOutlierDetector(save_reports=True)
    
    # Verify location configuration matches alpha
    alpha_config = detector.method_config['alpha']
    location_config = detector.method_config['location']
    
    print(f"\n🔍 CONFIGURATION VERIFICATION:")
    print(f"Alpha methods: {alpha_config['methods']}")
    print(f"Location methods: {location_config['methods']} ✓ MATCH")
    print(f"Alpha smoothing: {alpha_config['smoothing_method']}")
    print(f"Location smoothing: {location_config['smoothing_method']} ✓ MATCH")
    print(f"Alpha final smoothing: {alpha_config['final_smoothing_method']}")
    print(f"Location final smoothing: {location_config['final_smoothing_method']} ✓ MATCH")
    
    # Step 1: Run hyperparameter tuning if requested
    if run_tuning:
        print("\n🔍 STEP 1: HYPERPARAMETER TUNING FOR STATE")
        print("This will find optimal state parameter settings for your data...")
        
        best_config, best_score = detector.tune_state_parameters(
            base_dir=base_dir,
            num_test_files=8,  # Use 8 files for evaluation
            max_configs=40     # Test up to 40 parameter combinations
        )
        
        # Apply the optimized configuration
        detector.apply_optimized_state_config()
        
        # Save tuning results
        tuning_dir = os.path.join(base_dir, "hyperparameter_tuning_results")
        detector.save_tuning_results(tuning_dir)
        
        print(f"\n✅ Hyperparameter tuning complete! Best score: {best_score:.4f}")
        print("📋 Tuning results saved for analysis")
    else:
        print("\n⏭️ Skipping hyperparameter tuning - using current configuration")
    
    # Step 2: Process all files with optimized (or current) configuration
    print(f"\n🏭 STEP 2: PROCESSING ALL FILES WITH ALPHA-STYLE LOCATION")
    
    # Setup directories
    ideal_dir = os.path.join(base_dir, "ideal_inference")
    realistic_dir = os.path.join(base_dir, "realistic_inference")
    output_ideal_dir = os.path.join(base_dir, f"ideal_inference{output_suffix}")
    output_realistic_dir = os.path.join(base_dir, f"realistic_inference{output_suffix}")
    plots_dir = os.path.join(base_dir, f"{output_suffix}_15_analysis_inference")
    reports_dir = os.path.join(base_dir, f"{output_suffix}_15_reports_inference")
    
    # Create directories
    for dir_path in [output_ideal_dir, output_realistic_dir, plots_dir, reports_dir]:
        if not os.path.exists(dir_path):
            os.makedirs(dir_path)
            print(f"Created directory: {dir_path}")
    
    # Get common files
    ideal_files = set(os.listdir(ideal_dir))
    realistic_files = set(os.listdir(realistic_dir))
    common_files = sorted(list(ideal_files.intersection(realistic_files)))
    
    print(f"Found {len(common_files)} files to process")
    print(f"📊 Will generate PDF reports for first {pdf_limit} files")
    print(f"📁 All files will be processed and saved")
    
    files_processed = 0
    pdf_files_processed = 0
    
    # Create single PDF for the first batch
    pdf_path = os.path.join(plots_dir, f"smoothing_15_inference.pdf")
    
    with PdfPages(pdf_path) as pdf:
        for file_idx, filename in enumerate(tqdm(common_files, desc="Processing files")):
            
            ideal_file_path = os.path.join(ideal_dir, filename)
            realistic_file_path = os.path.join(realistic_dir, filename)
            output_ideal_path = os.path.join(output_ideal_dir, filename)
            output_realistic_path = os.path.join(output_realistic_dir, filename)
            
            try:
                # Process file pair with optimized configuration
                result = process_file_pair_advanced(
                    ideal_file_path, realistic_file_path, detector
                )
                
                # Save processed data
                with open(output_ideal_path, 'w') as f:
                    for value in result['processed_ideal']:
                        f.write(f"{value:.10f}\n")
                
                with open(output_realistic_path, 'w') as f:
                    for value in result['processed_realistic']:
                        f.write(f"{value:.10f}\n")
                
                files_processed += 1
                detector.global_stats['files_processed'] += 1
                
                # Generate plots only for first PDF_LIMIT files
                if file_idx < pdf_limit:
                    figs = create_advanced_outlier_plot(
                        filename, 
                        result['original_ideal'], result['original_realistic'],
                        result['processed_ideal'], result['processed_realistic'], 
                        result['file_report']
                    )
                    
                    # Save figures to PDF
                    for fig in figs:
                        pdf.savefig(fig, bbox_inches='tight', dpi=150)
                        plt.close(fig)
                    
                    pdf_files_processed += 1
                
                # Save detailed report for all files
                if detector.save_reports:
                    report_file = os.path.join(reports_dir, f"{filename.replace('.txt', '_alpha_style_location_report.json')}")
                    serializable_report = ensure_json_serializable(result['file_report'])
                    with open(report_file, 'w') as f:
                        json.dump(serializable_report, f, indent=2)
            
            except Exception as e:
                print(f"Error processing {filename}: {e}")
                continue
    
    # Save global statistics
    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
    stats_file = os.path.join(reports_dir, f"alpha_style_location_global_stats_{timestamp}.json")
    with open(stats_file, 'w') as f:
        json.dump(ensure_json_serializable(detector.global_stats), f, indent=2)
    
    print(f"\n🎉 ALPHA-STYLE LOCATION PROCESSING COMPLETE!")
    print(f"📁 Output directories:")
    print(f"  Ideal: {output_ideal_dir}")
    print(f"  Realistic: {output_realistic_dir}")
    print(f"📊 Analysis: {plots_dir}")
    print(f"📋 Reports: {reports_dir}")
    print(f"Files processed: {files_processed}")
    print(f"PDF reports generated: {pdf_files_processed}")
    
    print(f"\n🔄 LOCATION PARAMETER NOW USES:")
    print(f"  ✅ Same methods as Alpha: {location_config['methods']}")
    print(f"  ✅ Same smoothing as Alpha: {location_config['smoothing_method']}")
    print(f"  ✅ Same final smoothing as Alpha: {location_config['final_smoothing_method']}")
    print(f"  ✅ Same outlier percentage limit: {location_config['max_outlier_percentage']}%")
    
    if run_tuning and detector.best_state_config:
        print(f"\n🎯 OPTIMIZED STATE CONFIGURATION:")
        print(f"  Gaussian σ: {detector.best_state_config['gaussian_sigma_moderate']}")
        print(f"  Spline factor: {detector.best_state_config['spline_smoothing_factor']}")
        print(f"  EMA α: {detector.best_state_config['ema_alpha']}")
        print(f"  Final σ: {detector.best_state_config['final_gaussian_sigma']}")
        print(f"  LOF contamination: {detector.best_state_config['lof_contamination']}")
    
    return output_ideal_dir, output_realistic_dir, detector

# Main execution with tuning option
if __name__ == "__main__":
    base_dir = r"D:\allen\DLR\Original_data\inference_dataset"
    
    print("Choose processing mode:")
    print("1. Full optimization (run tuning + process all files with alpha-style location)")
    print("2. Process only (skip tuning, use current config with alpha-style location)")
    
    choice = input("Enter choice (1 or 2): ").strip()
    
    if choice == "1":
        print("\n🚀 Running full optimization with alpha-style location processing...")
        clean_ideal_dir, clean_realistic_dir, detector = process_all_files_enhanced_optimized_with_tuning(
            base_dir, 
            output_suffix="_method_15_inference",
            batch_size=50,
            pdf_limit=50,
            run_tuning=True
        )
        print(f"\n🎯 Next Steps:")
        print(f"1. Check hyperparameter_tuning_results/ for optimization details")
        print(f"2. Train your Gaussian Process models with {clean_ideal_dir} data")
        print(f"3. Location parameter now processed identically to Alpha!")
        print(f"4. Measure concordance improvement in your ML predictions!")
        
    elif choice == "2":
        print("\nâš¡ Processing with alpha-style location configuration...")
        clean_ideal_dir, clean_realistic_dir, detector = process_all_files_enhanced_optimized_with_tuning(
            base_dir, 
            output_suffix="_method_15_inference",
            batch_size=50,
            pdf_limit=50,
            run_tuning=False
        )
        print(f"\n✅ Location parameter now uses identical processing to Alpha!")
        print(f"💡 Consider running option 1 to also optimize state parameters!")
        
    else:
        print("Invalid choice. Please run again and select 1 or 2.")