GPopt.GeneralizationOpt

  1import numpy as np
  2import pandas as pd
  3from tqdm import tqdm
  4from scipy.stats.qmc import Sobol
  5from scipy.optimize import minimize as scipy_minimize
  6from scipy.stats import norm
  7from sklearn.model_selection import train_test_split
  8from sklearn.base import BaseEstimator, RegressorMixin
  9from SALib.analyze import sobol as sobol_analyze
 10from SALib.sample import sobol as sobol_sample
 11import warnings
 12
 13warnings.filterwarnings("ignore")
 14
 15
 16class GenericSurrogate(BaseEstimator, RegressorMixin):
 17    """
 18    Generic wrapper for any surrogate model with conformal prediction support
 19
 20    This class wraps any sklearn-compatible model and adds:
 21    - Conformal prediction intervals (distribution-free uncertainty)
 22    - Optional Bayesian uncertainty (if model supports it)
 23    - Unified prediction interface
 24
 25    Parameters:
 26    -----------
 27    base_model : object
 28        Any sklearn-compatible model with fit() and predict() methods
 29    conformal : bool
 30        Whether to use conformal prediction for uncertainty
 31    bayesian : bool
 32        Whether model provides native uncertainty (e.g., GP, Bayesian NN)
 33    """
 34
 35    def __init__(self, base_model, conformal=True, bayesian=False):
 36        self.base_model = base_model
 37        self.conformal = conformal
 38        self.bayesian = bayesian
 39        self.nonconformity_scores_ = None
 40        self.is_fitted_ = False
 41
 42    def fit(self, X, y, cal_ratio=0.2):
 43        """
 44        Fit surrogate with optional conformal calibration
 45
 46        Parameters:
 47        -----------
 48        X : array-like, shape (n_samples, n_features)
 49            Training features
 50        y : array-like, shape (n_samples,)
 51            Training targets
 52        cal_ratio : float
 53            Proportion of data for conformal calibration (if conformal=True)
 54
 55        Returns:
 56        --------
 57        self : object
 58        """
 59        if self.conformal and cal_ratio > 0:
 60            # Split conformal prediction: reserve data for calibration
 61            n_cal = int(len(X) * cal_ratio)
 62            X_train, X_cal = X[:-n_cal], X[-n_cal:]
 63            y_train, y_cal = y[:-n_cal], y[-n_cal:]
 64
 65            self.base_model.fit(X_train, y_train)
 66
 67            # Compute nonconformity scores on calibration set
 68            y_cal_pred = self.base_model.predict(X_cal)
 69            self.nonconformity_scores_ = np.abs(y_cal - y_cal_pred)
 70        else:
 71            # No conformal prediction: use all data
 72            self.base_model.fit(X, y)
 73
 74        self.is_fitted_ = True
 75        return self
 76
 77    def predict(self, X, return_std=False, confidence=0.95):
 78        """
 79        Unified prediction interface with optional uncertainty
 80
 81        Parameters:
 82        -----------
 83        X : array-like, shape (n_samples, n_features)
 84            Test features
 85        return_std : bool
 86            Whether to return uncertainty estimates
 87        confidence : float
 88            Confidence level for prediction intervals (if conformal=True)
 89
 90        Returns:
 91        --------
 92        predictions : array or tuple
 93            If return_std=False: predictions only
 94            If return_std=True: (predictions, lower_bound, upper_bound)
 95
 96        Notes:
 97        ------
 98        - Conformal: Distribution-free intervals with finite-sample guarantees
 99        - Bayesian: Model-based uncertainty (e.g., GP posterior std)
100        """
101        if not self.is_fitted_:
102            raise ValueError("Model must be fitted before prediction")
103
104        # Get point predictions
105        y_pred = self.base_model.predict(X)
106
107        if not return_std:
108            return y_pred
109
110        # Compute uncertainty estimates
111        if self.conformal and self.nonconformity_scores_ is not None:
112            # Conformal prediction intervals
113            alpha = 1 - confidence
114            n_scores = len(self.nonconformity_scores_)
115
116            # Calculate the effective quantile level based on (n+1) rule
117            # Cap the quantile level at 1.0 to prevent ValueError from np.quantile
118            quantile_level = min(
119                1.0, np.ceil((1 - alpha) * (n_scores + 1)) / n_scores
120            )
121
122            q = np.quantile(self.nonconformity_scores_, quantile_level)
123            y_lower = y_pred - q
124            y_upper = y_pred + q
125
126        elif self.bayesian and hasattr(self.base_model, "predict"):
127            # Try to get Bayesian uncertainty from model
128            try:
129                # For GP-like models
130                if (
131                    hasattr(self.base_model, "predict")
132                    and "return_std"
133                    in self.base_model.predict.__code__.co_varnames
134                ):
135                    _, y_std = self.base_model.predict(X, return_std=True)
136                    y_lower = y_pred - 2 * y_std  # 95% CI for Gaussian
137                    y_upper = y_pred + 2 * y_std
138                else:
139                    # Fallback: assume constant uncertainty
140                    y_std = np.std(y_pred) * np.ones_like(y_pred)
141                    y_lower = y_pred - 2 * y_std
142                    y_upper = y_pred + 2 * y_std
143            except:
144                # No uncertainty available
145                y_lower = y_pred
146                y_upper = y_pred
147        else:
148            # No uncertainty quantification available
149            y_lower = y_pred
150            y_upper = y_pred
151
152        return y_pred, y_lower, y_upper
153
154
155class GeneralizationOpt:
156    """
157    Generic framework for model generalization diagnostics
158
159    This class is agnostic to:
160    - The base model being tuned (LightGBM, XGBoost, neural nets, etc.)
161    - The surrogate model used (GP, RF, neural nets, etc.)
162    - The hyperparameter space structure
163
164    Example usage:
165    --------------
166    >>> from sklearn.gaussian_process import GaussianProcessRegressor
167    >>> from sklearn.ensemble import RandomForestRegressor
168    >>>
169    >>> # Define hyperparameter space
170    >>> hyperparams = {
171    >>>     'learning_rate': (0.01, 0.3, True),  # (min, max, log_scale)
172    >>>     'max_depth': (3, 10, False)
173    >>> }
174    >>>
175    >>> # Initialize framework
176    >>> diag = GeneralizationOpt(hyperparams, random_state=42)
177    >>>
178    >>> # Use GP surrogate with Bayesian uncertainty
179    >>> gp = GaussianProcessRegressor()
180    >>> surrogate = GenericSurrogate(gp, conformal=False, bayesian=True)
181    >>>
182    >>> # Or use RF with conformal prediction
183    >>> rf = RandomForestRegressor()
184    >>> surrogate = GenericSurrogate(rf, conformal=True, bayesian=False)
185    >>>
186    >>> # Run diagnostics
187    >>> results = diag.run_diagnostics(
188    >>>     X_train, X_test, y_train, y_test,
189    >>>     model_class=LGBMClassifier,
190    >>>     surrogate=surrogate,
191    >>>     n_samples=100
192    >>> )
193    """
194
195    def __init__(self, hyperparams, random_state=42):
196        """
197        Initialize diagnostics framework
198
199        Parameters:
200        -----------
201        hyperparams : dict
202            Hyperparameter space definition
203            Format: {param_name: (min_val, max_val, log_scale), ...}
204        random_state : int
205            Random seed for reproducibility
206        """
207        self.hyperparams = hyperparams
208        self.hyperparam_names = list(hyperparams.keys())
209        self.n_dims = len(self.hyperparam_names)
210        self.random_state = random_state
211
212        # Storage
213        self.sobol_samples = None
214        self.evaluation_results = None
215        self.surrogate = None
216        self.sobol_indices = None
217
218    def _denormalize_config(self, normalized_config):
219        """Convert normalized [0,1] to actual hyperparameter values"""
220        config = {}
221        for i, param_name in enumerate(self.hyperparam_names):
222            min_val, max_val, log_scale = self.hyperparams[param_name]
223            norm_val = normalized_config[i]
224
225            if log_scale:
226                log_min, log_max = np.log(min_val), np.log(max_val)
227                config[param_name] = np.exp(
228                    log_min + norm_val * (log_max - log_min)
229                )
230            else:
231                denormalized_value = min_val + norm_val * (max_val - min_val)
232                # Special handling for boolean parameters represented as (0, 1, False)
233                if (
234                    min_val == 0
235                    and max_val == 1
236                    and not log_scale
237                    and isinstance(min_val, int)
238                ):
239                    config[param_name] = bool(round(denormalized_value))
240                # Convert to int if needed (if min and max are integers and not log-scale)
241                elif (
242                    isinstance(min_val, int)
243                    and isinstance(max_val, int)
244                    and not log_scale
245                ):
246                    config[param_name] = int(round(denormalized_value))
247                else:
248                    config[param_name] = denormalized_value
249
250        return config
251
252    def generate_sobol_samples(self, n_samples=None):
253        """
254        Generate Sobol quasi-random samples in hyperparameter space
255
256        Parameters:
257        -----------
258        n_samples : int, optional
259            Number of samples (default: 75 * n_dims)
260
261        Returns:
262        --------
263        samples : ndarray, shape (n_samples, n_dims)
264            Normalized samples in [0, 1]^n_dims
265        """
266        if n_samples is None:
267            n_samples = 75 * self.n_dims
268
269        sobol = Sobol(d=self.n_dims, scramble=True, seed=self.random_state)
270        self.sobol_samples = sobol.random(n_samples)
271
272        return self.sobol_samples
273
274    def evaluate_configurations(
275        self,
276        X_train,
277        X_test,
278        y_train,
279        y_test,
280        model_class,
281        model_kwargs=None,
282        evaluation_fn=None,
283        n_folds=5,
284    ):
285        """
286        Evaluate model configurations and compute generalization gaps
287
288        Parameters:
289        -----------
290        X_train, X_test, y_train, y_test : arrays
291            Train/test data
292        model_class : class
293            Model class to instantiate (e.g., LGBMClassifier)
294        model_kwargs : dict, optional
295            Fixed kwargs for model (e.g., {'random_state': 42, 'verbose': -1})
296        evaluation_fn : callable, optional
297            Function to evaluate model: evaluation_fn(y_true, y_pred) -> score
298            Default: balanced_accuracy_score for classification
299        n_folds : int
300            Number of CV folds for training score estimation
301
302        Returns:
303        --------
304        results : DataFrame
305            Results with columns: hyperparams, cv_train, test_score, gap_abs, etc.
306        """
307        from sklearn.model_selection import StratifiedKFold, KFold
308        from sklearn.metrics import balanced_accuracy_score, r2_score
309
310        if self.sobol_samples is None:
311            raise ValueError("Must generate Sobol samples first")
312
313        if model_kwargs is None:
314            model_kwargs = {}
315
316        if evaluation_fn is None:
317            # Heuristic: check if classification or regression
318            if len(np.unique(y_train)) < 20:
319                evaluation_fn = balanced_accuracy_score
320                kfold = StratifiedKFold(
321                    n_splits=n_folds,
322                    shuffle=True,
323                    random_state=self.random_state,
324                )
325            else:
326                evaluation_fn = r2_score
327                kfold = KFold(
328                    n_splits=n_folds,
329                    shuffle=True,
330                    random_state=self.random_state,
331                )
332        else:
333            kfold = KFold(
334                n_splits=n_folds, shuffle=True, random_state=self.random_state
335            )
336
337        results = []
338        n_samples = len(self.sobol_samples)
339
340        print(
341            f"Evaluating {n_samples} configurations ({n_folds}-fold CV + test)..."
342        )
343
344        for idx, norm_config in tqdm(enumerate(self.sobol_samples)):
345            if (idx + 1) % 50 == 0:
346                print(
347                    f"  Progress: {idx + 1}/{n_samples} "
348                    f"({100*(idx+1)/n_samples:.1f}%)"
349                )
350
351            config = self._denormalize_config(norm_config)
352
353            # Cross-validation score
354            cv_scores = []
355            for train_idx, val_idx in kfold.split(X_train, y_train):
356                if hasattr(X_train, "iloc"):  # DataFrame
357                    X_tr, X_val = X_train.iloc[train_idx], X_train.iloc[val_idx]
358                    y_tr, y_val = y_train.iloc[train_idx], y_train.iloc[val_idx]
359                else:  # ndarray
360                    X_tr, X_val = X_train[train_idx], X_train[val_idx]
361                    y_tr, y_val = y_train[train_idx], y_train[val_idx]
362
363                model = model_class(**config, **model_kwargs)
364                model.fit(X_tr, y_tr)
365                y_pred = model.predict(X_val)
366                cv_scores.append(evaluation_fn(y_val, y_pred))
367
368            cv_train = np.mean(cv_scores)
369
370            # Test score
371            model = model_class(**config, **model_kwargs)
372            model.fit(X_train, y_train)
373            y_pred_test = model.predict(X_test)
374            test_score = evaluation_fn(y_test, y_pred_test)
375
376            # Generalization gaps
377            gap_abs = cv_train - test_score
378            gap_norm = (
379                (cv_train - test_score) / (1 - cv_train) if cv_train < 1 else 0
380            )
381            gap_mpe = (
382                ((cv_train - test_score) / cv_train * 100)
383                if cv_train > 0
384                else 0
385            )
386
387            results.append(
388                {
389                    **config,
390                    "cv_train": cv_train,
391                    "test_score": test_score,
392                    "gap_abs": gap_abs,
393                    "gap_norm": gap_norm,
394                    "gap_mpe": gap_mpe,
395                }
396            )
397
398        self.evaluation_results = pd.DataFrame(results)
399
400        print(f"\n✓ Completed {n_samples} evaluations")
401        print(
402            f"  CV Train Score: {self.evaluation_results['cv_train'].mean():.4f} \u00b1 "
403            f"{self.evaluation_results['cv_train'].std():.4f}"
404        )
405        print(
406            f"  Test Score: {self.evaluation_results['test_score'].mean():.4f} \u00b1 "
407            f"{self.evaluation_results['test_score'].std():.4f}"
408        )
409        print(
410            f"  Absolute Gap: {self.evaluation_results['gap_abs'].mean():.4f} \u00b1 "
411            f"{self.evaluation_results['gap_abs'].std():.4f}"
412        )
413
414        return self.evaluation_results
415
416    def fit_surrogate(self, surrogate, target_gap="gap_abs"):
417        """
418        Fit surrogate model to generalization gap data
419
420        Parameters:
421        -----------
422        surrogate : GenericSurrogate
423            Surrogate model instance
424        target_gap : str
425            Target variable to model ('gap_abs', 'gap_norm', 'gap_mpe')
426
427        Returns:
428        --------
429        surrogate : GenericSurrogate
430            Fitted surrogate model
431        metrics : dict
432            Performance metrics (R², MAE, coverage if applicable)
433        """
434        if self.evaluation_results is None:
435            raise ValueError("Must evaluate configurations first")
436
437        X = self.sobol_samples
438        y = self.evaluation_results[target_gap].values
439
440        # Train/validation split
441        split_idx = int(0.8 * len(X))
442        X_train, X_val = X[:split_idx], X[split_idx:]
443        y_train, y_val = y[:split_idx], y[split_idx:]
444
445        print(
446            f"Training surrogate on {len(X_train)} samples, "
447            f"validating on {len(X_val)} samples..."
448        )
449
450        # Fit surrogate
451        surrogate.fit(X_train, y_train)
452
453        # Evaluate
454        y_pred, y_lower, y_upper = surrogate.predict(X_val, return_std=True)
455
456        from sklearn.metrics import r2_score, mean_absolute_error
457
458        r2 = r2_score(y_val, y_pred)
459        mae = mean_absolute_error(y_val, y_pred)
460
461        metrics = {"r2": r2, "mae": mae}
462
463        # Compute coverage if intervals available
464        if not np.array_equal(y_lower, y_upper):
465            coverage = np.mean((y_val >= y_lower) & (y_val <= y_upper))
466            interval_width = np.mean(y_upper - y_lower)
467            metrics["coverage"] = coverage
468            metrics["interval_width"] = interval_width
469
470            print(f"\n✓ Surrogate Performance:")
471            print(f"  R² Score: {r2:.4f}")
472            print(f"  MAE: {mae:.6f}")
473            print(f"  Coverage: {coverage:.2%}")
474            print(f"  Mean Interval Width: {interval_width:.6f}")
475        else:
476            print(f"\n✓ Surrogate Performance:")
477            print(f"  R² Score: {r2:.4f}")
478            print(f"  MAE: {mae:.6f}")
479
480        self.surrogate = surrogate
481        return surrogate, metrics
482
483    def compute_sobol_indices(self, n_samples=10000):
484        """
485        Compute Sobol sensitivity indices for global sensitivity analysis
486
487        Parameters:
488        -----------
489        n_samples : int
490            Number of Monte Carlo samples for Sobol estimation
491
492        Returns:
493        --------
494        indices : DataFrame
495            Sobol indices with columns: Parameter, S1, ST, Interaction
496        """
497        if self.surrogate is None:
498            raise ValueError("Must fit surrogate first")
499
500        problem = {
501            "num_vars": self.n_dims,
502            "names": self.hyperparam_names,
503            "bounds": [[0, 1]] * self.n_dims,
504        }
505
506        param_values = sobol_sample.sample(
507            problem, n_samples, calc_second_order=False
508        )
509        Y = self.surrogate.predict(param_values)
510
511        Si = sobol_analyze.analyze(
512            problem, Y, calc_second_order=False, print_to_console=False
513        )
514
515        self.sobol_indices = pd.DataFrame(
516            {
517                "Parameter": self.hyperparam_names,
518                "S1": Si["S1"],
519                "ST": Si["ST"],
520                "Interaction": Si["ST"] - Si["S1"],
521            }
522        ).sort_values("S1", ascending=False)
523
524        print("\n" + "=" * 70)
525        print("Sobol Sensitivity Indices (Variance Decomposition)")
526        print("=" * 70)
527        print(self.sobol_indices.to_string(index=False))
528        print("\nS1 = First-order effect (main effect)")
529        print("ST = Total-order effect (main + interactions)")
530        print("Interaction = ST - S1 (pure interaction effects)")
531
532        return self.sobol_indices
533
534    def optimize_acquisition(
535        self,
536        acquisition="ei",
537        n_starts=10,
538        kappa=2.0,
539        xi=0.01,
540        minimize_gap=True,
541    ):
542        """
543        Optimize acquisition function to find best hyperparameter configuration
544
545        Parameters:
546        -----------
547        acquisition : str
548            Acquisition function: 'ei' (Expected Improvement),
549            'ucb' (Upper Confidence Bound), 'pi' (Probability of Improvement)
550        n_starts : int
551            Number of random restarts for multi-start optimization
552        kappa : float
553            Exploration parameter for UCB (higher = more exploration)
554        xi : float
555            Exploration parameter for EI/PI (higher = more exploration)
556        minimize_gap : bool
557            If True, minimize gap; if False, maximize test performance
558
559        Returns:
560        --------
561        results : dict
562            Optimization results including best config and expected improvement
563        """
564        if self.surrogate is None:
565            raise ValueError("Must fit surrogate first")
566
567        # Get current best
568        if minimize_gap:
569            y_best = self.evaluation_results["gap_abs"].min()
570            print(f"Current best gap: {y_best:.6f}")
571        else:
572            y_best = self.evaluation_results["test_score"].max()
573            print(f"Current best test score: {y_best:.6f}")
574
575        # Define acquisition functions
576        def expected_improvement(x):
577            x = x.reshape(1, -1)
578            y_pred, y_lower, y_upper = self.surrogate.predict(
579                x, return_std=True
580            )
581            mu, sigma = y_pred[0], (y_upper[0] - y_lower[0]) / 4  # Approx std
582
583            if sigma < 1e-10:
584                return 0.0
585
586            if minimize_gap:
587                imp = y_best - mu - xi
588            else:
589                imp = mu - y_best - xi
590
591            Z = imp / sigma
592            ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
593            return -ei
594
595        def upper_confidence_bound(x):
596            x = x.reshape(1, -1)
597            y_pred, y_lower, y_upper = self.surrogate.predict(
598                x, return_std=True
599            )
600            mu, sigma = y_pred[0], (y_upper[0] - y_lower[0]) / 4
601
602            if minimize_gap:
603                ucb = mu - kappa * sigma
604            else:
605                ucb = mu + kappa * sigma
606
607            return ucb if minimize_gap else -ucb
608
609        def probability_of_improvement(x):
610            x = x.reshape(1, -1)
611            y_pred, y_lower, y_upper = self.surrogate.predict(
612                x, return_std=True
613            )
614            mu, sigma = y_pred[0], (y_upper[0] - y_lower[0]) / 4
615
616            if sigma < 1e-10:
617                return 0.0
618
619            if minimize_gap:
620                Z = (y_best - mu - xi) / sigma
621            else:
622                Z = (mu - y_best - xi) / sigma
623
624            pi = norm.cdf(Z)
625            return -pi
626
627        acq_funcs = {
628            "ei": expected_improvement,
629            "ucb": upper_confidence_bound,
630            "pi": probability_of_improvement,
631        }
632        acq_func = acq_funcs[acquisition]
633
634        # Multi-start optimization
635        print(f"\nRunning {n_starts}-start L-BFGS-B optimization...")
636
637        best_result = None
638        best_acq_value = np.inf
639
640        np.random.seed(self.random_state)
641        starting_points = np.random.rand(n_starts, self.n_dims)
642
643        # Add best observed point as starting point
644        if len(self.sobol_samples) > 0:
645            best_idx = (
646                self.evaluation_results["gap_abs"].idxmin()
647                if minimize_gap
648                else self.evaluation_results["test_score"].idxmax()
649            )
650            starting_points[0] = self.sobol_samples[best_idx]
651
652        for i, x0 in enumerate(starting_points):
653            try:
654                result = scipy_minimize(
655                    acq_func,
656                    x0,
657                    method="L-BFGS-B",
658                    bounds=[(0, 1)] * self.n_dims,
659                    options={"maxiter": 1000, "ftol": 1e-9},
660                )
661
662                if result.success and result.fun < best_acq_value:
663                    best_acq_value = result.fun
664                    best_result = result
665            except Exception as e:
666                print(f"  Warning: Optimization {i+1} failed: {str(e)}")
667
668        if best_result is None:
669            raise RuntimeError("All optimization attempts failed")
670
671        # Extract best configuration
672        best_x_normalized = best_result.x
673        best_config = self._denormalize_config(best_x_normalized)
674
675        # Predict gap at optimum
676        y_pred, y_lower, y_upper = self.surrogate.predict(
677            best_x_normalized.reshape(1, -1), return_std=True
678        )
679        mu_best = y_pred[0]
680
681        print(f"\n✓ Optimization Complete")
682        print(f"  Acquisition: {acquisition.upper()}")
683        print(f"  Best acquisition value: {-best_acq_value:.6f}")
684        print(f"  Predicted gap at optimum: {mu_best:.6f}")
685        print(f"  95% CI: [{y_lower[0]:.6f}, {y_upper[0]:.6f}]")
686        print(f"\n  Recommended Configuration:")
687        for param, value in best_config.items():
688            print(
689                f"    {param}: {value if isinstance(value, int) else f'{value:.4f}'}"
690            )
691
692        improvement = y_best - mu_best if minimize_gap else mu_best - y_best
693        print(f"\n  Expected improvement: {improvement:.6f}")
694
695        return {
696            "best_config": best_config,
697            "best_config_normalized": best_x_normalized,
698            "predicted_gap_mean": mu_best,
699            "predicted_gap_lower": y_lower[0],
700            "predicted_gap_upper": y_upper[0],
701            "current_best": y_best,
702            "expected_improvement": improvement,
703            "acquisition_function": acquisition,
704        }
705
706    def run_diagnostics(
707        self,
708        X_train,
709        X_test,
710        y_train,
711        y_test,
712        model_class,
713        surrogate,
714        n_samples=None,
715        model_kwargs=None,
716        target_gap="gap_abs",
717    ):
718        """
719        Run complete diagnostic pipeline
720
721        Parameters:
722        -----------
723        X_train, X_test, y_train, y_test : arrays
724            Train/test data
725        model_class : class
726            Model class to tune
727        surrogate : GenericSurrogate
728            Surrogate model for gap modeling
729        n_samples : int, optional
730            Number of Sobol samples
731        model_kwargs : dict, optional
732            Fixed model kwargs
733        target_gap : str
734            Gap metric to model
735
736        Returns:
737        --------
738        results : dict
739            Complete diagnostic results
740        """
741        # Phase 1: Exploration
742        self.generate_sobol_samples(n_samples)
743        eval_results = self.evaluate_configurations(
744            X_train, X_test, y_train, y_test, model_class, model_kwargs
745        )
746
747        # Phase 2: Surrogate fitting
748        surrogate, metrics = self.fit_surrogate(surrogate, target_gap)
749
750        # Phase 3: Analysis
751        sobol_indices = self.compute_sobol_indices()
752
753        # Phase 4: Optimization
754        opt_results = self.optimize_acquisition()
755
756        return {
757            "evaluation_results": eval_results,
758            "surrogate_metrics": metrics,
759            "sobol_indices": sobol_indices,
760            "optimization": opt_results,
761        }