|
| 1 | +# Distributed under the MIT License. |
| 2 | +# See LICENSE.txt for details. |
| 3 | + |
| 4 | +""" |
| 5 | +Gaussian Process Regression machine learning diagnostic functions library. |
| 6 | +Contains all the functions necessary to validate and plot the GPR model used |
| 7 | +to predict better low-eccentricity orbital parameter initial guesses. |
| 8 | +""" |
| 9 | + |
| 10 | +import logging |
| 11 | +import os |
| 12 | + |
| 13 | +import matplotlib.pyplot as plt |
| 14 | +import numpy as np |
| 15 | +import torch |
| 16 | + |
| 17 | +os.environ["OMP_NUM_THREADS"] = "1" |
| 18 | +from multiprocessing import Pool |
| 19 | + |
| 20 | +from SimulationSupport.gpr import predict_with_gpr_model, train_gpr_model |
| 21 | + |
| 22 | +logger = logging.getLogger(__name__) |
| 23 | + |
| 24 | + |
| 25 | +# Leave-one-out parallelization set up |
| 26 | +def _loo_single(i, X, Y): |
| 27 | + """ |
| 28 | + Run a single LOO iteration for index i. |
| 29 | +
|
| 30 | + Trains a GPR on all points except index i, then predicts the held-out point. |
| 31 | + Called in parallel by loo_crossval. |
| 32 | +
|
| 33 | + Args: |
| 34 | + i (int): Index of the held-out point |
| 35 | + X (np.ndarray): Input features, with shape (N, D) |
| 36 | + Y (np.ndarray): Target variable, with shape (N, ) |
| 37 | +
|
| 38 | + Returns: |
| 39 | + i (int): Index of the held-out point |
| 40 | + pred_mean (float): Predicted mean for the held-out point |
| 41 | + pred_std (float): Predicted std dev for the held-out point |
| 42 | + """ |
| 43 | + N = len(Y) |
| 44 | + |
| 45 | + # Create train and test split |
| 46 | + # Boolean mask: all True except index i (held out point) |
| 47 | + train_mask = np.ones(N, dtype=bool) |
| 48 | + train_mask[i] = False |
| 49 | + |
| 50 | + X_train = X[train_mask] |
| 51 | + Y_train = Y[train_mask] |
| 52 | + # Slice preserves the 2D shape expected by the model |
| 53 | + X_test = X[i : i + 1] |
| 54 | + |
| 55 | + # Train and predict |
| 56 | + model_loo, likelihood_loo = train_gpr_model(X_train, Y_train) |
| 57 | + pred_mean, pred_std = predict_with_gpr_model( |
| 58 | + X_test, model_loo, likelihood_loo |
| 59 | + ) |
| 60 | + |
| 61 | + return i, pred_mean[0], pred_std[0] |
| 62 | + |
| 63 | + |
| 64 | +# Leave-one-out cross-validation |
| 65 | +def loo_crossval( |
| 66 | + X: np.ndarray, |
| 67 | + Y: np.ndarray, |
| 68 | + target_name="Target", |
| 69 | + n_jobs=-1, |
| 70 | +): |
| 71 | + """ |
| 72 | + Perform Leave-One-Out Cross-Validation for a GPR model. |
| 73 | +
|
| 74 | + Trains N models (each omitting one point), predicts the held-out point, |
| 75 | + collect predictions and uncertainties, and then computes and plots summary metrics. |
| 76 | + This gives an unbiased estimate of generalization performance. |
| 77 | +
|
| 78 | + Args: |
| 79 | + X (np.ndarray): Input features, with shape (N, D) |
| 80 | + Y (np.ndarray): Target variable, with shape (N, ) |
| 81 | + target_name (str): Label for plots and print outputs |
| 82 | + n_jobs (int): Number of parallel workers (-1 uses all available cores) |
| 83 | +
|
| 84 | + Returns: |
| 85 | + predictions_loo (np.ndarray): LOO predicted values, with shape (N, ) |
| 86 | + uncertainties_loo (np.ndarray): LOO predicted std devs, with shape (N, ) |
| 87 | + """ |
| 88 | + N = len(Y) |
| 89 | + predictions_loo = np.zeros_like(Y) |
| 90 | + uncertainties_loo = np.zeros_like(Y) |
| 91 | + |
| 92 | + # Force single worker when GPU is used as parallel processes can't share a GPU |
| 93 | + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") |
| 94 | + if device.type == "cuda": |
| 95 | + n_jobs = 1 |
| 96 | + |
| 97 | + # multiprocessing.Pool uses None to use all available cores |
| 98 | + pool_size = None if n_jobs == -1 else n_jobs |
| 99 | + |
| 100 | + logger.info( |
| 101 | + f"Processing {N} LOO iterations for {target_name} using" |
| 102 | + f" {pool_size or os.cpu_count()} parallel workers..." |
| 103 | + ) |
| 104 | + |
| 105 | + # Build argument tuples for each LOO fold |
| 106 | + # because pool.starmap needs one tuple of positional args per call |
| 107 | + args = [(i, X, Y) for i in range(N)] |
| 108 | + |
| 109 | + if n_jobs == 1: |
| 110 | + # Run sequentially when there is a single worker |
| 111 | + results = [_loo_single(*a) for a in args] |
| 112 | + else: |
| 113 | + # Run LOO folds in parallel with multiple workers |
| 114 | + with Pool(processes=pool_size) as pool: |
| 115 | + results = pool.starmap(_loo_single, args) |
| 116 | + |
| 117 | + # Unpack results back into the predictions and uncertainties arrays |
| 118 | + for i, pred_mean, pred_std in results: |
| 119 | + predictions_loo[i] = pred_mean |
| 120 | + uncertainties_loo[i] = pred_std |
| 121 | + |
| 122 | + return predictions_loo, uncertainties_loo |
| 123 | + |
| 124 | + |
| 125 | +def plot_loo_crossval(Y, predictions_loo, target_name="Target", plot=False): |
| 126 | + """ |
| 127 | + Compute summary statistics for Leave-One-Out Cross-Validation results, and optionally |
| 128 | + plot the correlation. |
| 129 | +
|
| 130 | + Args: |
| 131 | + Y (np.ndarray): True target values, with shape (N, ) |
| 132 | + predictions_loo (np.ndarray): LOO predicted values, with shape (N, ) |
| 133 | + target_name (str): Label for plots and log outputs |
| 134 | + plot (bool): Whether to produce a correlation plot. Default is False. |
| 135 | +
|
| 136 | + Returns: |
| 137 | + rmse_loo (float): Root mean squared error of the LOO predictions |
| 138 | + mae_loo (float): Mean absolute error of the LOO predictions |
| 139 | + r_squared_loo (float): R^2 computed from the Pearson correlation |
| 140 | + """ |
| 141 | + |
| 142 | + Y_loo = Y # Same as the original Y for the multi input case |
| 143 | + |
| 144 | + # Calculate metrics - always computed, regardless of whether plot is requested |
| 145 | + # R^2 is computed from the Pearson correlation coefficient - |
| 146 | + # equivalent to the coefficient of determination for a linear fit through the origin |
| 147 | + correlation = np.corrcoef(Y_loo, predictions_loo)[0, 1] |
| 148 | + r_squared_loo = correlation**2 |
| 149 | + |
| 150 | + # Metrics with goal values |
| 151 | + rmse_loo = np.sqrt(np.mean((Y_loo - predictions_loo) ** 2)) |
| 152 | + mae_loo = np.mean(np.abs(Y_loo - predictions_loo)) |
| 153 | + y_range = Y_loo.max() - Y_loo.min() |
| 154 | + |
| 155 | + # Plot correlation |
| 156 | + if plot: |
| 157 | + plt.figure(figsize=(8, 6)) |
| 158 | + plt.scatter(Y_loo, predictions_loo, alpha=0.6, s=20) |
| 159 | + |
| 160 | + # y = x reference line: perfect predictions would lie exactly on this line |
| 161 | + min_val = min(Y_loo.min(), predictions_loo.min()) |
| 162 | + max_val = max(Y_loo.max(), predictions_loo.max()) |
| 163 | + plt.plot( |
| 164 | + [min_val, max_val], |
| 165 | + [min_val, max_val], |
| 166 | + "r--", |
| 167 | + lw=2, |
| 168 | + label="Perfect Correlation", |
| 169 | + ) |
| 170 | + |
| 171 | + # Labels and formatting |
| 172 | + plt.xlabel(f"True Δ{target_name}", fontsize=12) |
| 173 | + plt.ylabel(f"LOO Predicted Δ{target_name}", fontsize=12) |
| 174 | + plt.title(f"LOO: GPR Predictions vs True ({target_name})", fontsize=14) |
| 175 | + plt.grid(True, alpha=0.3) |
| 176 | + plt.legend() |
| 177 | + |
| 178 | + # Display R^2 |
| 179 | + plt.text( |
| 180 | + 0.95, |
| 181 | + 0.95, |
| 182 | + f"R² = {r_squared_loo:.4f}", |
| 183 | + transform=plt.gca().transAxes, |
| 184 | + fontsize=12, |
| 185 | + bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), |
| 186 | + horizontalalignment="right", |
| 187 | + ) |
| 188 | + plt.tight_layout() |
| 189 | + plt.show() |
| 190 | + |
| 191 | + # Log metrics with goal values for quick analysis |
| 192 | + logger.info(f"Leave-one-out Cross Validation Results ({target_name})") |
| 193 | + logger.info( |
| 194 | + f"RMSE: {rmse_loo:.6f} ({100 * rmse_loo / y_range:.6f}% of target" |
| 195 | + " range; goal: < 1 % of target range, lower is better)" |
| 196 | + ) |
| 197 | + logger.info( |
| 198 | + f"MAE: {mae_loo:.6f} ({100 * mae_loo / y_range:.6f}% of target range;" |
| 199 | + " goal: < 1 % of target range, lower is better)" |
| 200 | + ) |
| 201 | + logger.info( |
| 202 | + f"R²: {r_squared_loo:.4f} (goal: > 0.95 excellent, > 0.90 good, < 0.70" |
| 203 | + " poor)" |
| 204 | + ) |
| 205 | + logger.info(f"\n Dataset size: {len(Y_loo)} points") |
| 206 | + logger.info( |
| 207 | + f"Each model is trained on {len(Y_loo)-1} points, and tested on 1 point" |
| 208 | + ) |
| 209 | + logger.info("This provides an unbiased generalization estimate.") |
| 210 | + |
| 211 | + return rmse_loo, mae_loo, r_squared_loo |
| 212 | + |
| 213 | + |
| 214 | +# Leave-one-out residual computation and plotting |
| 215 | +def plot_loo_residuals(Y_loo, predictions_loo, target_name="Target", show=True): |
| 216 | + """ |
| 217 | + Calculate LOO prediction residuals, plot a histogram, and print summary statistics. |
| 218 | +
|
| 219 | + Args: |
| 220 | + Y_loo (np.ndarray): True target values from LOO cross validation |
| 221 | + predictions_loo (np.ndarray): Predicted values from LOO cross validation |
| 222 | + target_name (str): Name of target variable |
| 223 | + show (bool): Whether to produce a residual plot. Default is True. |
| 224 | +
|
| 225 | + Returns: |
| 226 | + residuals_loo (np.ndarray): Per-point residuals (true - predicted), with shape (N,) |
| 227 | +
|
| 228 | + """ |
| 229 | + |
| 230 | + # Compute residuals: LOO prediction error |
| 231 | + # Residual = true - predicted |
| 232 | + residuals_loo = Y_loo - predictions_loo |
| 233 | + |
| 234 | + # Make histogram |
| 235 | + plt.figure(figsize=(8, 5)) |
| 236 | + plt.hist(residuals_loo, bins=20, color="skyblue", edgecolor="k", alpha=0.8) |
| 237 | + # Vertical line at zero: residuals centered here indicate no systematic bias; |
| 238 | + # ideally the histogram is centered on this line; a shifted distribution suggests |
| 239 | + # the model is over- or under-predicting |
| 240 | + plt.axvline(0, color="r", linestyle="--", label="Zero Error") |
| 241 | + |
| 242 | + plt.title(f"LOO Residuals Histogram for {target_name}") |
| 243 | + plt.xlabel(" Residuals", fontsize=16) |
| 244 | + plt.ylabel("Count", fontsize=16) |
| 245 | + plt.tick_params(axis="both", which="major", labelsize=14) |
| 246 | + plt.grid(True, alpha=0.3) |
| 247 | + plt.legend(fontsize=14) |
| 248 | + plt.tight_layout() |
| 249 | + if show: |
| 250 | + plt.show() |
| 251 | + |
| 252 | + # Print statistics |
| 253 | + logger.info(f"Residual statistics for {target_name}:") |
| 254 | + logger.info(f"Mean residual: {np.mean(residuals_loo):.4e}") |
| 255 | + logger.info(f"Std of residuals: {np.std (residuals_loo):.4e}") |
| 256 | + logger.info(f"Max residual: {np.max (residuals_loo):.4e}") |
| 257 | + logger.info(f"Min residual: {np.min (residuals_loo):.4e}") |
| 258 | + |
| 259 | + return residuals_loo |
0 commit comments