Logging and Analysis

# Main logging and plotting functions available in cma module
import cma

def disp(iteration=None):
    """
    Display current optimization state and progress.
    
    Parameters:
    -----------
    iteration : int, optional
        Iteration number to display (default: current).
    """
    pass

def plot(name_prefix='outcmaes', **kwargs):
    """
    Plot CMA-ES data from logged files.
    
    Parameters:
    -----------
    name_prefix : str, optional
        File prefix for data files to plot (default 'outcmaes').
    **kwargs : dict
        Additional plotting options.
    """
    pass

def plot_zip(name_prefix='outcmaes', **kwargs):
    """
    Plot multiple CMA-ES runs from zipped data files.
    
    Parameters:
    -----------
    name_prefix : str, optional
        File prefix pattern for multiple runs.
    """
    pass

class CMADataLogger:
    """
    Data logger for CMAEvolutionStrategy with automatic file management.
    
    Logs optimization data to files and provides loading, visualization,
    and analysis capabilities. Data is written to disk in real-time and
    can be reloaded for post-processing and analysis.
    """
    
    default_prefix = 'outcmaes/'  # Default file prefix
    
    def __init__(self, name_prefix=None, modulo=1, append=False, expensive_modulo=1):
        """
        Initialize CMA-ES data logger.
        
        Parameters:
        -----------
        name_prefix : str, optional
            File name prefix for all data files. If None, uses default 'outcmaes/'.
            Can be a directory path ending with '/' for organized storage.
            
        modulo : int, optional
            Log data every modulo iterations (default 1 = every iteration).
            
        append : bool, optional
            Whether to append to existing files (default False = overwrite).
            
        expensive_modulo : int, optional
            Frequency for expensive operations like eigendecomposition logging
            (default 1). Set higher for large dimensions to reduce overhead.
            
        Examples:
        ---------
        >>> import cma
        >>> 
        >>> # Basic logging
        >>> logger = cma.CMADataLogger()
        >>> es = cma.CMAEvolutionStrategy([0, 0, 0], 0.5)
        >>> logger.register(es)  # Connect logger to optimizer
        >>> 
        >>> while not es.stop():
        ...     solutions = es.ask()
        ...     fitness_values = [sum(x**2) for x in solutions]
        ...     es.tell(solutions, fitness_values)
        ...     logger.add()  # Log current state
        >>> 
        >>> # Custom file prefix
        >>> logger = cma.CMADataLogger('my_experiment/run_001')
        >>> 
        >>> # Reduced logging frequency for large problems
        >>> logger = cma.CMADataLogger(modulo=10, expensive_modulo=50)
        """
        pass
    
    def register(self, es, append=None, modulo=None):
        """
        Register evolution strategy for automatic logging.
        
        Parameters:
        -----------
        es : CMAEvolutionStrategy
            Evolution strategy instance to log.
            
        append : bool, optional
            Override constructor append setting.
            
        modulo : int, optional
            Override constructor modulo setting.
            
        Returns:
        --------
        CMADataLogger
            Self for method chaining.
            
        Examples:
        ---------
        >>> import cma
        >>> 
        >>> es = cma.CMAEvolutionStrategy([0, 0, 0], 0.5)
        >>> logger = cma.CMADataLogger().register(es)
        >>> 
        >>> # Method chaining
        >>> logger = cma.CMADataLogger('experiment1/').register(es, modulo=5)
        """
        pass
    
    def add(self, es=None, more_data=[], modulo=None):
        """
        Add current evolution strategy state to log.
        
        Parameters:
        -----------
        es : CMAEvolutionStrategy, optional
            Evolution strategy to log. If None, uses registered ES.
            
        more_data : list, optional
            Additional data columns to append to fit file.
            
        modulo : int, optional
            Override default logging frequency.
            
        Examples:
        ---------
        >>> import cma
        >>> 
        >>> es = cma.CMAEvolutionStrategy([0, 0, 0], 0.5)
        >>> logger = cma.CMADataLogger().register(es)
        >>> 
        >>> while not es.stop():
        ...     solutions = es.ask()
        ...     fitness_values = [sum(x**2) for x in solutions]
        ...     es.tell(solutions, fitness_values)
        ...     
        ...     # Log with additional data
        ...     custom_metric = es.sigma * es.result.evaluations
        ...     logger.add(more_data=[custom_metric])
        """
        pass
    
    def load(self, filenameprefix=None):
        """
        Load logged data from files.
        
        Parameters:
        -----------
        filenameprefix : str, optional
            File prefix to load. If None, uses instance prefix.
            
        Returns:
        --------
        CMADataLogger
            Self with loaded data in attributes.
            
        Examples:
        ---------
        >>> import cma
        >>> 
        >>> # Load previously logged data
        >>> logger = cma.CMADataLogger().load()
        >>> 
        >>> # Load specific experiment
        >>> logger = cma.CMADataLogger().load('experiment1/run_005')
        >>> 
        >>> # Access loaded data
        >>> print(f"Iterations: {len(logger.f)}")
        >>> print(f"Final fitness: {logger.f[-1, 5]}")
        >>> print(f"Final mean: {logger.xmean[-1]}")
        """
        pass
    
    def save_to(self, filenameprefix):
        """
        Save current data to different file location.
        
        Parameters:
        -----------
        filenameprefix : str
            New file prefix for saving data.
        """
        pass
    
    @property
    def data(self):
        """
        Dictionary of loaded data arrays.
        
        Returns:
        --------
        dict
            Data dictionary with keys: 'xmean', 'xrecent', 'std', 'f', 'D', 'corrspec'
            corresponding to mean solutions, recent best, standard deviations,
            fitness values, axis lengths, and correlation spectrum.
        """
        pass

# Data files and their contents
data_files = {
    'fit': """
    Fitness and evaluation data (outcmaesfit.dat):
    Columns: [iteration, evaluations, sigma, axis_ratio, best_fitness, 
              median_fitness, worst_fitness, further_data...]
    
    Example access:
    >>> logger = cma.CMADataLogger().load()
    >>> iterations = logger.f[:, 0]
    >>> evaluations = logger.f[:, 1] 
    >>> best_fitness = logger.f[:, 4]
    >>> median_fitness = logger.f[:, 5]
    """,
    
    'xmean': """
    Mean solution evolution (outcmaesxmean.dat):
    Columns: [iteration, evaluations, sigma, axis_ratio, mean_x1, mean_x2, ...]
    
    Example access:
    >>> logger = cma.CMADataLogger().load()
    >>> mean_solutions = logger.xmean[:, 4:]  # Skip iteration info
    >>> final_mean = logger.xmean[-1, 4:]     # Final mean solution
    """,
    
    'xrecentbest': """
    Recent best solutions (outcmaesxrecentbest.dat):
    Columns: [iteration, evaluations, sigma, axis_ratio, x1, x2, ...]
    
    Example access:
    >>> logger = cma.CMADataLogger().load() 
    >>> recent_best = logger.xrecent[:, 4:]   # Recent best solutions
    >>> best_solution = logger.xrecent[-1, 4:]  # Final best solution
    """,
    
    'stddev': """
    Standard deviations (outcmaesstddev.dat):
    Columns: [iteration, evaluations, sigma, max_std/min_std, std1, std2, ...]
    
    Example access:
    >>> logger = cma.CMADataLogger().load()
    >>> std_devs = logger.std[:, 4:]          # Standard deviations per coordinate
    >>> condition_estimate = logger.std[:, 3] # Max/min std ratio
    """,
    
    'axlen': """
    Principal axis lengths (outcmaesaxlen.dat):
    Columns: [iteration, evaluations, sigma, condition, len1, len2, ...]
    Principal axes ordered by length (largest first).
    
    Example access:
    >>> logger = cma.CMADataLogger().load()
    >>> axis_lengths = logger.D[:, 4:]        # Principal axis lengths
    >>> condition_numbers = logger.D[:, 3]    # Condition numbers
    """,
    
    'axlencorr': """
    Correlation between consecutive axis length vectors (outcmaesaxlencorr.dat):
    Indicates stability of principal directions.
    """
}

# Example: Comprehensive data access
def access_logged_data_example():
    """Example of accessing all types of logged data."""
    
    import cma
    import numpy as np
    
    # Run optimization with logging
    es = cma.CMAEvolutionStrategy(5 * [0.5], 0.3)
    logger = cma.CMADataLogger('example_run/').register(es)
    
    while not es.stop():
        solutions = es.ask()
        fitness_values = [sum((x - np.array([1, 2, 3, 4, 5]))**2) for x in solutions]
        es.tell(solutions, fitness_values)
        logger.add()
        
        if es.countiter % 50 == 0:
            es.disp()
    
    # Reload and analyze data
    logger.load()
    
    # Fitness evolution
    iterations = logger.f[:, 0]
    evaluations = logger.f[:, 1] 
    best_fitness = logger.f[:, 4]
    median_fitness = logger.f[:, 5]
    
    print(f"Optimization completed in {iterations[-1]} iterations")
    print(f"Total evaluations: {evaluations[-1]}")
    print(f"Final best fitness: {best_fitness[-1]:.2e}")
    
    # Solution evolution
    mean_evolution = logger.xmean[:, 4:]
    final_mean = mean_evolution[-1]
    target_solution = np.array([1, 2, 3, 4, 5])
    
    print(f"Final mean solution: {final_mean}")
    print(f"Distance to target: {np.linalg.norm(final_mean - target_solution):.4f}")
    
    # Step-size and conditioning evolution
    sigmas = logger.f[:, 2]
    axis_ratios = logger.f[:, 3]
    
    print(f"Initial sigma: {sigmas[0]:.4f}")
    print(f"Final sigma: {sigmas[-1]:.4f}")
    print(f"Final axis ratio (condition): {axis_ratios[-1]:.2e}")
    
    # Standard deviation evolution per coordinate
    std_evolution = logger.std[:, 4:]
    final_stds = std_evolution[-1]
    
    print(f"Final standard deviations: {final_stds}")
    
    return logger

# access_logged_data_example()

def plot(filenameprefix=None, fig=None, iteridx=None, **kwargs):
    """
    Plot logged CMA-ES data with comprehensive visualization.
    
    Parameters:
    -----------
    filenameprefix : str, optional
        File prefix of logged data. If None, uses default 'outcmaes'.
        
    fig : matplotlib.figure.Figure, optional
        Figure to plot into. If None, creates new figure.
        
    iteridx : slice or array-like, optional
        Iteration indices to plot (e.g., slice(100, 500) for iterations 100-500).
        
    **kwargs : dict
        Additional plotting options:
        - iabscissa : int (0=iteration, 1=evaluation count for x-axis)
        - plot_mean : bool (whether to plot mean solution evolution)
        - foffset : float (offset for fitness values to handle zero/negative values)
        - x_opt : array-like (known optimum for reference)
        - fontsize : int (font size for labels)
        - xsemilog : bool (semi-log x-axis)
        
    Returns:
    --------
    matplotlib.figure.Figure
        Figure with CMA-ES diagnostic plots.
        
    Examples:
    ---------
    >>> import cma
    >>> 
    >>> # Run optimization with default logging
    >>> x, es = cma.fmin2(cma.ff.sphere, [1, 2, 3], 0.5)
    >>> 
    >>> # Plot results
    >>> fig = cma.plot()  # Plot data from default location
    >>> 
    >>> # Plot specific experiment
    >>> fig = cma.plot('my_experiment/run_001')
    >>> 
    >>> # Plot subset of iterations
    >>> fig = cma.plot(iteridx=slice(50, 200))  # Iterations 50-200
    >>> 
    >>> # Plot with evaluation count on x-axis
    >>> fig = cma.plot(iabscissa=1)
    >>> 
    >>> # Plot with known optimum reference
    >>> fig = cma.plot(x_opt=[0, 0, 0])  # Show distance to [0,0,0]
    """
    pass

def plot_zip(filenameprefix=None, **kwargs):
    """
    Plot data from compressed (.zip) CMA-ES log files.
    
    Useful for plotting archived optimization runs without extracting files.
    
    Parameters:
    -----------
    filenameprefix : str, optional
        Prefix of zip file containing logged data.
        
    **kwargs : dict
        Same options as plot() function.
        
    Examples:
    ---------
    >>> import cma
    >>> 
    >>> # Plot from compressed log files
    >>> fig = cma.plot_zip('archived_runs/experiment_001.zip')
    """
    pass

def disp(filenameprefix=None):
    """
    Display text summary of logged CMA-ES data.
    
    Parameters:
    -----------
    filenameprefix : str, optional
        File prefix of logged data to display.
        
    Examples:
    ---------
    >>> import cma
    >>> 
    >>> # Display summary of default logged data
    >>> cma.disp()
    >>> 
    >>> # Display specific experiment
    >>> cma.disp('experiments/run_042')
    """
    pass

class CMADataLogger:
    """Extended plotting capabilities of CMADataLogger."""
    
    def plot(self, fig=None, iabscissa=0, iteridx=None, **kwargs):
        """
        Create comprehensive diagnostic plots for CMA-ES optimization.
        
        Generates multi-panel plot showing:
        1. Fitness evolution (best, median, worst)
        2. Step-size (sigma) evolution  
        3. Principal axis lengths (eigenvalues of C)
        4. Standard deviations per coordinate
        5. Mean solution evolution (optional)
        6. Recent best solution evolution (optional)
        
        Parameters:
        -----------
        fig : matplotlib.figure.Figure, optional
            Figure to plot into. Creates new if None.
            
        iabscissa : int, optional
            X-axis type: 0=iterations (default), 1=evaluation count.
            
        iteridx : slice or array, optional
            Iteration range to plot (e.g., slice(100, 500)).
            
        **kwargs : dict
            Additional options:
            - plot_mean : bool (plot mean evolution, default False)
            - foffset : float (fitness offset for log scale, default 1e-19)
            - x_opt : array (known optimum for reference)
            - fontsize : int (font size, default 7)
            - xsemilog : bool (semi-log x-axis, default False)
            - downsample_to : int (max points to plot, default 1e7)
            
        Examples:
        ---------
        >>> import cma
        >>> 
        >>> # Run optimization
        >>> es = cma.CMAEvolutionStrategy([1, 2, 3], 0.5)
        >>> logger = cma.CMADataLogger().register(es)
        >>> 
        >>> while not es.stop():
        ...     solutions = es.ask()
        ...     fitness_values = [sum(x**2) for x in solutions]
        ...     es.tell(solutions, fitness_values)
        ...     logger.add()
        >>> 
        >>> # Create diagnostic plots
        >>> fig = logger.plot()
        >>> 
        >>> # Plot with custom options
        >>> fig = logger.plot(
        ...     plot_mean=True,           # Show mean evolution
        ...     x_opt=[0, 0, 0],         # Reference optimum
        ...     iabscissa=1,             # X-axis: evaluation count
        ...     iteridx=slice(50, None)  # Skip first 50 iterations
        ... )
        """
        pass
    
    def disp(self):
        """
        Display text summary of logged optimization data.
        
        Shows:
        - Problem dimension and optimization settings
        - Iteration and evaluation count
        - Best, worst, and current fitness values
        - Convergence indicators
        - Step-size and conditioning information
        
        Examples:
        ---------
        >>> logger = cma.CMADataLogger().load()
        >>> logger.disp()
        Loaded data from outcmaes* files
        Dimension: 5
        Iterations: 150
        Evaluations: 1350
        Best fitness: 1.23e-08
        Final sigma: 0.0045
        Condition number: 2.34e+03
        """
        pass

import cma
import numpy as np
import matplotlib.pyplot as plt

# Pattern 1: Real-time plotting during optimization
def realtime_plotting_example():
    """Example of real-time plotting during optimization."""
    
    es = cma.CMAEvolutionStrategy(5 * [0.5], 0.3, {'verb_disp': 0})
    logger = cma.CMADataLogger().register(es)
    
    # Create figure for real-time updates
    plt.ion()  # Interactive mode
    fig, axes = plt.subplots(2, 2, figsize=(10, 8))
    
    iteration = 0
    while not es.stop() and iteration < 200:
        solutions = es.ask()
        fitness_values = [sum((x - 1)**2) for x in solutions]
        es.tell(solutions, fitness_values)
        logger.add()
        
        # Update plots every 10 iterations
        if iteration % 10 == 0:
            plt.clf()  # Clear figure
            
            # Create updated diagnostic plots
            logger.plot(fig=fig)
            plt.draw()
            plt.pause(0.1)  # Brief pause for display
            
        iteration += 1
    
    plt.ioff()  # Turn off interactive mode
    plt.show()
    
    return logger

# Pattern 2: Custom plotting with logged data
def custom_plotting_example():
    """Create custom plots using logged CMA-ES data."""
    
    # Load logged data
    logger = cma.CMADataLogger().load()
    
    # Custom multi-panel plot
    fig, axes = plt.subplots(3, 2, figsize=(12, 10))
    
    # 1. Fitness convergence
    axes[0, 0].semilogy(logger.f[:, 0], logger.f[:, 4], 'b-', label='Best')
    axes[0, 0].semilogy(logger.f[:, 0], logger.f[:, 5], 'r--', label='Median')  
    axes[0, 0].set_xlabel('Iteration')
    axes[0, 0].set_ylabel('Fitness')
    axes[0, 0].legend()
    axes[0, 0].set_title('Fitness Evolution')
    
    # 2. Step-size evolution
    axes[0, 1].semilogy(logger.f[:, 0], logger.f[:, 2])
    axes[0, 1].set_xlabel('Iteration')
    axes[0, 1].set_ylabel('Sigma')
    axes[0, 1].set_title('Step-size Evolution')
    
    # 3. Condition number evolution
    axes[1, 0].semilogy(logger.f[:, 0], logger.f[:, 3])
    axes[1, 0].set_xlabel('Iteration')
    axes[1, 0].set_ylabel('Axis Ratio')
    axes[1, 0].set_title('Condition Number')
    
    # 4. Standard deviations by coordinate
    for i in range(min(5, logger.std.shape[1] - 4)):  # Plot first 5 coordinates
        axes[1, 1].semilogy(logger.std[:, 0], logger.std[:, 4 + i], 
                           label=f'x[{i}]')
    axes[1, 1].set_xlabel('Iteration')
    axes[1, 1].set_ylabel('Standard Deviation')
    axes[1, 1].legend()
    axes[1, 1].set_title('Standard Deviations')
    
    # 5. Principal axis lengths (if available)
    if hasattr(logger, 'D') and logger.D is not None:
        for i in range(min(3, logger.D.shape[1] - 4)):
            axes[2, 0].semilogy(logger.D[:, 0], logger.D[:, 4 + i],
                               label=f'Axis {i+1}')
        axes[2, 0].set_xlabel('Iteration')
        axes[2, 0].set_ylabel('Axis Length')
        axes[2, 0].legend()
        axes[2, 0].set_title('Principal Axis Lengths')
    
    # 6. Solution trajectory (2D projection)
    if logger.xmean.shape[1] >= 6:  # At least 2D problem
        axes[2, 1].plot(logger.xmean[:, 4], logger.xmean[:, 5], 'b-o', markersize=2)
        axes[2, 1].plot(logger.xmean[0, 4], logger.xmean[0, 5], 'go', markersize=8, label='Start')
        axes[2, 1].plot(logger.xmean[-1, 4], logger.xmean[-1, 5], 'ro', markersize=8, label='End')
        axes[2, 1].set_xlabel('x[0]')
        axes[2, 1].set_ylabel('x[1]')
        axes[2, 1].legend()
        axes[2, 1].set_title('Solution Trajectory (2D)')
    
    plt.tight_layout()
    plt.show()
    
    return fig

# Pattern 3: Comparative plotting of multiple runs
def comparative_plotting_example():
    """Compare multiple optimization runs."""
    
    # Simulate multiple runs (in practice, load from different files)
    run_data = {}
    
    for run_id, sigma0 in enumerate([0.1, 0.3, 1.0]):
        logger = cma.CMADataLogger(f'temp_run_{run_id}/')
        es = cma.CMAEvolutionStrategy(3 * [0.5], sigma0, {'verbose': -9})
        logger.register(es)
        
        while not es.stop():
            solutions = es.ask()
            fitness_values = [sum(x**2) for x in solutions]
            es.tell(solutions, fitness_values)
            logger.add()
        
        logger.load()
        run_data[f'sigma_{sigma0}'] = logger
    
    # Comparative plot
    plt.figure(figsize=(12, 8))
    
    # Subplot 1: Fitness comparison
    plt.subplot(2, 2, 1)
    for label, logger in run_data.items():
        plt.semilogy(logger.f[:, 1], logger.f[:, 4], label=label)  # evals vs best fitness
    plt.xlabel('Evaluations')
    plt.ylabel('Best Fitness')
    plt.legend()
    plt.title('Fitness Convergence Comparison')
    
    # Subplot 2: Step-size comparison
    plt.subplot(2, 2, 2)
    for label, logger in run_data.items():
        plt.semilogy(logger.f[:, 0], logger.f[:, 2], label=label)  # iteration vs sigma
    plt.xlabel('Iteration')
    plt.ylabel('Sigma')
    plt.legend()
    plt.title('Step-size Evolution Comparison')
    
    # Subplot 3: Convergence speed
    plt.subplot(2, 2, 3)
    convergence_evals = []
    sigma_values = []
    
    for label, logger in run_data.items():
        # Find evaluation where fitness < 1e-6
        converged_idx = np.where(logger.f[:, 4] < 1e-6)[0]
        if len(converged_idx) > 0:
            convergence_evals.append(logger.f[converged_idx[0], 1])
            sigma_values.append(float(label.split('_')[1]))
        else:
            convergence_evals.append(logger.f[-1, 1])  # Final evaluation
            sigma_values.append(float(label.split('_')[1]))
    
    plt.bar([f'σ={s}' for s in sigma_values], convergence_evals)
    plt.xlabel('Initial Step-size')
    plt.ylabel('Evaluations to Convergence')
    plt.title('Convergence Speed')
    
    # Subplot 4: Final condition numbers
    plt.subplot(2, 2, 4)
    final_conditions = [logger.f[-1, 3] for logger in run_data.values()]
    plt.bar([f'σ={s}' for s in sigma_values], final_conditions)
    plt.xlabel('Initial Step-size')
    plt.ylabel('Final Condition Number')
    plt.yscale('log')
    plt.title('Final Conditioning')
    
    plt.tight_layout()
    plt.show()
    
    return run_data

# Pattern 4: Publication-quality plots
def publication_quality_plots():
    """Create publication-quality plots from CMA-ES data."""
    
    # Load data
    logger = cma.CMADataLogger().load()
    
    # Set publication style
    plt.style.use('classic')
    plt.rcParams.update({
        'font.size': 12,
        'font.family': 'serif',
        'axes.linewidth': 1.5,
        'axes.grid': True,
        'grid.alpha': 0.3,
        'lines.linewidth': 2,
        'figure.dpi': 300
    })
    
    fig, ax = plt.subplots(figsize=(8, 6))
    
    # Plot fitness evolution with confidence bands (if multiple runs available)
    iterations = logger.f[:, 0]
    best_fitness = logger.f[:, 4]
    median_fitness = logger.f[:, 5]
    
    ax.semilogy(iterations, best_fitness, 'b-', linewidth=2, label='Best fitness')
    ax.semilogy(iterations, median_fitness, 'r--', linewidth=2, label='Median fitness')
    
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Function value')
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_title('CMA-ES Convergence on Sphere Function')
    
    # Add annotations
    if len(best_fitness) > 0:
        # Mark significant milestones
        target_values = [1e-2, 1e-6, 1e-10]
        for target in target_values:
            achieved_idx = np.where(best_fitness <= target)[0]
            if len(achieved_idx) > 0:
                iter_achieved = iterations[achieved_idx[0]]
                ax.annotate(f'10^{{{int(np.log10(target))}}} at iter {iter_achieved}',
                           xy=(iter_achieved, target),
                           xytext=(iter_achieved + 20, target * 10),
                           arrowprops=dict(arrowstyle='->', color='gray'))
    
    plt.tight_layout()
    plt.savefig('cmaes_convergence.pdf', dpi=300, bbox_inches='tight')
    plt.show()
    
    return fig

# Run plotting examples (comment out as needed)
# realtime_plotting_example()
# custom_plotting_example()
# comparative_plotting_example()
# publication_quality_plots()

def convergence_analysis(logger_or_prefix=None):
    """
    Comprehensive convergence analysis of CMA-ES optimization.
    
    Parameters:
    -----------
    logger_or_prefix : CMADataLogger or str, optional
        Logger instance or file prefix. If None, loads default data.
        
    Returns:
    --------
    dict
        Analysis results including convergence rates, bottlenecks, and diagnostics.
        
    Examples:
    ---------
    >>> import cma
    >>> 
    >>> # Run optimization
    >>> x, es = cma.fmin2(cma.ff.elli, [1, 2, 3], 0.5)
    >>> 
    >>> # Analyze convergence
    >>> analysis = convergence_analysis()
    >>> 
    >>> print(f"Linear convergence rate: {analysis['linear_rate']:.3f}")
    >>> print(f"Bottleneck phase: {analysis['bottleneck_phase']}")
    >>> print(f"Final condition: {analysis['final_condition']:.2e}")
    """
    
    # Load data if needed
    if isinstance(logger_or_prefix, str):
        logger = cma.CMADataLogger().load(logger_or_prefix)
    elif logger_or_prefix is None:
        logger = cma.CMADataLogger().load()
    else:
        logger = logger_or_prefix
    
    analysis = {}
    
    # Basic statistics
    analysis['total_iterations'] = int(logger.f[-1, 0])
    analysis['total_evaluations'] = int(logger.f[-1, 1])
    analysis['initial_fitness'] = logger.f[0, 4]
    analysis['final_fitness'] = logger.f[-1, 4]
    analysis['fitness_improvement'] = np.log10(analysis['initial_fitness'] / analysis['final_fitness'])
    
    # Convergence rate analysis
    log_fitness = np.log(logger.f[:, 4])
    iterations = logger.f[:, 0]
    
    # Linear convergence rate (slope of log(fitness) vs iteration)
    if len(log_fitness) > 10:
        # Use latter half for stable convergence rate
        start_idx = len(log_fitness) // 2
        coeffs = np.polyfit(iterations[start_idx:], log_fitness[start_idx:], 1)
        analysis['linear_rate'] = -coeffs[0]  # Negative because we want decrease
    else:
        analysis['linear_rate'] = 0
    
    # Identify convergence phases
    fitness_diff = np.diff(log_fitness)
    
    # Find where improvement slows down significantly  
    slow_improvement = fitness_diff > -1e-3  # Very slow improvement
    if np.any(slow_improvement):
        analysis['bottleneck_start'] = int(iterations[np.where(slow_improvement)[0][0]])
    else:
        analysis['bottleneck_start'] = None
    
    # Step-size analysis
    sigmas = logger.f[:, 2]
    analysis['initial_sigma'] = sigmas[0]
    analysis['final_sigma'] = sigmas[-1]
    analysis['sigma_reduction'] = sigmas[0] / sigmas[-1]
    
    # Conditioning analysis
    axis_ratios = logger.f[:, 3]
    analysis['initial_condition'] = axis_ratios[0]
    analysis['final_condition'] = axis_ratios[-1]
    analysis['max_condition'] = np.max(axis_ratios)
    
    # Stagnation detection
    fitness_window = 20  # Window for stagnation check
    if len(logger.f) > fitness_window:
        recent_fitness = logger.f[-fitness_window:, 4]
        fitness_variation = (np.max(recent_fitness) - np.min(recent_fitness)) / np.mean(recent_fitness)
        analysis['recent_stagnation'] = fitness_variation < 1e-12
    else:
        analysis['recent_stagnation'] = False
    
    # Efficiency metrics
    analysis['evaluations_per_iteration'] = analysis['total_evaluations'] / analysis['total_iterations']
    if analysis['fitness_improvement'] > 0:
        analysis['evals_per_order_magnitude'] = analysis['total_evaluations'] / analysis['fitness_improvement']
    else:
        analysis['evals_per_order_magnitude'] = float('inf')
    
    return analysis

def diagnostic_summary(logger_or_prefix=None):
    """
    Generate diagnostic summary for CMA-ES optimization.
    
    Identifies common issues and provides recommendations.
    """
    
    analysis = convergence_analysis(logger_or_prefix)
    
    print("CMA-ES Diagnostic Summary")
    print("=" * 40)
    
    # Basic performance
    print(f"Total iterations: {analysis['total_iterations']}")
    print(f"Total evaluations: {analysis['total_evaluations']}")
    print(f"Fitness improvement: {analysis['fitness_improvement']:.1f} orders of magnitude")
    print(f"Convergence rate: {analysis['linear_rate']:.3f} (log units per iteration)")
    
    # Efficiency assessment
    print(f"\nEfficiency Metrics:")
    print(f"  Evaluations per iteration: {analysis['evaluations_per_iteration']:.1f}")
    if analysis['evals_per_order_magnitude'] < float('inf'):
        print(f"  Evaluations per order of magnitude: {analysis['evals_per_order_magnitude']:.0f}")
    
    # Conditioning assessment
    print(f"\nConditioning Analysis:")
    print(f"  Initial condition number: {analysis['initial_condition']:.2e}")
    print(f"  Final condition number: {analysis['final_condition']:.2e}")
    print(f"  Maximum condition number: {analysis['max_condition']:.2e}")
    
    # Issue detection and recommendations
    print(f"\nDiagnostic Findings:")
    
    issues = []
    recommendations = []
    
    # Check for poor conditioning
    if analysis['final_condition'] > 1e12:
        issues.append("Very high condition number (> 1e12)")
        recommendations.append("Consider coordinate scaling or preconditioning")
    
    # Check for stagnation
    if analysis['recent_stagnation']:
        issues.append("Recent fitness stagnation detected")
        recommendations.append("May need restart or different termination criteria")
    
    # Check convergence rate
    if analysis['linear_rate'] < 0.01:
        issues.append("Slow convergence rate (< 0.01)")
        recommendations.append("Consider larger population or different step-size")
    
    # Check step-size reduction
    if analysis['sigma_reduction'] > 1e6:
        issues.append("Excessive step-size reduction")
        recommendations.append("Initial step-size may have been too large")
    
    # Check efficiency
    expected_evals = analysis['total_iterations'] * (4 + 3 * np.log(5))  # Rough estimate for 5D
    if analysis['total_evaluations'] > 2 * expected_evals:
        issues.append("Higher than expected evaluation count")
        recommendations.append("Check for expensive function evaluations or large population")
    
    if not issues:
        print("  No significant issues detected - optimization appears healthy")
    else:
        for i, issue in enumerate(issues):
            print(f"  Issue {i+1}: {issue}")
        
        print(f"\nRecommendations:")
        for i, rec in enumerate(recommendations):
            print(f"  {i+1}. {rec}")
    
    return analysis

# Example usage
def run_diagnostic_example():
    """Example of comprehensive diagnostic analysis."""
    
    import cma
    
    # Run optimization with logging
    def objective(x):
        # Ill-conditioned ellipsoid 
        scales = [10**(i/2) for i in range(len(x))]
        return sum(s * xi**2 for s, xi in zip(scales, x))
    
    x, es = cma.fmin2(objective, 5 * [1], 0.5, 
                      options={'maxfevals': 5000, 'verbose': -1})
    
    # Run diagnostics
    analysis = diagnostic_summary()
    
    # Custom analysis
    print(f"\nCustom Analysis:")
    logger = cma.CMADataLogger().load()
    
    # Check for premature convergence
    final_sigma = logger.f[-1, 2]
    dimension = 5
    if final_sigma < 1e-12:
        print("  Warning: Very small final step-size - possible premature convergence")
    
    # Check population diversity
    if hasattr(logger, 'std') and logger.std is not None:
        final_stds = logger.std[-1, 4:]
        std_ratio = np.max(final_stds) / np.min(final_stds)
        print(f"  Final coordinate std ratio: {std_ratio:.2e}")
        
        if std_ratio > 1e6:
            print("  Warning: Very different coordinate scales - consider rescaling")
    
    return analysis

# run_diagnostic_example()

def performance_metrics(results_list):
    """
    Compute standard performance metrics for optimization results.
    
    Parameters:
    -----------
    results_list : list
        List of (x_best, es) tuples from multiple optimization runs.
        
    Returns:
    --------
    dict
        Performance metrics including success rates, efficiency measures.
    """
    
    metrics = {
        'success_rate': 0,
        'mean_evaluations': 0,
        'median_evaluations': 0,
        'mean_final_fitness': 0,
        'evaluation_counts': [],
        'final_fitness_values': [],
        'successful_runs': 0,
        'total_runs': len(results_list)
    }
    
    success_threshold = 1e-8
    
    for x_best, es in results_list:
        final_fitness = es.result.fbest
        evaluations = es.result.evaluations
        
        metrics['evaluation_counts'].append(evaluations)
        metrics['final_fitness_values'].append(final_fitness)
        
        if final_fitness < success_threshold:
            metrics['successful_runs'] += 1
    
    # Compute derived metrics
    metrics['success_rate'] = metrics['successful_runs'] / metrics['total_runs']
    metrics['mean_evaluations'] = np.mean(metrics['evaluation_counts'])
    metrics['median_evaluations'] = np.median(metrics['evaluation_counts'])
    metrics['mean_final_fitness'] = np.mean(metrics['final_fitness_values'])
    
    # Success-conditional metrics
    successful_evals = [evals for evals, (_, es) in zip(metrics['evaluation_counts'], results_list)
                       if es.result.fbest < success_threshold]
    
    if successful_evals:
        metrics['mean_successful_evaluations'] = np.mean(successful_evals)
        metrics['median_successful_evaluations'] = np.median(successful_evals)
    else:
        metrics['mean_successful_evaluations'] = None
        metrics['median_successful_evaluations'] = None
    
    return metrics

def benchmark_suite_analysis():
    """Run CMA-ES on standard benchmark suite and analyze performance."""
    
    import cma
    
    # Test functions
    test_functions = {
        'sphere': lambda x: sum(x**2),
        'ellipsoid': lambda x: sum(10**(6*i/(len(x)-1)) * xi**2 for i, xi in enumerate(x)),
        'rosenbrock': lambda x: sum(100*(x[i+1] - x[i]**2)**2 + (1 - x[i])**2 for i in range(len(x)-1)),
        'rastrigin': lambda x: sum(xi**2 - 10*np.cos(2*np.pi*xi) + 10 for xi in x)
    }
    
    dimensions = [2, 5, 10]
    repetitions = 5
    
    results = {}
    
    for func_name, func in test_functions.items():
        results[func_name] = {}
        
        for dim in dimensions:
            print(f"Testing {func_name} in {dim}D...")
            
            run_results = []
            
            for rep in range(repetitions):
                try:
                    x, es = cma.fmin2(func, dim * [0.5], 0.3,
                                    options={'maxfevals': 1000 * dim**2, 'verbose': -9})
                    run_results.append((x, es))
                except Exception as e:
                    print(f"  Run {rep} failed: {e}")
            
            if run_results:
                metrics = performance_metrics(run_results)
                results[func_name][dim] = metrics
    
    # Print summary
    print("\nBenchmark Results Summary")
    print("=" * 50)
    
    for func_name in results:
        print(f"\nFunction: {func_name}")
        print("Dim  Success%   Mean Evals   Median Evals   Mean Final")
        print("-" * 55)
        
        for dim in sorted(results[func_name].keys()):
            m = results[func_name][dim]
            print(f"{dim:3d}    {m['success_rate']:5.1%}      {m['mean_evaluations']:8.0f}     "
                  f"{m['median_evaluations']:8.0f}     {m['mean_final_fitness']:8.2e}")
    
    return results

# benchmark_suite_analysis()