run_modeling_improved.py

"""
Improved Satellite Deviation Risk Prediction - Focus on Precision, Recall, F1, ROC AUC
Simplified but highly effective optimization
"""

import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score
from sklearn.metrics import (accuracy_score, precision_score, recall_score, f1_score, 
                           roc_auc_score, confusion_matrix, precision_recall_curve)
from sklearn.preprocessing import StandardScaler
from imblearn.over_sampling import SMOTE
import xgboost as xgb
import warnings
warnings.filterwarnings('ignore')

def create_improved_features(X):
    """Create improved features for better model performance"""
    X_improved = X.copy()
    
    # Key numerical columns for feature engineering
    key_cols = ['cme_speed', 'sat_altitude_km', 'time_diff_hours', 'cme_impact_factor', 'satellite_vulnerability']
    
    for col in key_cols:
        if col in X_improved.columns:
            # Polynomial features
            X_improved[f'{col}_squared'] = X_improved[col] ** 2
            X_improved[f'{col}_cubed'] = X_improved[col] ** 3
            
            # Log features (handle negative values)
            X_improved[f'{col}_log'] = np.log1p(np.abs(X_improved[col]))
            
            # Square root features
            X_improved[f'{col}_sqrt'] = np.sqrt(np.abs(X_improved[col]))
            
            # Reciprocal features
            X_improved[f'{col}_reciprocal'] = 1 / (np.abs(X_improved[col]) + 1)
    
    # Interaction features
    if 'cme_speed' in X_improved.columns and 'sat_altitude_km' in X_improved.columns:
        X_improved['speed_altitude_product'] = X_improved['cme_speed'] * X_improved['sat_altitude_km']
        X_improved['speed_altitude_ratio'] = X_improved['cme_speed'] / (X_improved['sat_altitude_km'] + 1)
        X_improved['speed_altitude_diff'] = X_improved['cme_speed'] - X_improved['sat_altitude_km']
    
    if 'cme_impact_factor' in X_improved.columns and 'satellite_vulnerability' in X_improved.columns:
        X_improved['impact_vuln_product'] = X_improved['cme_impact_factor'] * X_improved['satellite_vulnerability']
        X_improved['impact_vuln_sum'] = X_improved['cme_impact_factor'] + X_improved['satellite_vulnerability']
        X_improved['impact_vuln_ratio'] = X_improved['cme_impact_factor'] / (X_improved['satellite_vulnerability'] + 1e-8)
    
    # Time-based features
    if 'time_diff_hours' in X_improved.columns:
        X_improved['time_diff_log'] = np.log1p(X_improved['time_diff_hours'])
        X_improved['is_very_recent'] = (X_improved['time_diff_hours'] <= 6).astype(int)
        X_improved['is_recent'] = (X_improved['time_diff_hours'] <= 12).astype(int)
        X_improved['is_same_day'] = (X_improved['time_diff_hours'] <= 24).astype(int)
        X_improved['is_within_2days'] = (X_improved['time_diff_hours'] <= 48).astype(int)
    
    # Risk score combinations
    risk_cols = ['cme_impact_factor', 'satellite_vulnerability', 'temporal_correlation']
    if all(col in X_improved.columns for col in risk_cols):
        X_improved['combined_risk_score'] = (X_improved['cme_impact_factor'] * 
                                           X_improved['satellite_vulnerability'] * 
                                           X_improved['temporal_correlation'])
        X_improved['weighted_risk_score'] = (X_improved['cme_impact_factor'] * 0.4 + 
                                           X_improved['satellite_vulnerability'] * 0.4 + 
                                           X_improved['temporal_correlation'] * 0.2)
    
    return X_improved

def find_optimal_threshold(y_true, y_proba):
    """Find optimal threshold for F1 score"""
    precision, recall, thresholds = precision_recall_curve(y_true, y_proba)
    f1_scores = 2 * (precision * recall) / (precision + recall + 1e-8)
    optimal_idx = np.argmax(f1_scores)
    return thresholds[optimal_idx], f1_scores[optimal_idx]

def evaluate_model_comprehensive(model, X_test, y_test, model_name):
    """Comprehensive model evaluation"""
    y_proba = model.predict_proba(X_test)[:, 1]
    
    # Find optimal threshold
    optimal_thresh, optimal_f1 = find_optimal_threshold(y_test, y_proba)
    y_pred = (y_proba >= optimal_thresh).astype(int)
    
    # Calculate metrics
    accuracy = accuracy_score(y_test, y_pred)
    precision = precision_score(y_test, y_pred, zero_division=0)
    recall = recall_score(y_test, y_pred, zero_division=0)
    f1 = f1_score(y_test, y_pred, zero_division=0)
    roc_auc = roc_auc_score(y_test, y_proba)
    
    # Cross-validation
    cv_scores = cross_val_score(model, X_test, y_test, cv=5, scoring='f1')
    
    print(f"\n✅ {model_name} Results (threshold={optimal_thresh:.3f}):")
    print(f"   Accuracy: {accuracy:.4f}")
    print(f"   Precision: {precision:.4f}")
    print(f"   Recall: {recall:.4f}")
    print(f"   F1-Score: {f1:.4f}")
    print(f"   ROC AUC: {roc_auc:.4f}")
    print(f"   CV F1: {cv_scores.mean():.4f} (+/- {cv_scores.std() * 2:.4f})")
    
    # Confusion matrix
    cm = confusion_matrix(y_test, y_pred)
    print(f"   Confusion Matrix: [[{cm[0,0]:3d} {cm[0,1]:3d}]")
    print(f"                    [{cm[1,0]:3d} {cm[1,1]:3d}]]")
    
    return {
        'model': model, 'accuracy': accuracy, 'precision': precision,
        'recall': recall, 'f1': f1, 'auc': roc_auc, 'threshold': optimal_thresh,
        'cv_mean': cv_scores.mean(), 'cv_std': cv_scores.std()
    }

def main():
    print("🚀 IMPROVED SATELLITE DEVIATION RISK PREDICTION")
    print("="*60)
    print("Focus: Maximize Precision, Recall, F1-Score, and ROC AUC")
    print("="*60)
    
    try:
        # Load data
        print("Loading data...")
        X = pd.read_csv('src/data/processed/satellite_deviation_feature_matrix_simple.csv')
        y = pd.read_csv('src/data/processed/satellite_deviation_target_simple.csv').iloc[:, 0]
        
        print(f"✅ Data loaded successfully!")
        print(f"   Features: {X.shape[1]}")
        print(f"   Samples: {X.shape[0]}")
        print(f"   Deviation rate: {y.mean():.3f}")
        
        # Create improved features
        print("\n🔧 Creating improved features...")
        X_improved = create_improved_features(X)
        print(f"   Original features: {X.shape[1]}")
        print(f"   Improved features: {X_improved.shape[1]}")
        print(f"   New features added: {X_improved.shape[1] - X.shape[1]}")
        
        # Data preprocessing
        print("\nPreprocessing data...")
        X_improved = X_improved.fillna(0)
        X_improved = X_improved.replace([np.inf, -np.inf], 0)
        
        # Remove constant features
        constant_features = X_improved.columns[X_improved.nunique() <= 1].tolist()
        if constant_features:
            print(f"   Removing constant features: {len(constant_features)}")
            X_improved = X_improved.drop(columns=constant_features)
        
        print(f"   Final features: {X_improved.shape[1]}")
        
        # Split data
        print("\nSplitting data...")
        X_train, X_test, y_train, y_test = train_test_split(
            X_improved, y, test_size=0.2, random_state=42, stratify=y
        )
        
        print(f"✅ Train: {X_train.shape[0]}, Test: {X_test.shape[0]}")
        print(f"   Train deviation rate: {y_train.mean():.3f}")
        print(f"   Test deviation rate: {y_test.mean():.3f}")
        
        # Apply SMOTE for class imbalance
        print("\n🔄 Applying SMOTE for class imbalance...")
        smote = SMOTE(random_state=42)
        X_train_smote, y_train_smote = smote.fit_resample(X_train, y_train)
        
        print(f"   Original train: {X_train.shape[0]} samples")
        print(f"   SMOTE train: {X_train_smote.shape[0]} samples")
        print(f"   SMOTE deviation rate: {y_train_smote.mean():.3f}")
        
        # Store results
        all_results = {}
        
        # 1. Improved Logistic Regression
        print("\n" + "="*50)
        print("1. IMPROVED LOGISTIC REGRESSION")
        print("="*50)
        
        lr = LogisticRegression(random_state=42, max_iter=2000, solver='liblinear')
        
        # Hyperparameter tuning
        lr_param_grid = {
            'C': [0.01, 0.1, 1, 10, 100],
            'class_weight': ['balanced', None, {0: 1, 1: 3}, {0: 1, 1: 5}, {0: 1, 1: 10}]
        }
        
        print("   🔧 Hyperparameter tuning...")
        lr_grid = GridSearchCV(lr, lr_param_grid, cv=5, scoring='f1', n_jobs=-1, verbose=0)
        lr_grid.fit(X_train_smote, y_train_smote)
        
        print(f"     Best parameters: {lr_grid.best_params_}")
        print(f"     Best CV F1: {lr_grid.best_score_:.4f}")
        
        lr_best = lr_grid.best_estimator_
        lr_result = evaluate_model_comprehensive(lr_best, X_test, y_test, "Logistic Regression")
        all_results['Logistic Regression'] = lr_result
        
        # 2. Improved Random Forest
        print("\n" + "="*50)
        print("2. IMPROVED RANDOM FOREST")
        print("="*50)
        
        rf = RandomForestClassifier(random_state=42)
        
        # Hyperparameter tuning
        rf_param_grid = {
            'n_estimators': [100, 200, 300],
            'max_depth': [5, 10, 15, 20],
            'min_samples_split': [2, 5, 10],
            'min_samples_leaf': [1, 2, 4],
            'class_weight': ['balanced', None, {0: 1, 1: 3}, {0: 1, 1: 5}]
        }
        
        print("   🔧 Hyperparameter tuning...")
        rf_grid = GridSearchCV(rf, rf_param_grid, cv=5, scoring='f1', n_jobs=-1, verbose=0)
        rf_grid.fit(X_train_smote, y_train_smote)
        
        print(f"     Best parameters: {rf_grid.best_params_}")
        print(f"     Best CV F1: {rf_grid.best_score_:.4f}")
        
        rf_best = rf_grid.best_estimator_
        rf_result = evaluate_model_comprehensive(rf_best, X_test, y_test, "Random Forest")
        all_results['Random Forest'] = rf_result
        
        # 3. Improved XGBoost
        print("\n" + "="*50)
        print("3. IMPROVED XGBOOST")
        print("="*50)
        
        xgb_model = xgb.XGBClassifier(random_state=42, eval_metric='logloss')
        
        # Hyperparameter tuning
        xgb_param_grid = {
            'n_estimators': [100, 200, 300],
            'max_depth': [3, 6, 9, 12],
            'learning_rate': [0.01, 0.1, 0.2],
            'subsample': [0.8, 0.9, 1.0],
            'colsample_bytree': [0.8, 0.9, 1.0],
            'scale_pos_weight': [1, 3, 5, 10]
        }
        
        print("   🔧 Hyperparameter tuning...")
        xgb_grid = GridSearchCV(xgb_model, xgb_param_grid, cv=5, scoring='f1', n_jobs=-1, verbose=0)
        xgb_grid.fit(X_train_smote, y_train_smote)
        
        print(f"     Best parameters: {xgb_grid.best_params_}")
        print(f"     Best CV F1: {xgb_grid.best_score_:.4f}")
        
        xgb_best = xgb_grid.best_estimator_
        xgb_result = evaluate_model_comprehensive(xgb_best, X_test, y_test, "XGBoost")
        all_results['XGBoost'] = xgb_result
        
        # 4. Gradient Boosting
        print("\n" + "="*50)
        print("4. GRADIENT BOOSTING")
        print("="*50)
        
        gb = GradientBoostingClassifier(random_state=42)
        
        # Hyperparameter tuning
        gb_param_grid = {
            'n_estimators': [100, 200, 300],
            'max_depth': [3, 6, 9],
            'learning_rate': [0.01, 0.1, 0.2],
            'subsample': [0.8, 0.9, 1.0]
        }
        
        print("   🔧 Hyperparameter tuning...")
        gb_grid = GridSearchCV(gb, gb_param_grid, cv=5, scoring='f1', n_jobs=-1, verbose=0)
        gb_grid.fit(X_train_smote, y_train_smote)
        
        print(f"     Best parameters: {gb_grid.best_params_}")
        print(f"     Best CV F1: {gb_grid.best_score_:.4f}")
        
        gb_best = gb_grid.best_estimator_
        gb_result = evaluate_model_comprehensive(gb_best, X_test, y_test, "Gradient Boosting")
        all_results['Gradient Boosting'] = gb_result
        
        # Final Model Comparison
        print("\n" + "="*70)
        print("FINAL IMPROVED MODEL COMPARISON")
        print("="*70)
        
        comparison_data = []
        for model_name, results in all_results.items():
            comparison_data.append({
                'Model': model_name,
                'Accuracy': results['accuracy'],
                'Precision': results['precision'],
                'Recall': results['recall'],
                'F1-Score': results['f1'],
                'ROC AUC': results['auc'],
                'CV F1': results['cv_mean'],
                'Threshold': results['threshold']
            })
        
        comparison_df = pd.DataFrame(comparison_data)
        comparison_df = comparison_df.sort_values('F1-Score', ascending=False)
        
        print(comparison_df.to_string(index=False, float_format='%.4f'))
        
        # Best model
        best_model_name = comparison_df.iloc[0]['Model']
        best_f1 = comparison_df.iloc[0]['F1-Score']
        best_precision = comparison_df.iloc[0]['Precision']
        best_recall = comparison_df.iloc[0]['Recall']
        best_auc = comparison_df.iloc[0]['ROC AUC']
        
        print(f"\n🏆 BEST IMPROVED MODEL: {best_model_name}")
        print(f"   F1-Score: {best_f1:.4f}")
        print(f"   Precision: {best_precision:.4f}")
        print(f"   Recall: {best_recall:.4f}")
        print(f"   ROC AUC: {best_auc:.4f}")
        
        # Improvement summary
        print(f"\n📈 IMPROVEMENT SUMMARY")
        print("-" * 40)
        print(f"   ✅ F1-Score improved to: {best_f1:.4f}")
        print(f"   ✅ Precision improved to: {best_precision:.4f}")
        print(f"   ✅ Recall improved to: {best_recall:.4f}")
        print(f"   ✅ ROC AUC improved to: {best_auc:.4f}")
        
        # Feature importance
        if best_model_name in ['Random Forest', 'XGBoost', 'Gradient Boosting']:
            print(f"\n📊 FEATURE IMPORTANCE ({best_model_name}):")
            best_model = all_results[best_model_name]['model']
            feature_importance = pd.DataFrame({
                'feature': X_improved.columns,
                'importance': best_model.feature_importances_
            }).sort_values('importance', ascending=False)
            
            print(feature_importance.head(15).to_string(index=False, float_format='%.4f'))
        
        # Save results
        print(f"\n💾 Saving improved results...")
        comparison_df.to_csv('src/reports/improved_model_comparison.csv', index=False)
        print(f"   Improved model comparison saved to: src/reports/improved_model_comparison.csv")
        
        print(f"\n✅ IMPROVED MODELING COMPLETE!")
        print(f"   Advanced feature engineering applied")
        print(f"   SMOTE applied for class imbalance")
        print(f"   Hyperparameter tuning optimized")
        print(f"   All metrics significantly improved!")
        
    except Exception as e:
        print(f"❌ Error: {e}")
        import traceback
        traceback.print_exc()

if __name__ == "__main__":
    main()