refactor: Major reorganization of optimization_engine module structure

BREAKING CHANGE: Module paths have been reorganized for better maintainability.
Backwards compatibility aliases with deprecation warnings are provided.

New Structure:
- core/           - Optimization runners (runner, intelligent_optimizer, etc.)
- processors/     - Data processing
  - surrogates/   - Neural network surrogates
- nx/             - NX/Nastran integration (solver, updater, session_manager)
- study/          - Study management (creator, wizard, state, reset)
- reporting/      - Reports and analysis (visualizer, report_generator)
- config/         - Configuration management (manager, builder)
- utils/          - Utilities (logger, auto_doc, etc.)
- future/         - Research/experimental code

Migration:
- ~200 import changes across 125 files
- All __init__.py files use lazy loading to avoid circular imports
- Backwards compatibility layer supports old import paths with warnings
- All existing functionality preserved

To migrate existing code:
  OLD: from optimization_engine.nx_solver import NXSolver
  NEW: from optimization_engine.nx.solver import NXSolver

  OLD: from optimization_engine.runner import OptimizationRunner
  NEW: from optimization_engine.core.runner import OptimizationRunner

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

optimization_engine/reporting/__init__.py

@@ -0,0 +1,44 @@
+"""
+Reporting & Analysis
+====================
+
+Report generation and results analysis.
+
+Modules:
+- report_generator: HTML/PDF report generation
+- markdown_report: Markdown report format
+- results_analyzer: Comprehensive results analysis
+- visualizer: Plotting and visualization
+- landscape_analyzer: Design space analysis
+"""
+
+# Lazy imports to avoid import errors
+def __getattr__(name):
+    if name == 'generate_optimization_report':
+        from .report_generator import generate_optimization_report
+        return generate_optimization_report
+    elif name == 'generate_markdown_report':
+        from .markdown_report import generate_markdown_report
+        return generate_markdown_report
+    elif name == 'MarkdownReportGenerator':
+        from .markdown_report import MarkdownReportGenerator
+        return MarkdownReportGenerator
+    elif name == 'ResultsAnalyzer':
+        from .results_analyzer import ResultsAnalyzer
+        return ResultsAnalyzer
+    elif name == 'Visualizer':
+        from .visualizer import Visualizer
+        return Visualizer
+    elif name == 'LandscapeAnalyzer':
+        from .landscape_analyzer import LandscapeAnalyzer
+        return LandscapeAnalyzer
+    raise AttributeError(f"module 'optimization_engine.reporting' has no attribute '{name}'")
+
+__all__ = [
+    'generate_optimization_report',
+    'generate_markdown_report',
+    'MarkdownReportGenerator',
+    'ResultsAnalyzer',
+    'Visualizer',
+    'LandscapeAnalyzer',
+]

optimization_engine/reporting/landscape_analyzer.py

@@ -0,0 +1,386 @@
+"""
+Landscape Analyzer - Automatic optimization problem characterization.
+
+This module analyzes the characteristics of an optimization landscape to inform
+intelligent strategy selection. It computes metrics like smoothness, multimodality,
+parameter correlation, and noise level.
+
+Part of Protocol 10: Intelligent Multi-Strategy Optimization (IMSO)
+"""
+
+import numpy as np
+from typing import Dict, List, Optional
+from scipy.stats import spearmanr, variation
+from scipy.spatial.distance import pdist, squareform
+from sklearn.cluster import DBSCAN
+import optuna
+
+
+class LandscapeAnalyzer:
+    """Analyzes optimization landscape characteristics from trial history."""
+
+    def __init__(self, min_trials_for_analysis: int = 10, verbose: bool = True):
+        """
+        Args:
+            min_trials_for_analysis: Minimum trials needed for reliable analysis
+            verbose: Whether to print diagnostic messages
+        """
+        self.min_trials = min_trials_for_analysis
+        self.verbose = verbose
+
+    def analyze(self, study: optuna.Study) -> Dict:
+        """
+        Analyze optimization landscape characteristics.
+
+        STUDY-AWARE: Uses study.trials directly for analysis.
+
+        Args:
+            study: Optuna study with completed trials
+
+        Returns:
+            Dictionary with landscape characteristics:
+            - smoothness: 0-1, how smooth the objective landscape is
+            - multimodal: boolean, multiple local optima detected
+            - n_modes: estimated number of local optima
+            - parameter_correlation: dict of correlation scores
+            - noise_level: estimated noise in evaluations
+            - dimensionality: number of design variables
+            - landscape_type: classification (smooth/rugged/multimodal)
+        """
+        # Get completed trials
+        completed_trials = [t for t in study.trials if t.state == optuna.trial.TrialState.COMPLETE]
+
+        if len(completed_trials) < self.min_trials:
+            return {
+                'ready': False,
+                'total_trials': len(completed_trials),
+                'message': f'Need {self.min_trials - len(completed_trials)} more trials for landscape analysis'
+            }
+
+        # Check if this is a multi-objective study
+        # Multi-objective studies have trial.values (plural), not trial.value
+        is_multi_objective = len(study.directions) > 1
+
+        if is_multi_objective:
+            return {
+                'ready': False,
+                'total_trials': len(completed_trials),
+                'message': 'Landscape analysis not supported for multi-objective optimization'
+            }
+
+        # Extract data
+        X = []  # Parameter values
+        y = []  # Objective values
+        param_names = []
+
+        for trial in completed_trials:
+            X.append(list(trial.params.values()))
+            y.append(trial.value)
+            if not param_names:
+                param_names = list(trial.params.keys())
+
+        X = np.array(X)
+        y = np.array(y)
+
+        # Compute characteristics
+        smoothness = self._compute_smoothness(X, y)
+        multimodal, n_modes = self._detect_multimodality(X, y)
+        correlation_scores = self._compute_parameter_correlation(X, y, param_names)
+        noise_level = self._estimate_noise(X, y)
+        landscape_type = self._classify_landscape(smoothness, multimodal, noise_level, n_modes)
+
+        # Compute parameter ranges for coverage metrics
+        param_ranges = self._compute_parameter_ranges(completed_trials)
+
+        return {
+            'ready': True,
+            'total_trials': len(completed_trials),
+            'dimensionality': X.shape[1],
+            'parameter_names': param_names,
+            'smoothness': smoothness,
+            'multimodal': multimodal,
+            'n_modes': n_modes,
+            'parameter_correlation': correlation_scores,
+            'noise_level': noise_level,
+            'landscape_type': landscape_type,
+            'parameter_ranges': param_ranges,
+            'objective_statistics': {
+                'mean': float(np.mean(y)),
+                'std': float(np.std(y)),
+                'min': float(np.min(y)),
+                'max': float(np.max(y)),
+                'range': float(np.max(y) - np.min(y))
+            }
+        }
+
+    def _compute_smoothness(self, X: np.ndarray, y: np.ndarray) -> float:
+        """
+        Compute landscape smoothness score.
+
+        High smoothness = nearby points have similar objective values
+        Low smoothness = nearby points have very different values (rugged)
+
+        Method: Compare objective differences vs parameter distances
+        """
+        if len(y) < 3:
+            return 0.5  # Unknown
+
+        # Normalize parameters to [0, 1] for fair distance computation
+        X_norm = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0) + 1e-10)
+
+        # Compute pairwise distances in parameter space
+        param_distances = pdist(X_norm, metric='euclidean')
+
+        # Compute pairwise differences in objective space
+        objective_diffs = pdist(y.reshape(-1, 1), metric='euclidean')
+
+        # Smoothness = correlation between distance and objective difference
+        # Smooth landscape: nearby points → similar objectives (high correlation)
+        # Rugged landscape: nearby points → very different objectives (low correlation)
+
+        if len(param_distances) > 0 and len(objective_diffs) > 0:
+            # Filter out zero distances to avoid division issues
+            mask = param_distances > 1e-6
+            if np.sum(mask) > 5:
+                param_distances = param_distances[mask]
+                objective_diffs = objective_diffs[mask]
+
+                # Compute correlation
+                correlation, _ = spearmanr(param_distances, objective_diffs)
+
+                # Convert to smoothness score: high correlation = smooth
+                # Handle NaN from constant arrays
+                if np.isnan(correlation):
+                    smoothness = 0.5
+                else:
+                    smoothness = max(0.0, min(1.0, (correlation + 1.0) / 2.0))
+            else:
+                smoothness = 0.5
+        else:
+            smoothness = 0.5
+
+        return smoothness
+
+    def _detect_multimodality(self, X: np.ndarray, y: np.ndarray) -> tuple:
+        """
+        Detect multiple local optima using clustering.
+
+        Returns:
+            (is_multimodal, n_modes)
+        """
+        if len(y) < 10:
+            return False, 1
+
+        # Find good trials (bottom 30%)
+        threshold = np.percentile(y, 30)
+        good_trials_mask = y <= threshold
+
+        if np.sum(good_trials_mask) < 3:
+            return False, 1
+
+        X_good = X[good_trials_mask]
+
+        # Normalize for clustering
+        X_norm = (X_good - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0) + 1e-10)
+
+        # Use DBSCAN to find clusters of good solutions
+        # If they're spread across multiple regions → multimodal
+        try:
+            clustering = DBSCAN(eps=0.2, min_samples=2).fit(X_norm)
+            n_clusters = len(set(clustering.labels_)) - (1 if -1 in clustering.labels_ else 0)
+
+            is_multimodal = n_clusters > 1
+            n_modes = max(1, n_clusters)
+        except:
+            is_multimodal = False
+            n_modes = 1
+
+        return is_multimodal, n_modes
+
+    def _compute_parameter_correlation(self, X: np.ndarray, y: np.ndarray, param_names: List[str]) -> Dict:
+        """
+        Compute correlation between each parameter and objective.
+
+        Returns dict of {param_name: correlation_score}
+        High absolute correlation → parameter strongly affects objective
+        """
+        correlations = {}
+
+        for i, param_name in enumerate(param_names):
+            param_values = X[:, i]
+
+            # Spearman correlation (handles nonlinearity)
+            corr, p_value = spearmanr(param_values, y)
+
+            if np.isnan(corr):
+                corr = 0.0
+
+            correlations[param_name] = {
+                'correlation': float(corr),
+                'abs_correlation': float(abs(corr)),
+                'p_value': float(p_value) if not np.isnan(p_value) else 1.0
+            }
+
+        # Compute overall correlation strength
+        avg_abs_corr = np.mean([v['abs_correlation'] for v in correlations.values()])
+
+        correlations['overall_strength'] = float(avg_abs_corr)
+
+        return correlations
+
+    def _estimate_noise(self, X: np.ndarray, y: np.ndarray) -> float:
+        """
+        Estimate noise level in objective evaluations.
+
+        For deterministic FEA simulations, this should be very low.
+        High noise would suggest numerical issues or simulation instability.
+
+        Method: Look at local variations - similar inputs should give similar outputs.
+        Wide exploration range (high CV) is NOT noise.
+        """
+        if len(y) < 10:
+            return 0.0
+
+        # Calculate pairwise distances in parameter space
+        from scipy.spatial.distance import pdist, squareform
+
+        # Normalize X to [0,1] for distance calculation
+        X_norm = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0) + 1e-10)
+
+        # Compute pairwise distances
+        param_distances = squareform(pdist(X_norm, 'euclidean'))
+        objective_diffs = np.abs(y[:, np.newaxis] - y[np.newaxis, :])
+
+        # Find pairs that are close in parameter space (distance < 0.1)
+        close_pairs_mask = (param_distances > 1e-6) & (param_distances < 0.1)
+
+        if np.sum(close_pairs_mask) < 5:
+            # Not enough close pairs to assess noise
+            return 0.0
+
+        # For close pairs, measure objective variation
+        # True noise: close inputs give very different outputs
+        # Smooth function: close inputs give similar outputs
+        close_objective_diffs = objective_diffs[close_pairs_mask]
+        close_param_dists = param_distances[close_pairs_mask]
+
+        # Normalize by expected difference based on smoothness
+        # Noise = unexpected variation for nearby points
+        expected_diff = np.median(close_objective_diffs / (close_param_dists + 1e-10))
+        actual_std = np.std(close_objective_diffs / (close_param_dists + 1e-10))
+
+        # Coefficient of variation of local gradients
+        if expected_diff > 1e-6:
+            local_cv = actual_std / expected_diff
+            noise_score = min(1.0, local_cv / 2.0)
+        else:
+            noise_score = 0.0
+
+        return float(noise_score)
+
+    def _classify_landscape(self, smoothness: float, multimodal: bool, noise: float, n_modes: int = 1) -> str:
+        """
+        Classify landscape type for strategy selection.
+
+        Args:
+            smoothness: Smoothness score (0-1)
+            multimodal: Whether multiple modes detected
+            noise: Noise level (0-1)
+            n_modes: Number of modes detected
+
+        Returns one of:
+        - 'smooth_unimodal': Single smooth bowl (best for CMA-ES, GP-BO)
+        - 'smooth_multimodal': Multiple smooth regions (good for GP-BO, TPE)
+        - 'rugged_unimodal': Single rugged region (TPE, hybrid)
+        - 'rugged_multimodal': Multiple rugged regions (TPE, evolutionary)
+        - 'noisy': High noise level (robust methods)
+        """
+        # IMPROVEMENT: Detect false multimodality from smooth continuous manifolds
+        # If only 2 modes detected with high smoothness and low noise,
+        # it's likely a continuous smooth surface, not true multimodality
+        if multimodal and n_modes == 2 and smoothness > 0.6 and noise < 0.2:
+            if self.verbose:
+                print(f"[LANDSCAPE] Reclassifying: 2 modes with smoothness={smoothness:.2f}, noise={noise:.2f}")
+                print(f"[LANDSCAPE] This appears to be a smooth continuous manifold, not true multimodality")
+            multimodal = False  # Override: treat as unimodal
+
+        if noise > 0.5:
+            return 'noisy'
+
+        if smoothness > 0.6:
+            if multimodal:
+                return 'smooth_multimodal'
+            else:
+                return 'smooth_unimodal'
+        else:
+            if multimodal:
+                return 'rugged_multimodal'
+            else:
+                return 'rugged_unimodal'
+
+    def _compute_parameter_ranges(self, trials: List) -> Dict:
+        """Compute explored parameter ranges."""
+        if not trials:
+            return {}
+
+        param_names = list(trials[0].params.keys())
+        ranges = {}
+
+        for param in param_names:
+            values = [t.params[param] for t in trials]
+            distribution = trials[0].distributions[param]
+
+            ranges[param] = {
+                'explored_min': float(np.min(values)),
+                'explored_max': float(np.max(values)),
+                'explored_range': float(np.max(values) - np.min(values)),
+                'bounds_min': float(distribution.low),
+                'bounds_max': float(distribution.high),
+                'bounds_range': float(distribution.high - distribution.low),
+                'coverage': float((np.max(values) - np.min(values)) / (distribution.high - distribution.low))
+            }
+
+        return ranges
+
+
+def print_landscape_report(landscape: Dict, verbose: bool = True):
+    """Print formatted landscape analysis report."""
+    if not verbose:
+        return
+
+    # Handle None (multi-objective studies)
+    if landscape is None:
+        print(f"\n  [LANDSCAPE ANALYSIS] Skipped for multi-objective optimization")
+        return
+
+    if not landscape.get('ready', False):
+        print(f"\n  [LANDSCAPE ANALYSIS] {landscape.get('message', 'Not ready')}")
+        return
+
+    print(f"\n{'='*70}")
+    print(f"  LANDSCAPE ANALYSIS REPORT")
+    print(f"{'='*70}")
+    print(f"  Total Trials Analyzed: {landscape['total_trials']}")
+    print(f"  Dimensionality: {landscape['dimensionality']} parameters")
+    print(f"\n  LANDSCAPE CHARACTERISTICS:")
+    print(f"    Type: {landscape['landscape_type'].upper()}")
+    print(f"    Smoothness: {landscape['smoothness']:.2f} {'(smooth)' if landscape['smoothness'] > 0.6 else '(rugged)'}")
+    print(f"    Multimodal: {'YES' if landscape['multimodal'] else 'NO'} ({landscape['n_modes']} modes)")
+    print(f"    Noise Level: {landscape['noise_level']:.2f} {'(low)' if landscape['noise_level'] < 0.3 else '(high)'}")
+
+    print(f"\n  PARAMETER CORRELATIONS:")
+    for param, info in landscape['parameter_correlation'].items():
+        if param != 'overall_strength':
+            corr = info['correlation']
+            strength = 'strong' if abs(corr) > 0.5 else 'moderate' if abs(corr) > 0.3 else 'weak'
+            direction = 'positive' if corr > 0 else 'negative'
+            print(f"    {param}: {corr:+.3f} ({strength} {direction})")
+
+    print(f"\n  OBJECTIVE STATISTICS:")
+    stats = landscape['objective_statistics']
+    print(f"    Best: {stats['min']:.6f}")
+    print(f"    Mean: {stats['mean']:.6f}")
+    print(f"    Std:  {stats['std']:.6f}")
+    print(f"    Range: {stats['range']:.6f}")
+
+    print(f"{'='*70}\n")

optimization_engine/reporting/markdown_report.py

@@ -0,0 +1,569 @@
+"""
+Generate comprehensive markdown optimization reports with graphs.
+Uses Optuna's built-in visualization library for professional-quality plots.
+"""
+
+import json
+import sys
+from pathlib import Path
+from typing import Dict, Any, List, Optional
+import numpy as np
+import matplotlib
+matplotlib.use('Agg')  # Non-interactive backend
+import matplotlib.pyplot as plt
+import optuna
+from optuna.visualization import (
+    plot_optimization_history,
+    plot_parallel_coordinate,
+    plot_param_importances,
+    plot_slice,
+    plot_contour
+)
+
+
+def create_confidence_progression_plot(confidence_history: List[Dict], phase_transitions: List[Dict], output_dir: Path) -> Optional[str]:
+    """Create confidence progression plot showing confidence metrics over trials."""
+    if not confidence_history:
+        return None
+
+    trial_numbers = [c['trial_number'] for c in confidence_history]
+    overall = [c['confidence_metrics']['overall_confidence'] for c in confidence_history]
+    convergence = [c['confidence_metrics']['convergence_score'] for c in confidence_history]
+    coverage = [c['confidence_metrics']['exploration_coverage'] for c in confidence_history]
+    stability = [c['confidence_metrics']['prediction_stability'] for c in confidence_history]
+
+    plt.figure(figsize=(12, 7))
+    plt.plot(trial_numbers, overall, 'b-', linewidth=2.5, label='Overall Confidence')
+    plt.plot(trial_numbers, convergence, 'g--', alpha=0.7, label='Convergence Score')
+    plt.plot(trial_numbers, coverage, 'orange', linestyle='--', alpha=0.7, label='Exploration Coverage')
+    plt.plot(trial_numbers, stability, 'purple', linestyle='--', alpha=0.7, label='Prediction Stability')
+
+    # Mark phase transitions
+    for transition in phase_transitions:
+        trial_num = transition['trial_number']
+        plt.axvline(x=trial_num, color='red', linestyle='-', linewidth=2, alpha=0.8)
+        plt.text(trial_num, 0.95, f'  Exploitation Phase', rotation=90,
+                verticalalignment='top', fontsize=10, color='red', fontweight='bold')
+
+    # Mark confidence threshold
+    plt.axhline(y=0.65, color='gray', linestyle=':', linewidth=1.5, alpha=0.6, label='Confidence Threshold (65%)')
+
+    plt.xlabel('Trial Number', fontsize=11)
+    plt.ylabel('Confidence Score (0-1)', fontsize=11)
+    plt.title('Surrogate Confidence Progression', fontsize=13, fontweight='bold')
+    plt.legend(loc='lower right', fontsize=9)
+    plt.grid(True, alpha=0.3)
+    plt.ylim(0, 1.05)
+    plt.tight_layout()
+
+    plot_file = output_dir / 'confidence_progression.png'
+    plt.savefig(plot_file, dpi=150)
+    plt.close()
+
+    return plot_file.name
+
+
+def create_convergence_plot(history: List[Dict], target: Optional[float], output_dir: Path) -> str:
+    """Create convergence plot showing best objective over trials."""
+    trial_numbers = [t['trial_number'] for t in history]
+    objectives = [t['objective'] for t in history]
+
+    # Calculate cumulative best
+    cumulative_best = []
+    current_best = float('inf')
+    for obj in objectives:
+        current_best = min(current_best, obj)
+        cumulative_best.append(current_best)
+
+    plt.figure(figsize=(10, 6))
+    plt.plot(trial_numbers, objectives, 'o-', alpha=0.5, label='Trial objective')
+    plt.plot(trial_numbers, cumulative_best, 'r-', linewidth=2, label='Best so far')
+
+    if target is not None:
+        plt.axhline(y=0, color='g', linestyle='--', linewidth=2, label=f'Target (error = 0)')
+
+    plt.xlabel('Trial Number')
+    plt.ylabel('Objective Value (Error from Target)')
+    plt.title('Optimization Convergence')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.tight_layout()
+
+    plot_file = output_dir / 'convergence_plot.png'
+    plt.savefig(plot_file, dpi=150)
+    plt.close()
+
+    return plot_file.name
+
+
+def create_design_space_plot(history: List[Dict], output_dir: Path) -> str:
+    """Create 2D design space exploration plot."""
+    first_trial = history[0]
+    var_names = list(first_trial['design_variables'].keys())
+
+    if len(var_names) != 2:
+        return None  # Only works for 2D problems
+
+    var1_name, var2_name = var_names
+    var1_values = [t['design_variables'][var1_name] for t in history]
+    var2_values = [t['design_variables'][var2_name] for t in history]
+    objectives = [t['objective'] for t in history]
+
+    plt.figure(figsize=(10, 8))
+    scatter = plt.scatter(var1_values, var2_values, c=objectives, s=100,
+                          cmap='viridis', alpha=0.6, edgecolors='black')
+    plt.colorbar(scatter, label='Objective Value')
+    plt.xlabel(var1_name.replace('_', ' ').title())
+    plt.ylabel(var2_name.replace('_', ' ').title())
+    plt.title('Design Space Exploration')
+    plt.grid(True, alpha=0.3)
+    plt.tight_layout()
+
+    plot_file = output_dir / 'design_space_plot.png'
+    plt.savefig(plot_file, dpi=150)
+    plt.close()
+
+    return plot_file.name
+
+
+def create_parameter_sensitivity_plot(history: List[Dict], output_dir: Path) -> str:
+    """Create parameter sensitivity plots."""
+    first_trial = history[0]
+    var_names = list(first_trial['design_variables'].keys())
+
+    fig, axes = plt.subplots(1, len(var_names), figsize=(6*len(var_names), 5))
+    if len(var_names) == 1:
+        axes = [axes]
+
+    for idx, var_name in enumerate(var_names):
+        var_values = [t['design_variables'][var_name] for t in history]
+        objectives = [t['objective'] for t in history]
+
+        axes[idx].scatter(var_values, objectives, alpha=0.6, s=50)
+        axes[idx].set_xlabel(var_name.replace('_', ' ').title())
+        axes[idx].set_ylabel('Objective Value')
+        axes[idx].set_title(f'Sensitivity to {var_name.replace("_", " ").title()}')
+        axes[idx].grid(True, alpha=0.3)
+
+    plt.tight_layout()
+    plot_file = output_dir / 'parameter_sensitivity.png'
+    plt.savefig(plot_file, dpi=150)
+    plt.close()
+
+    return plot_file.name
+
+
+def create_optuna_plots(study: optuna.Study, output_dir: Path) -> Dict[str, str]:
+    """
+    Create professional Optuna visualization plots.
+
+    Args:
+        study: Optuna study object
+        output_dir: Directory to save plots
+
+    Returns:
+        Dictionary mapping plot names to filenames
+    """
+    plots = {}
+
+    try:
+        # 1. Parallel Coordinate Plot - shows parameter interactions
+        fig = plot_parallel_coordinate(study)
+        if fig is not None:
+            plot_file = output_dir / 'optuna_parallel_coordinate.png'
+            fig.write_image(str(plot_file), width=1200, height=600)
+            plots['parallel_coordinate'] = plot_file.name
+    except Exception as e:
+        print(f"Warning: Could not create parallel coordinate plot: {e}")
+
+    try:
+        # 2. Optimization History - convergence over trials
+        fig = plot_optimization_history(study)
+        if fig is not None:
+            plot_file = output_dir / 'optuna_optimization_history.png'
+            fig.write_image(str(plot_file), width=1000, height=600)
+            plots['optimization_history'] = plot_file.name
+    except Exception as e:
+        print(f"Warning: Could not create optimization history plot: {e}")
+
+    try:
+        # 3. Parameter Importances - which parameters matter most
+        fig = plot_param_importances(study)
+        if fig is not None:
+            plot_file = output_dir / 'optuna_param_importances.png'
+            fig.write_image(str(plot_file), width=800, height=500)
+            plots['param_importances'] = plot_file.name
+    except Exception as e:
+        print(f"Warning: Could not create parameter importance plot: {e}")
+
+    try:
+        # 4. Slice Plot - individual parameter effects
+        fig = plot_slice(study)
+        if fig is not None:
+            plot_file = output_dir / 'optuna_slice.png'
+            fig.write_image(str(plot_file), width=1000, height=600)
+            plots['slice'] = plot_file.name
+    except Exception as e:
+        print(f"Warning: Could not create slice plot: {e}")
+
+    try:
+        # 5. Contour Plot - parameter interaction heatmap (2D only)
+        if len(study.best_params) == 2:
+            fig = plot_contour(study)
+            if fig is not None:
+                plot_file = output_dir / 'optuna_contour.png'
+                fig.write_image(str(plot_file), width=800, height=800)
+                plots['contour'] = plot_file.name
+    except Exception as e:
+        print(f"Warning: Could not create contour plot: {e}")
+
+    return plots
+
+
+def generate_markdown_report(history_file: Path, target_value: Optional[float] = None,
+                             tolerance: float = 0.1, reports_dir: Optional[Path] = None,
+                             study: Optional[optuna.Study] = None) -> str:
+    """Generate comprehensive markdown optimization report with graphs."""
+
+    # Load history
+    with open(history_file) as f:
+        history = json.load(f)
+
+    if not history:
+        return "# Optimization Report\n\nNo optimization history found."
+
+    # Graphs should be saved to 3_reports/ folder (same as markdown file)
+    study_dir = history_file.parent.parent
+    study_name = study_dir.name
+
+    if reports_dir is None:
+        reports_dir = study_dir / "3_reports"
+    reports_dir.mkdir(parents=True, exist_ok=True)
+
+    # Load phase transition and confidence history if available
+    results_dir = study_dir / "2_results"
+    phase_transitions = []
+    confidence_history = []
+    phase_transition_file = results_dir / "phase_transitions.json"
+    confidence_history_file = results_dir / "confidence_history.json"
+
+    if phase_transition_file.exists():
+        try:
+            with open(phase_transition_file) as f:
+                phase_transitions = json.load(f)
+        except Exception:
+            pass
+
+    if confidence_history_file.exists():
+        try:
+            with open(confidence_history_file) as f:
+                confidence_history = json.load(f)
+        except Exception:
+            pass
+
+    # Generate plots in reports folder
+    convergence_plot = create_convergence_plot(history, target_value, reports_dir)
+    design_space_plot = create_design_space_plot(history, reports_dir)
+    sensitivity_plot = create_parameter_sensitivity_plot(history, reports_dir)
+
+    # Generate confidence progression plot if data available
+    confidence_plot = None
+    if confidence_history:
+        print("  Generating confidence progression plot...")
+        confidence_plot = create_confidence_progression_plot(confidence_history, phase_transitions, reports_dir)
+
+    # Generate Optuna plots if study object provided
+    optuna_plots = {}
+    if study is not None:
+        print("  Generating Optuna visualization plots...")
+        optuna_plots = create_optuna_plots(study, reports_dir)
+        print(f"  Generated {len(optuna_plots)} Optuna plots")
+
+    # Build markdown report
+    lines = []
+    lines.append(f"# {study_name.replace('_', ' ').title()} - Optimization Report")
+    lines.append("")
+    lines.append(f"**Total Trials**: {len(history)}")
+    lines.append("")
+
+    # Study information
+    lines.append("## Study Information")
+    lines.append("")
+    first_trial = history[0]
+    design_vars = list(first_trial['design_variables'].keys())
+    lines.append(f"- **Design Variables**: {', '.join([v.replace('_', ' ').title() for v in design_vars])}")
+    lines.append(f"- **Number of Trials**: {len(history)}")
+    lines.append("")
+
+    # Adaptive optimization strategy information
+    if phase_transitions or confidence_history:
+        lines.append("## Adaptive Optimization Strategy")
+        lines.append("")
+        lines.append("This study used adaptive surrogate-based optimization with confidence-driven phase transitions.")
+        lines.append("")
+
+        if phase_transitions:
+            lines.append("### Phase Transitions")
+            lines.append("")
+            for transition in phase_transitions:
+                trial_num = transition['trial_number']
+                conf = transition['confidence_metrics']['overall_confidence']
+                lines.append(f"- **Trial #{trial_num}**: EXPLORATION → EXPLOITATION")
+                lines.append(f"  - Confidence at transition: {conf:.1%}")
+                lines.append(f"  - Convergence score: {transition['confidence_metrics']['convergence_score']:.1%}")
+                lines.append(f"  - Exploration coverage: {transition['confidence_metrics']['exploration_coverage']:.1%}")
+                lines.append(f"  - Prediction stability: {transition['confidence_metrics']['prediction_stability']:.1%}")
+                lines.append("")
+        else:
+            lines.append("### Phase Transitions")
+            lines.append("")
+            lines.append("No phase transitions occurred - optimization remained in exploration phase.")
+            lines.append("This may indicate:")
+            lines.append("- Insufficient trials to build surrogate confidence")
+            lines.append("- Poor exploration coverage of the design space")
+            lines.append("- Unstable convergence behavior")
+            lines.append("")
+
+        if confidence_plot:
+            lines.append("### Confidence Progression")
+            lines.append("")
+            lines.append(f"![Confidence Progression]({confidence_plot})")
+            lines.append("")
+            lines.append("This plot shows how the surrogate model confidence evolved over the optimization.")
+            lines.append("The red vertical line (if present) marks the transition to exploitation phase.")
+            lines.append("")
+        lines.append("")
+
+    # Best result
+    objectives = [t['objective'] for t in history]
+    best_idx = np.argmin(objectives)
+    best_trial = history[best_idx]
+
+    lines.append("## Best Result")
+    lines.append("")
+    lines.append(f"- **Trial**: #{best_trial['trial_number']}")
+    lines.append("")
+
+    # Show actual results FIRST (what the client cares about)
+    lines.append("### Achieved Performance")
+    for result, value in best_trial['results'].items():
+        metric_name = result.replace('_', ' ').title()
+        lines.append(f"- **{metric_name}**: {value:.4f}")
+
+        # Show target comparison if available
+        if target_value is not None and 'frequency' in result.lower():
+            error = abs(value - target_value)
+            lines.append(f"  - Target: {target_value:.4f}")
+            lines.append(f"  - Error: {error:.4f} ({(error/target_value*100):.2f}%)")
+    lines.append("")
+
+    # Then design parameters that achieved it
+    lines.append("### Design Parameters")
+    for var, value in best_trial['design_variables'].items():
+        lines.append(f"- **{var.replace('_', ' ').title()}**: {value:.4f}")
+    lines.append("")
+
+    # Technical objective last (for engineers)
+    lines.append("<details>")
+    lines.append("<summary>Technical Details (Objective Function)</summary>")
+    lines.append("")
+    lines.append(f"- **Objective Value (Error)**: {best_trial['objective']:.6f}")
+    lines.append("")
+    lines.append("</details>")
+    lines.append("")
+
+    # Success assessment
+    if target_value is not None:
+        lines.append("## Success Assessment")
+        lines.append("")
+        best_objective = min(objectives)
+
+        if best_objective <= tolerance:
+            lines.append(f"### ✅ TARGET ACHIEVED")
+            lines.append("")
+            lines.append(f"Target value {target_value} was achieved within tolerance {tolerance}!")
+            lines.append(f"- **Best Error**: {best_objective:.6f}")
+        else:
+            lines.append(f"### ⚠️ TARGET NOT YET ACHIEVED")
+            lines.append("")
+            lines.append(f"Target value {target_value} not achieved within tolerance {tolerance}")
+            lines.append(f"- **Best Error**: {best_objective:.6f}")
+            lines.append(f"- **Required Improvement**: {best_objective - tolerance:.6f}")
+            lines.append(f"- **Recommendation**: Continue optimization with more trials")
+        lines.append("")
+
+    # Top 5 trials - show ACTUAL METRICS not just objective
+    lines.append("## Top 5 Trials")
+    lines.append("")
+    sorted_history = sorted(history, key=lambda x: x['objective'])
+
+    # Extract result column names (e.g., "first_frequency")
+    result_cols = list(sorted_history[0]['results'].keys())
+    result_col_names = [r.replace('_', ' ').title() for r in result_cols]
+
+    # Build header with results AND design vars
+    header_cols = ["Rank", "Trial"] + result_col_names + [v.replace('_', ' ').title() for v in design_vars]
+    lines.append("| " + " | ".join(header_cols) + " |")
+    lines.append("|" + "|".join(["-"*max(6, len(c)) for c in header_cols]) + "|")
+
+    for i, trial in enumerate(sorted_history[:5], 1):
+        result_vals = [f"{trial['results'][r]:.2f}" for r in result_cols]
+        var_vals = [f"{trial['design_variables'][v]:.2f}" for v in design_vars]
+        row_data = [str(i), f"#{trial['trial_number']}"] + result_vals + var_vals
+        lines.append("| " + " | ".join(row_data) + " |")
+    lines.append("")
+
+    # Statistics
+    lines.append("## Statistics")
+    lines.append("")
+    lines.append(f"- **Mean Objective**: {np.mean(objectives):.6f}")
+    lines.append(f"- **Std Deviation**: {np.std(objectives):.6f}")
+    lines.append(f"- **Best Objective**: {np.min(objectives):.6f}")
+    lines.append(f"- **Worst Objective**: {np.max(objectives):.6f}")
+    lines.append("")
+
+    # Design variable ranges
+    lines.append("### Design Variable Ranges")
+    lines.append("")
+    for var in design_vars:
+        values = [t['design_variables'][var] for t in history]
+        lines.append(f"**{var.replace('_', ' ').title()}**:")
+        lines.append(f"- Min: {min(values):.6f}")
+        lines.append(f"- Max: {max(values):.6f}")
+        lines.append(f"- Mean: {np.mean(values):.6f}")
+        lines.append("")
+
+    # Convergence plot
+    lines.append("## Convergence Plot")
+    lines.append("")
+    lines.append(f"![Convergence Plot]({convergence_plot})")
+    lines.append("")
+    lines.append("This plot shows how the optimization converged over time. The blue line shows each trial's objective value, while the red line shows the best objective found so far.")
+    lines.append("")
+
+    # Design space plot
+    if design_space_plot:
+        lines.append("## Design Space Exploration")
+        lines.append("")
+        lines.append(f"![Design Space Plot]({design_space_plot})")
+        lines.append("")
+        lines.append("This plot shows which regions of the design space were explored. Darker colors indicate better objective values.")
+        lines.append("")
+
+    # Sensitivity plot
+    lines.append("## Parameter Sensitivity")
+    lines.append("")
+    lines.append(f"![Parameter Sensitivity]({sensitivity_plot})")
+    lines.append("")
+    lines.append("These plots show how each design variable affects the objective value. Steeper slopes indicate higher sensitivity.")
+    lines.append("")
+
+    # Optuna Advanced Visualizations
+    if optuna_plots:
+        lines.append("## Advanced Optimization Analysis (Optuna)")
+        lines.append("")
+        lines.append("The following plots leverage Optuna's professional visualization library to provide deeper insights into the optimization process.")
+        lines.append("")
+
+        # Parallel Coordinate Plot
+        if 'parallel_coordinate' in optuna_plots:
+            lines.append("### Parallel Coordinate Plot")
+            lines.append("")
+            lines.append(f"![Parallel Coordinate]({optuna_plots['parallel_coordinate']})")
+            lines.append("")
+            lines.append("This interactive plot shows how different parameter combinations lead to different objective values. Each line represents one trial, colored by objective value. You can see parameter interactions and identify promising regions.")
+            lines.append("")
+
+        # Optimization History
+        if 'optimization_history' in optuna_plots:
+            lines.append("### Optimization History")
+            lines.append("")
+            lines.append(f"![Optimization History]({optuna_plots['optimization_history']})")
+            lines.append("")
+            lines.append("Professional visualization of convergence over trials, showing both individual trial performance and best value progression.")
+            lines.append("")
+
+        # Parameter Importance
+        if 'param_importances' in optuna_plots:
+            lines.append("### Parameter Importance Analysis")
+            lines.append("")
+            lines.append(f"![Parameter Importance]({optuna_plots['param_importances']})")
+            lines.append("")
+            lines.append("This analysis quantifies which design variables have the most impact on the objective. Based on fANOVA (functional ANOVA) or other importance metrics.")
+            lines.append("")
+
+        # Slice Plot
+        if 'slice' in optuna_plots:
+            lines.append("### Parameter Slice Analysis")
+            lines.append("")
+            lines.append(f"![Slice Plot]({optuna_plots['slice']})")
+            lines.append("")
+            lines.append("Shows how changing each parameter individually affects the objective value, with other parameters held constant.")
+            lines.append("")
+
+        # Contour Plot
+        if 'contour' in optuna_plots:
+            lines.append("### Parameter Interaction Contour")
+            lines.append("")
+            lines.append(f"![Contour Plot]({optuna_plots['contour']})")
+            lines.append("")
+            lines.append("2D heatmap showing how combinations of two parameters affect the objective. Reveals interaction effects and optimal regions.")
+            lines.append("")
+
+    # Trial history table - show actual results
+    lines.append("## Complete Trial History")
+    lines.append("")
+    lines.append("<details>")
+    lines.append("<summary>Click to expand full trial history</summary>")
+    lines.append("")
+
+    # Build complete history table with results
+    history_header = ["Trial"] + result_col_names + [v.replace('_', ' ').title() for v in design_vars]
+    lines.append("| " + " | ".join(history_header) + " |")
+    lines.append("|" + "|".join(["-"*max(6, len(c)) for c in history_header]) + "|")
+
+    for trial in history:
+        result_vals = [f"{trial['results'][r]:.2f}" for r in result_cols]
+        var_vals = [f"{trial['design_variables'][v]:.2f}" for v in design_vars]
+        row_data = [f"#{trial['trial_number']}"] + result_vals + var_vals
+        lines.append("| " + " | ".join(row_data) + " |")
+
+    lines.append("")
+    lines.append("</details>")
+    lines.append("")
+
+    # Footer
+    lines.append("---")
+    lines.append("")
+    lines.append(f"*Report generated automatically by Atomizer optimization system*")
+
+    return '\n'.join(lines)
+
+
+def main():
+    """Command-line interface."""
+    if len(sys.argv) < 2:
+        print("Usage: python generate_report_markdown.py <history_file> [target_value] [tolerance]")
+        sys.exit(1)
+
+    history_file = Path(sys.argv[1])
+    if not history_file.exists():
+        print(f"Error: History file not found: {history_file}")
+        sys.exit(1)
+
+    target_value = float(sys.argv[2]) if len(sys.argv) > 2 else None
+    tolerance = float(sys.argv[3]) if len(sys.argv) > 3 else 0.1
+
+    # Generate report
+    report = generate_markdown_report(history_file, target_value, tolerance)
+
+    # Save report
+    report_file = history_file.parent / 'OPTIMIZATION_REPORT.md'
+    with open(report_file, 'w', encoding='utf-8') as f:
+        f.write(report)
+
+    print(f"Report saved to: {report_file}")
+
+
+if __name__ == '__main__':
+    main()

optimization_engine/reporting/report_generator.py

@@ -0,0 +1,152 @@
+"""
+Generate human-readable optimization reports from incremental history JSON.
+
+This script should be run automatically at the end of optimization, or manually
+to generate a report for any completed optimization study.
+"""
+
+import json
+import sys
+from pathlib import Path
+from typing import Dict, Any, List
+import numpy as np
+
+
+def generate_optimization_report(history_file: Path, target_value: float = None, tolerance: float = 0.1) -> str:
+    """
+    Generate a comprehensive human-readable optimization report.
+
+    Args:
+        history_file: Path to optimization_history_incremental.json
+        target_value: Target objective value (if applicable)
+        tolerance: Acceptable tolerance for success (default 0.1)
+
+    Returns:
+        Report text as a string
+    """
+    # Load history
+    with open(history_file) as f:
+        history = json.load(f)
+
+    if not history:
+        return "No optimization history found."
+
+    report = []
+    report.append('=' * 80)
+    report.append('OPTIMIZATION REPORT')
+    report.append('=' * 80)
+    report.append('')
+
+    # Study information
+    study_dir = history_file.parent.parent.parent
+    study_name = study_dir.name
+    report.append('STUDY INFORMATION')
+    report.append('-' * 80)
+    report.append(f'Study: {study_name}')
+    report.append(f'Total trials: {len(history)}')
+    report.append('')
+
+    # Design variables
+    first_trial = history[0]
+    design_vars = list(first_trial['design_variables'].keys())
+    report.append('DESIGN VARIABLES')
+    report.append('-' * 80)
+    for var in design_vars:
+        values = [t['design_variables'][var] for t in history]
+        report.append(f'  {var}:')
+        report.append(f'    Range: {min(values):.4f} - {max(values):.4f}')
+        report.append(f'    Mean: {np.mean(values):.4f}')
+    report.append('')
+
+    # Objective results
+    results = list(first_trial['results'].keys())
+    report.append('OBJECTIVE RESULTS')
+    report.append('-' * 80)
+    for result in results:
+        values = [t['results'][result] for t in history]
+        report.append(f'  {result}:')
+        report.append(f'    Range: {min(values):.4f} - {max(values):.4f}')
+        report.append(f'    Mean: {np.mean(values):.4f}')
+        report.append(f'    Std dev: {np.std(values):.4f}')
+    report.append('')
+
+    # Best trial
+    objectives = [t['objective'] for t in history]
+    best_trial = history[np.argmin(objectives)]
+
+    report.append('BEST TRIAL')
+    report.append('-' * 80)
+    report.append(f'Trial #{best_trial["trial_number"]}')
+    report.append(f'  Objective value: {best_trial["objective"]:.4f}')
+    report.append('  Design variables:')
+    for var, value in best_trial['design_variables'].items():
+        report.append(f'    {var}: {value:.4f}')
+    report.append('  Results:')
+    for result, value in best_trial['results'].items():
+        report.append(f'    {result}: {value:.4f}')
+    report.append('')
+
+    # Top 5 trials
+    report.append('TOP 5 TRIALS (by objective value)')
+    report.append('-' * 80)
+    sorted_history = sorted(history, key=lambda x: x['objective'])
+    for i, trial in enumerate(sorted_history[:5], 1):
+        report.append(f'{i}. Trial #{trial["trial_number"]}: Objective = {trial["objective"]:.4f}')
+        vars_str = ', '.join([f'{k}={v:.2f}' for k, v in trial['design_variables'].items()])
+        report.append(f'   {vars_str}')
+    report.append('')
+
+    # Success assessment (if target provided)
+    if target_value is not None:
+        report.append('SUCCESS ASSESSMENT')
+        report.append('-' * 80)
+        best_objective = min(objectives)
+        error = abs(best_objective - target_value)
+
+        if error <= tolerance:
+            report.append(f'[SUCCESS] Target {target_value} achieved within tolerance {tolerance}!')
+            report.append(f'  Best objective: {best_objective:.4f}')
+            report.append(f'  Error: {error:.4f}')
+        else:
+            report.append(f'[INCOMPLETE] Target {target_value} not achieved')
+            report.append(f'  Best objective: {best_objective:.4f}')
+            report.append(f'  Error: {error:.4f}')
+            report.append(f'  Need {error - tolerance:.4f} improvement')
+        report.append('')
+
+    report.append('=' * 80)
+
+    return '\n'.join(report)
+
+
+def main():
+    """Command-line interface for report generation."""
+    if len(sys.argv) < 2:
+        print("Usage: python generate_report.py <history_file> [target_value] [tolerance]")
+        print("Example: python generate_report.py studies/my_study/2_substudies/results/optimization_history_incremental.json 115.0 0.1")
+        sys.exit(1)
+
+    history_file = Path(sys.argv[1])
+    if not history_file.exists():
+        print(f"Error: History file not found: {history_file}")
+        sys.exit(1)
+
+    target_value = float(sys.argv[2]) if len(sys.argv) > 2 else None
+    tolerance = float(sys.argv[3]) if len(sys.argv) > 3 else 0.1
+
+    # Generate report
+    report = generate_optimization_report(history_file, target_value, tolerance)
+
+    # Save report
+    report_file = history_file.parent / 'OPTIMIZATION_REPORT.txt'
+    with open(report_file, 'w') as f:
+        f.write(report)
+
+    # Print to console
+    print(report)
+    print()
+    print(f"Report saved to: {report_file}")
+
+
+if __name__ == '__main__':
+    main()

optimization_engine/reporting/results_analyzer.py

@@ -0,0 +1,393 @@
+"""
+Comprehensive Results Analyzer
+
+Performs thorough introspection of OP2, F06, and other Nastran output files
+to discover ALL available results, not just what we expect.
+
+This helps ensure we don't miss important data that's actually in the output files.
+"""
+
+from pathlib import Path
+from typing import Dict, Any, List, Optional
+import json
+from dataclasses import dataclass, asdict
+from pyNastran.op2.op2 import OP2
+
+
+@dataclass
+class OP2Contents:
+    """Complete inventory of OP2 file contents."""
+    file_path: str
+    subcases: List[int]
+
+    # Displacement results
+    displacement_available: bool
+    displacement_subcases: List[int]
+
+    # Stress results (by element type)
+    stress_results: Dict[str, List[int]]  # element_type -> [subcases]
+
+    # Strain results (by element type)
+    strain_results: Dict[str, List[int]]
+
+    # Force results
+    force_results: Dict[str, List[int]]
+
+    # Other results
+    other_results: Dict[str, Any]
+
+    # Grid point forces/stresses
+    grid_point_forces: List[int]
+    spc_forces: List[int]
+    mpc_forces: List[int]
+
+    # Summary
+    total_result_types: int
+    element_types_with_results: List[str]
+
+
+@dataclass
+class F06Contents:
+    """Complete inventory of F06 file contents."""
+    file_path: str
+    has_displacement: bool
+    has_stress: bool
+    has_strain: bool
+    has_forces: bool
+    element_types_found: List[str]
+    error_messages: List[str]
+    warning_messages: List[str]
+
+
+class ComprehensiveResultsAnalyzer:
+    """
+    Analyzes ALL Nastran output files to discover available results.
+
+    This is much more thorough than just checking expected results.
+    """
+
+    def __init__(self, output_dir: Path):
+        """
+        Initialize analyzer.
+
+        Args:
+            output_dir: Directory containing Nastran output files (.op2, .f06, etc.)
+        """
+        self.output_dir = Path(output_dir)
+
+    def analyze_op2(self, op2_file: Path) -> OP2Contents:
+        """
+        Comprehensively analyze OP2 file contents.
+
+        Args:
+            op2_file: Path to OP2 file
+
+        Returns:
+            OP2Contents with complete inventory
+        """
+        print(f"\n[OP2 ANALYSIS] Reading: {op2_file.name}")
+
+        model = OP2()
+        model.read_op2(str(op2_file))
+
+        # Discover all subcases
+        all_subcases = set()
+
+        # Check displacement
+        displacement_available = hasattr(model, 'displacements') and len(model.displacements) > 0
+        displacement_subcases = list(model.displacements.keys()) if displacement_available else []
+        all_subcases.update(displacement_subcases)
+
+        print(f"  Displacement: {'YES' if displacement_available else 'NO'}")
+        if displacement_subcases:
+            print(f"    Subcases: {displacement_subcases}")
+
+        # Check ALL stress results by scanning attributes
+        stress_results = {}
+        element_types_with_stress = []
+
+        # List of known stress attribute names (safer than scanning all attributes)
+        stress_attrs = [
+            'cquad4_stress', 'ctria3_stress', 'ctetra_stress', 'chexa_stress', 'cpenta_stress',
+            'cbar_stress', 'cbeam_stress', 'crod_stress', 'conrod_stress', 'ctube_stress',
+            'cshear_stress', 'cbush_stress', 'cgap_stress', 'celas1_stress', 'celas2_stress',
+            'celas3_stress', 'celas4_stress'
+        ]
+
+        for attr_name in stress_attrs:
+            if hasattr(model, attr_name):
+                try:
+                    stress_obj = getattr(model, attr_name)
+                    if isinstance(stress_obj, dict) and len(stress_obj) > 0:
+                        element_type = attr_name.replace('_stress', '')
+                        subcases = list(stress_obj.keys())
+                        stress_results[element_type] = subcases
+                        element_types_with_stress.append(element_type)
+                        all_subcases.update(subcases)
+                        print(f"  Stress [{element_type}]: YES")
+                        print(f"    Subcases: {subcases}")
+                except Exception as e:
+                    # Skip attributes that cause errors
+                    pass
+
+        if not stress_results:
+            print(f"  Stress: NO stress results found")
+
+        # Check ALL strain results
+        strain_results = {}
+        strain_attrs = [attr.replace('_stress', '_strain') for attr in stress_attrs]
+
+        for attr_name in strain_attrs:
+            if hasattr(model, attr_name):
+                try:
+                    strain_obj = getattr(model, attr_name)
+                    if isinstance(strain_obj, dict) and len(strain_obj) > 0:
+                        element_type = attr_name.replace('_strain', '')
+                        subcases = list(strain_obj.keys())
+                        strain_results[element_type] = subcases
+                        all_subcases.update(subcases)
+                        print(f"  Strain [{element_type}]: YES")
+                        print(f"    Subcases: {subcases}")
+                except Exception as e:
+                    pass
+
+        if not strain_results:
+            print(f"  Strain: NO strain results found")
+
+        # Check ALL force results
+        force_results = {}
+        force_attrs = [attr.replace('_stress', '_force') for attr in stress_attrs]
+
+        for attr_name in force_attrs:
+            if hasattr(model, attr_name):
+                try:
+                    force_obj = getattr(model, attr_name)
+                    if isinstance(force_obj, dict) and len(force_obj) > 0:
+                        element_type = attr_name.replace('_force', '')
+                        subcases = list(force_obj.keys())
+                        force_results[element_type] = subcases
+                        all_subcases.update(subcases)
+                        print(f"  Force [{element_type}]: YES")
+                        print(f"    Subcases: {subcases}")
+                except Exception as e:
+                    pass
+
+        if not force_results:
+            print(f"  Force: NO force results found")
+
+        # Check grid point forces
+        grid_point_forces = list(model.grid_point_forces.keys()) if hasattr(model, 'grid_point_forces') else []
+        if grid_point_forces:
+            print(f"  Grid Point Forces: YES")
+            print(f"    Subcases: {grid_point_forces}")
+            all_subcases.update(grid_point_forces)
+
+        # Check SPC/MPC forces
+        spc_forces = list(model.spc_forces.keys()) if hasattr(model, 'spc_forces') else []
+        mpc_forces = list(model.mpc_forces.keys()) if hasattr(model, 'mpc_forces') else []
+
+        if spc_forces:
+            print(f"  SPC Forces: YES")
+            print(f"    Subcases: {spc_forces}")
+            all_subcases.update(spc_forces)
+
+        if mpc_forces:
+            print(f"  MPC Forces: YES")
+            print(f"    Subcases: {mpc_forces}")
+            all_subcases.update(mpc_forces)
+
+        # Check for other interesting results
+        other_results = {}
+        interesting_attrs = ['eigenvalues', 'eigenvectors', 'thermal_load_vectors',
+                           'load_vectors', 'contact', 'glue', 'slide_lines']
+
+        for attr_name in interesting_attrs:
+            if hasattr(model, attr_name):
+                obj = getattr(model, attr_name)
+                if obj and (isinstance(obj, dict) and len(obj) > 0) or (not isinstance(obj, dict)):
+                    other_results[attr_name] = str(type(obj))
+                    print(f"  {attr_name}: YES")
+
+        # Collect all element types that have any results
+        all_element_types = set()
+        all_element_types.update(stress_results.keys())
+        all_element_types.update(strain_results.keys())
+        all_element_types.update(force_results.keys())
+
+        total_result_types = (
+            len(stress_results) +
+            len(strain_results) +
+            len(force_results) +
+            (1 if displacement_available else 0) +
+            (1 if grid_point_forces else 0) +
+            (1 if spc_forces else 0) +
+            (1 if mpc_forces else 0) +
+            len(other_results)
+        )
+
+        print(f"\n  SUMMARY:")
+        print(f"    Total subcases: {len(all_subcases)}")
+        print(f"    Total result types: {total_result_types}")
+        print(f"    Element types with results: {sorted(all_element_types)}")
+
+        return OP2Contents(
+            file_path=str(op2_file),
+            subcases=sorted(all_subcases),
+            displacement_available=displacement_available,
+            displacement_subcases=displacement_subcases,
+            stress_results=stress_results,
+            strain_results=strain_results,
+            force_results=force_results,
+            other_results=other_results,
+            grid_point_forces=grid_point_forces,
+            spc_forces=spc_forces,
+            mpc_forces=mpc_forces,
+            total_result_types=total_result_types,
+            element_types_with_results=sorted(all_element_types)
+        )
+
+    def analyze_f06(self, f06_file: Path) -> F06Contents:
+        """
+        Analyze F06 file for available results.
+
+        Args:
+            f06_file: Path to F06 file
+
+        Returns:
+            F06Contents with inventory
+        """
+        print(f"\n[F06 ANALYSIS] Reading: {f06_file.name}")
+
+        if not f06_file.exists():
+            print(f"  F06 file not found")
+            return F06Contents(
+                file_path=str(f06_file),
+                has_displacement=False,
+                has_stress=False,
+                has_strain=False,
+                has_forces=False,
+                element_types_found=[],
+                error_messages=[],
+                warning_messages=[]
+            )
+
+        # Read F06 file
+        with open(f06_file, 'r', encoding='latin-1', errors='ignore') as f:
+            content = f.read()
+
+        # Search for key sections
+        has_displacement = 'D I S P L A C E M E N T' in content
+        has_stress = 'S T R E S S E S' in content
+        has_strain = 'S T R A I N S' in content
+        has_forces = 'F O R C E S' in content
+
+        print(f"  Displacement: {'YES' if has_displacement else 'NO'}")
+        print(f"  Stress: {'YES' if has_stress else 'NO'}")
+        print(f"  Strain: {'YES' if has_strain else 'NO'}")
+        print(f"  Forces: {'YES' if has_forces else 'NO'}")
+
+        # Find element types mentioned
+        element_keywords = ['CQUAD4', 'CTRIA3', 'CTETRA', 'CHEXA', 'CPENTA', 'CBAR', 'CBEAM', 'CROD']
+        element_types_found = []
+
+        for elem_type in element_keywords:
+            if elem_type in content:
+                element_types_found.append(elem_type)
+
+        if element_types_found:
+            print(f"  Element types: {element_types_found}")
+
+        # Extract errors and warnings
+        error_messages = []
+        warning_messages = []
+
+        for line in content.split('\n'):
+            line_upper = line.upper()
+            if 'ERROR' in line_upper or 'FATAL' in line_upper:
+                error_messages.append(line.strip())
+            elif 'WARNING' in line_upper or 'WARN' in line_upper:
+                warning_messages.append(line.strip())
+
+        if error_messages:
+            print(f"  Errors found: {len(error_messages)}")
+            for err in error_messages[:5]:  # Show first 5
+                print(f"    {err}")
+
+        if warning_messages:
+            print(f"  Warnings found: {len(warning_messages)}")
+            for warn in warning_messages[:5]:  # Show first 5
+                print(f"    {warn}")
+
+        return F06Contents(
+            file_path=str(f06_file),
+            has_displacement=has_displacement,
+            has_stress=has_stress,
+            has_strain=has_strain,
+            has_forces=has_forces,
+            element_types_found=element_types_found,
+            error_messages=error_messages[:20],  # Keep first 20
+            warning_messages=warning_messages[:20]
+        )
+
+    def analyze_all(self, op2_pattern: str = "*.op2", f06_pattern: str = "*.f06") -> Dict[str, Any]:
+        """
+        Analyze all OP2 and F06 files in directory.
+
+        Args:
+            op2_pattern: Glob pattern for OP2 files
+            f06_pattern: Glob pattern for F06 files
+
+        Returns:
+            Dict with complete analysis results
+        """
+        print("="*80)
+        print("COMPREHENSIVE NASTRAN RESULTS ANALYSIS")
+        print("="*80)
+        print(f"\nDirectory: {self.output_dir}")
+
+        results = {
+            'directory': str(self.output_dir),
+            'op2_files': [],
+            'f06_files': []
+        }
+
+        # Find and analyze all OP2 files
+        op2_files = list(self.output_dir.glob(op2_pattern))
+        print(f"\nFound {len(op2_files)} OP2 file(s)")
+
+        for op2_file in op2_files:
+            op2_contents = self.analyze_op2(op2_file)
+            results['op2_files'].append(asdict(op2_contents))
+
+        # Find and analyze all F06 files
+        f06_files = list(self.output_dir.glob(f06_pattern))
+        print(f"\nFound {len(f06_files)} F06 file(s)")
+
+        for f06_file in f06_files:
+            f06_contents = self.analyze_f06(f06_file)
+            results['f06_files'].append(asdict(f06_contents))
+
+        print("\n" + "="*80)
+        print("ANALYSIS COMPLETE")
+        print("="*80)
+
+        return results
+
+
+if __name__ == '__main__':
+    import sys
+
+    if len(sys.argv) > 1:
+        output_dir = Path(sys.argv[1])
+    else:
+        output_dir = Path.cwd()
+
+    analyzer = ComprehensiveResultsAnalyzer(output_dir)
+    results = analyzer.analyze_all()
+
+    # Save results to JSON
+    output_file = output_dir / "comprehensive_results_analysis.json"
+    with open(output_file, 'w', encoding='utf-8') as f:
+        json.dump(results, f, indent=2)
+
+    print(f"\nResults saved to: {output_file}")

optimization_engine/reporting/visualizer.py

@@ -0,0 +1,555 @@
+"""
+Optimization Visualization System
+
+Generates publication-quality plots for optimization results:
+- Convergence plots
+- Design space exploration
+- Parallel coordinate plots
+- Parameter sensitivity heatmaps
+- Constraint violation tracking
+"""
+
+from pathlib import Path
+from typing import Dict, List, Any, Optional
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+from matplotlib.figure import Figure
+import pandas as pd
+from datetime import datetime
+
+# Configure matplotlib for publication quality
+mpl.rcParams['figure.dpi'] = 150
+mpl.rcParams['savefig.dpi'] = 300
+mpl.rcParams['font.size'] = 10
+mpl.rcParams['font.family'] = 'sans-serif'
+mpl.rcParams['axes.labelsize'] = 10
+mpl.rcParams['axes.titlesize'] = 11
+mpl.rcParams['xtick.labelsize'] = 9
+mpl.rcParams['ytick.labelsize'] = 9
+mpl.rcParams['legend.fontsize'] = 9
+
+
+class OptimizationVisualizer:
+    """
+    Generate comprehensive visualizations for optimization studies.
+
+    Automatically creates:
+    - Convergence plot (objective vs trials)
+    - Design space exploration (parameter evolution)
+    - Parallel coordinate plot (high-dimensional view)
+    - Sensitivity heatmap (correlations)
+    - Constraint violation tracking
+    """
+
+    def __init__(self, substudy_dir: Path):
+        """
+        Initialize visualizer for a substudy.
+
+        Args:
+            substudy_dir: Path to substudy directory containing history.json
+        """
+        self.substudy_dir = Path(substudy_dir)
+        self.plots_dir = self.substudy_dir / 'plots'
+        self.plots_dir.mkdir(exist_ok=True)
+
+        # Load data
+        self.history = self._load_history()
+        self.config = self._load_config()
+        self.df = self._history_to_dataframe()
+
+    def _load_history(self) -> List[Dict]:
+        """Load optimization history from JSON."""
+        history_file = self.substudy_dir / 'history.json'
+        if not history_file.exists():
+            raise FileNotFoundError(f"History file not found: {history_file}")
+
+        with open(history_file, 'r') as f:
+            return json.load(f)
+
+    def _load_config(self) -> Dict:
+        """Load optimization configuration."""
+        # Try to find config in parent directories
+        for parent in [self.substudy_dir, self.substudy_dir.parent, self.substudy_dir.parent.parent]:
+            config_files = list(parent.glob('*config.json'))
+            if config_files:
+                with open(config_files[0], 'r') as f:
+                    return json.load(f)
+
+        # Return minimal config if not found
+        return {'design_variables': {}, 'objectives': [], 'constraints': []}
+
+    def _history_to_dataframe(self) -> pd.DataFrame:
+        """Convert history to flat DataFrame for analysis."""
+        rows = []
+        for entry in self.history:
+            row = {
+                'trial': entry.get('trial_number'),
+                'timestamp': entry.get('timestamp'),
+                'total_objective': entry.get('total_objective')
+            }
+
+            # Add design variables
+            for var, val in entry.get('design_variables', {}).items():
+                row[f'dv_{var}'] = val
+
+            # Add objectives
+            for obj, val in entry.get('objectives', {}).items():
+                row[f'obj_{obj}'] = val
+
+            # Add constraints
+            for const, val in entry.get('constraints', {}).items():
+                row[f'const_{const}'] = val
+
+            rows.append(row)
+
+        return pd.DataFrame(rows)
+
+    def generate_all_plots(self, save_formats: List[str] = ['png', 'pdf']) -> Dict[str, List[Path]]:
+        """
+        Generate all visualization plots.
+
+        Args:
+            save_formats: List of formats to save plots in (png, pdf, svg)
+
+        Returns:
+            Dictionary mapping plot type to list of saved file paths
+        """
+        saved_files = {}
+
+        print(f"Generating plots in: {self.plots_dir}")
+
+        # 1. Convergence plot
+        print("  - Generating convergence plot...")
+        saved_files['convergence'] = self.plot_convergence(save_formats)
+
+        # 2. Design space exploration
+        print("  - Generating design space exploration...")
+        saved_files['design_space'] = self.plot_design_space(save_formats)
+
+        # 3. Parallel coordinate plot
+        print("  - Generating parallel coordinate plot...")
+        saved_files['parallel_coords'] = self.plot_parallel_coordinates(save_formats)
+
+        # 4. Sensitivity heatmap
+        print("  - Generating sensitivity heatmap...")
+        saved_files['sensitivity'] = self.plot_sensitivity_heatmap(save_formats)
+
+        # 5. Constraint violations (if constraints exist)
+        if any('const_' in col for col in self.df.columns):
+            print("  - Generating constraint violation plot...")
+            saved_files['constraints'] = self.plot_constraint_violations(save_formats)
+
+        # 6. Objective breakdown (if multi-objective)
+        obj_cols = [col for col in self.df.columns if col.startswith('obj_')]
+        if len(obj_cols) > 1:
+            print("  - Generating objective breakdown...")
+            saved_files['objectives'] = self.plot_objective_breakdown(save_formats)
+
+        print(f"SUCCESS: All plots saved to: {self.plots_dir}")
+        return saved_files
+
+    def plot_convergence(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Plot optimization convergence: objective value vs trial number.
+        Shows both individual trials and running best.
+        """
+        fig, ax = plt.subplots(figsize=(10, 6))
+
+        trials = self.df['trial'].values
+        objectives = self.df['total_objective'].values
+
+        # Calculate running best
+        running_best = np.minimum.accumulate(objectives)
+
+        # Plot individual trials
+        ax.scatter(trials, objectives, alpha=0.6, s=30, color='steelblue',
+                   label='Trial objective', zorder=2)
+
+        # Plot running best
+        ax.plot(trials, running_best, color='darkred', linewidth=2,
+                label='Running best', zorder=3)
+
+        # Highlight best trial
+        best_idx = np.argmin(objectives)
+        ax.scatter(trials[best_idx], objectives[best_idx],
+                   color='gold', s=200, marker='*', edgecolors='black',
+                   linewidths=1.5, label='Best trial', zorder=4)
+
+        ax.set_xlabel('Trial Number')
+        ax.set_ylabel('Total Objective Value')
+        ax.set_title('Optimization Convergence')
+        ax.legend(loc='best')
+        ax.grid(True, alpha=0.3)
+
+        # Add improvement annotation
+        improvement = (objectives[0] - objectives[best_idx]) / objectives[0] * 100
+        ax.text(0.02, 0.98, f'Improvement: {improvement:.1f}%\nBest trial: {trials[best_idx]}',
+                transform=ax.transAxes, verticalalignment='top',
+                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
+
+        plt.tight_layout()
+        return self._save_figure(fig, 'convergence', save_formats)
+
+    def plot_design_space(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Plot design variable evolution over trials.
+        Shows how parameters change during optimization.
+        """
+        dv_cols = [col for col in self.df.columns if col.startswith('dv_')]
+        n_vars = len(dv_cols)
+
+        if n_vars == 0:
+            print("  Warning: No design variables found, skipping design space plot")
+            return []
+
+        # Create subplots
+        fig, axes = plt.subplots(n_vars, 1, figsize=(10, 3*n_vars), sharex=True)
+        if n_vars == 1:
+            axes = [axes]
+
+        trials = self.df['trial'].values
+        objectives = self.df['total_objective'].values
+        best_idx = np.argmin(objectives)
+
+        for idx, col in enumerate(dv_cols):
+            ax = axes[idx]
+            var_name = col.replace('dv_', '')
+            values = self.df[col].values
+
+            # Color points by objective value (normalized)
+            norm = mpl.colors.Normalize(vmin=objectives.min(), vmax=objectives.max())
+            colors = plt.cm.viridis_r(norm(objectives))  # reversed so better = darker
+
+            # Plot evolution
+            scatter = ax.scatter(trials, values, c=colors, s=40, alpha=0.7,
+                                edgecolors='black', linewidths=0.5)
+
+            # Highlight best trial
+            ax.scatter(trials[best_idx], values[best_idx],
+                       color='gold', s=200, marker='*', edgecolors='black',
+                       linewidths=1.5, zorder=10)
+
+            # Get units from config
+            units = self.config.get('design_variables', {}).get(var_name, {}).get('units', '')
+            ylabel = f'{var_name}'
+            if units:
+                ylabel += f' [{units}]'
+
+            ax.set_ylabel(ylabel)
+            ax.grid(True, alpha=0.3)
+
+            # Add colorbar for first subplot
+            if idx == 0:
+                cbar = plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap='viridis_r'),
+                                    ax=ax, orientation='horizontal', pad=0.1)
+                cbar.set_label('Objective Value (darker = better)')
+
+        axes[-1].set_xlabel('Trial Number')
+        fig.suptitle('Design Space Exploration', fontsize=12, y=1.0)
+        plt.tight_layout()
+
+        return self._save_figure(fig, 'design_space_evolution', save_formats)
+
+    def plot_parallel_coordinates(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Parallel coordinate plot showing high-dimensional design space.
+        Each line represents one trial, colored by objective value.
+        """
+        # Get design variables and objective
+        dv_cols = [col for col in self.df.columns if col.startswith('dv_')]
+
+        if len(dv_cols) == 0:
+            print("  Warning: No design variables found, skipping parallel coordinates plot")
+            return []
+
+        # Prepare data: normalize all columns to [0, 1]
+        plot_data = self.df[dv_cols + ['total_objective']].copy()
+
+        # Normalize each column
+        normalized = pd.DataFrame()
+        for col in plot_data.columns:
+            col_min = plot_data[col].min()
+            col_max = plot_data[col].max()
+            if col_max > col_min:
+                normalized[col] = (plot_data[col] - col_min) / (col_max - col_min)
+            else:
+                normalized[col] = 0.5  # If constant, put in middle
+
+        # Create figure
+        fig, ax = plt.subplots(figsize=(12, 6))
+
+        # Setup x-axis
+        n_vars = len(normalized.columns)
+        x_positions = np.arange(n_vars)
+
+        # Color by objective value
+        objectives = self.df['total_objective'].values
+        norm = mpl.colors.Normalize(vmin=objectives.min(), vmax=objectives.max())
+        colormap = plt.cm.viridis_r
+
+        # Plot each trial as a line
+        for idx in range(len(normalized)):
+            values = normalized.iloc[idx].values
+            color = colormap(norm(objectives[idx]))
+            ax.plot(x_positions, values, color=color, alpha=0.3, linewidth=1)
+
+        # Highlight best trial
+        best_idx = np.argmin(objectives)
+        best_values = normalized.iloc[best_idx].values
+        ax.plot(x_positions, best_values, color='gold', linewidth=3,
+                label='Best trial', zorder=10, marker='o', markersize=8,
+                markeredgecolor='black', markeredgewidth=1.5)
+
+        # Setup axes
+        ax.set_xticks(x_positions)
+        labels = [col.replace('dv_', '').replace('_', '\n') for col in dv_cols] + ['Objective']
+        ax.set_xticklabels(labels, rotation=0, ha='center')
+        ax.set_ylabel('Normalized Value [0-1]')
+        ax.set_title('Parallel Coordinate Plot - Design Space Overview')
+        ax.set_ylim(-0.05, 1.05)
+        ax.grid(True, alpha=0.3, axis='y')
+        ax.legend(loc='best')
+
+        # Add colorbar
+        sm = mpl.cm.ScalarMappable(cmap=colormap, norm=norm)
+        sm.set_array([])
+        cbar = plt.colorbar(sm, ax=ax, orientation='vertical', pad=0.02)
+        cbar.set_label('Objective Value (darker = better)')
+
+        plt.tight_layout()
+        return self._save_figure(fig, 'parallel_coordinates', save_formats)
+
+    def plot_sensitivity_heatmap(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Correlation heatmap showing sensitivity between design variables and objectives.
+        """
+        # Get numeric columns
+        dv_cols = [col for col in self.df.columns if col.startswith('dv_')]
+        obj_cols = [col for col in self.df.columns if col.startswith('obj_')]
+
+        if not dv_cols or not obj_cols:
+            print("  Warning: Insufficient data for sensitivity heatmap, skipping")
+            return []
+
+        # Calculate correlation matrix
+        analysis_cols = dv_cols + obj_cols + ['total_objective']
+        corr_matrix = self.df[analysis_cols].corr()
+
+        # Extract DV vs Objective correlations
+        sensitivity = corr_matrix.loc[dv_cols, obj_cols + ['total_objective']]
+
+        # Create heatmap
+        fig, ax = plt.subplots(figsize=(10, max(6, len(dv_cols) * 0.6)))
+
+        im = ax.imshow(sensitivity.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
+
+        # Set ticks
+        ax.set_xticks(np.arange(len(sensitivity.columns)))
+        ax.set_yticks(np.arange(len(sensitivity.index)))
+
+        # Labels
+        x_labels = [col.replace('obj_', '').replace('_', ' ') for col in sensitivity.columns]
+        y_labels = [col.replace('dv_', '').replace('_', ' ') for col in sensitivity.index]
+        ax.set_xticklabels(x_labels, rotation=45, ha='right')
+        ax.set_yticklabels(y_labels)
+
+        # Add correlation values as text
+        for i in range(len(sensitivity.index)):
+            for j in range(len(sensitivity.columns)):
+                value = sensitivity.values[i, j]
+                color = 'white' if abs(value) > 0.5 else 'black'
+                ax.text(j, i, f'{value:.2f}', ha='center', va='center',
+                       color=color, fontsize=9)
+
+        ax.set_title('Parameter Sensitivity Analysis\n(Correlation: Design Variables vs Objectives)')
+
+        # Colorbar
+        cbar = plt.colorbar(im, ax=ax)
+        cbar.set_label('Correlation Coefficient', rotation=270, labelpad=20)
+
+        plt.tight_layout()
+        return self._save_figure(fig, 'sensitivity_heatmap', save_formats)
+
+    def plot_constraint_violations(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Plot constraint violations over trials.
+        """
+        const_cols = [col for col in self.df.columns if col.startswith('const_')]
+
+        if not const_cols:
+            return []
+
+        fig, ax = plt.subplots(figsize=(10, 6))
+
+        trials = self.df['trial'].values
+
+        for col in const_cols:
+            const_name = col.replace('const_', '').replace('_', ' ')
+            values = self.df[col].values
+
+            # Plot constraint value
+            ax.plot(trials, values, marker='o', markersize=4,
+                   label=const_name, alpha=0.7, linewidth=1.5)
+
+        ax.axhline(y=0, color='red', linestyle='--', linewidth=2,
+                   label='Feasible threshold', zorder=1)
+
+        ax.set_xlabel('Trial Number')
+        ax.set_ylabel('Constraint Value (< 0 = satisfied)')
+        ax.set_title('Constraint Violations Over Trials')
+        ax.legend(loc='best')
+        ax.grid(True, alpha=0.3)
+
+        plt.tight_layout()
+        return self._save_figure(fig, 'constraint_violations', save_formats)
+
+    def plot_objective_breakdown(self, save_formats: List[str] = ['png']) -> List[Path]:
+        """
+        Stacked area plot showing individual objective contributions.
+        """
+        obj_cols = [col for col in self.df.columns if col.startswith('obj_')]
+
+        if len(obj_cols) < 2:
+            return []
+
+        fig, ax = plt.subplots(figsize=(10, 6))
+
+        trials = self.df['trial'].values
+
+        # Normalize objectives for stacking
+        obj_data = self.df[obj_cols].values.T
+
+        ax.stackplot(trials, *obj_data,
+                     labels=[col.replace('obj_', '').replace('_', ' ') for col in obj_cols],
+                     alpha=0.7)
+
+        # Also plot total
+        ax.plot(trials, self.df['total_objective'].values,
+               color='black', linewidth=2, linestyle='--',
+               label='Total objective', zorder=10)
+
+        ax.set_xlabel('Trial Number')
+        ax.set_ylabel('Objective Value')
+        ax.set_title('Multi-Objective Breakdown')
+        ax.legend(loc='best')
+        ax.grid(True, alpha=0.3)
+
+        plt.tight_layout()
+        return self._save_figure(fig, 'objective_breakdown', save_formats)
+
+    def _save_figure(self, fig: Figure, name: str, formats: List[str]) -> List[Path]:
+        """
+        Save figure in multiple formats.
+
+        Args:
+            fig: Matplotlib figure
+            name: Base filename (without extension)
+            formats: List of file formats (png, pdf, svg)
+
+        Returns:
+            List of saved file paths
+        """
+        saved_paths = []
+        for fmt in formats:
+            filepath = self.plots_dir / f'{name}.{fmt}'
+            fig.savefig(filepath, bbox_inches='tight')
+            saved_paths.append(filepath)
+
+        plt.close(fig)
+        return saved_paths
+
+    def generate_plot_summary(self) -> Dict[str, Any]:
+        """
+        Generate summary statistics for inclusion in reports.
+
+        Returns:
+            Dictionary with key statistics and insights
+        """
+        objectives = self.df['total_objective'].values
+        trials = self.df['trial'].values
+
+        best_idx = np.argmin(objectives)
+        best_trial = int(trials[best_idx])
+        best_value = float(objectives[best_idx])
+        initial_value = float(objectives[0])
+        improvement_pct = (initial_value - best_value) / initial_value * 100
+
+        # Convergence metrics
+        running_best = np.minimum.accumulate(objectives)
+        improvements = np.diff(running_best)
+        significant_improvements = np.sum(improvements < -0.01 * initial_value)  # >1% improvement
+
+        # Design variable ranges
+        dv_cols = [col for col in self.df.columns if col.startswith('dv_')]
+        dv_exploration = {}
+        for col in dv_cols:
+            var_name = col.replace('dv_', '')
+            values = self.df[col].values
+            dv_exploration[var_name] = {
+                'min_explored': float(values.min()),
+                'max_explored': float(values.max()),
+                'best_value': float(values[best_idx]),
+                'range_coverage': float((values.max() - values.min()))
+            }
+
+        summary = {
+            'total_trials': int(len(trials)),
+            'best_trial': best_trial,
+            'best_objective': best_value,
+            'initial_objective': initial_value,
+            'improvement_percent': improvement_pct,
+            'significant_improvements': int(significant_improvements),
+            'design_variable_exploration': dv_exploration,
+            'convergence_rate': float(np.mean(np.abs(improvements[:10]))) if len(improvements) > 10 else 0.0,
+            'timestamp': datetime.now().isoformat()
+        }
+
+        # Save summary
+        summary_file = self.plots_dir / 'plot_summary.json'
+        with open(summary_file, 'w') as f:
+            json.dump(summary, f, indent=2)
+
+        return summary
+
+
+def generate_plots_for_substudy(substudy_dir: Path, formats: List[str] = ['png', 'pdf']):
+    """
+    Convenience function to generate all plots for a substudy.
+
+    Args:
+        substudy_dir: Path to substudy directory
+        formats: List of save formats
+
+    Returns:
+        OptimizationVisualizer instance
+    """
+    visualizer = OptimizationVisualizer(substudy_dir)
+    visualizer.generate_all_plots(save_formats=formats)
+    summary = visualizer.generate_plot_summary()
+
+    print(f"\n{'='*60}")
+    print(f"VISUALIZATION SUMMARY")
+    print(f"{'='*60}")
+    print(f"Total trials: {summary['total_trials']}")
+    print(f"Best trial: {summary['best_trial']}")
+    print(f"Improvement: {summary['improvement_percent']:.2f}%")
+    print(f"Plots saved to: {visualizer.plots_dir}")
+    print(f"{'='*60}\n")
+
+    return visualizer
+
+
+if __name__ == '__main__':
+    import sys
+
+    if len(sys.argv) < 2:
+        print("Usage: python visualizer.py <substudy_directory> [formats...]")
+        print("Example: python visualizer.py studies/beam/substudies/opt1 png pdf")
+        sys.exit(1)
+
+    substudy_path = Path(sys.argv[1])
+    formats = sys.argv[2:] if len(sys.argv) > 2 else ['png', 'pdf']
+
+    generate_plots_for_substudy(substudy_path, formats)