From 4b4e13ebd21142d6b6526015549d9f700a2a398b Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 27 Aug 2025 10:45:21 -0700 Subject: [PATCH 1/8] abs - removing constraint on variance parameter --- modules/performUQ/common/gp_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/performUQ/common/gp_model.py b/modules/performUQ/common/gp_model.py index 37804d4e1..4ce051e40 100644 --- a/modules/performUQ/common/gp_model.py +++ b/modules/performUQ/common/gp_model.py @@ -169,8 +169,8 @@ def _set_kernel_hyperparameter_bounds(self, kernel): # For standardized inputs, set reasonable bounds if hasattr(kernel, 'lengthscale'): kernel.lengthscale.constrain_bounded(0.01, 10.0) - if hasattr(kernel, 'variance'): - kernel.variance.constrain_bounded(0.01, 100.0) + # if hasattr(kernel, 'variance'): + # kernel.variance.constrain_bounded(0.01, 100.0) # TODO (ABS): For unscaled inputs, use default bounds or wider ranges def _fit_pca_and_create_models(self, *, reoptimize=True, num_random_restarts=10): From 10cf00ea43bd321885b3773f6059761f2c4df36b Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 27 Aug 2025 10:46:43 -0700 Subject: [PATCH 2/8] abs - using unique samples and updating exploitation criteria calculation --- modules/performUQ/common/adaptive_doe.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index 55f1b2d61..df2acb33d 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -141,6 +141,13 @@ def select_training_points( # 1. Fixed candidate pool if use_mse_w: + # get unique MCI samples to avoid redundant computations + mci_samples, weights = remove_duplicate_inputs( + mci_samples, + weights + if weights is not None + else np.ones((mci_samples.shape[0], 1)), + ) candidate_pool = mci_samples.copy() else: bounds = compute_lhs_bounds(x_train, mci_samples, padding=0) @@ -164,10 +171,8 @@ def select_training_points( # 5. Vectorized acquisition if use_mse_w: - # Just return the average variance (same for all candidates) - acquisition_values = np.mean(pred_var) * np.ones( - (scaled_candidates.shape[0], 1) - ) + # Just return the predictive variance + acquisition_values = pred_var else: # Covariance matrix: shape (n_mci, n_candidates) k_matrix = self.gp_model_for_doe.kern.K( From c328a1e2ab0e3cb2dbce8077ac113675060eb34c Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 27 Aug 2025 14:11:53 -0700 Subject: [PATCH 3/8] abs - increasing the upper bound on lengthscale --- modules/performUQ/common/gp_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/performUQ/common/gp_model.py b/modules/performUQ/common/gp_model.py index 4ce051e40..e08f0f461 100644 --- a/modules/performUQ/common/gp_model.py +++ b/modules/performUQ/common/gp_model.py @@ -168,7 +168,7 @@ def _set_kernel_hyperparameter_bounds(self, kernel): if self.scale_inputs: # For standardized inputs, set reasonable bounds if hasattr(kernel, 'lengthscale'): - kernel.lengthscale.constrain_bounded(0.01, 10.0) + kernel.lengthscale.constrain_bounded(0.01, 100.0) # if hasattr(kernel, 'variance'): # kernel.variance.constrain_bounded(0.01, 100.0) # TODO (ABS): For unscaled inputs, use default bounds or wider ranges From 3fe4c0092f04b245d702d1e91827f93a312cc855 Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 27 Aug 2025 14:12:45 -0700 Subject: [PATCH 4/8] abs - maintaining consistency in indices for mse criteria --- modules/performUQ/common/adaptive_doe.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index df2acb33d..39b161089 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -200,7 +200,12 @@ def select_training_points( # 7. Update training set and candidate pool x_train = np.vstack([x_train, next_point]) candidate_pool = np.delete(candidate_pool, idx, axis=0) - + if use_mse_w: + # For MSE: maintain alignment by deleting from mci_samples too + mci_samples = np.delete(mci_samples, idx, axis=0) + if weights is not None: + weights = np.delete(weights, idx, axis=0) + # For IMSEw: keep mci_samples constant (integration domain unchanged) return np.array(selected_points) def write_gp_for_doe_to_json(self, filepath: str | Path): From 3e6e82969701c6d167481572f6022ae86e3acb5d Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Tue, 2 Sep 2025 10:37:32 -0700 Subject: [PATCH 5/8] abs - fixing correlation calculation in imsew criterion --- modules/performUQ/common/adaptive_doe.py | 72 ++++++++++++++++-------- 1 file changed, 47 insertions(+), 25 deletions(-) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index 39b161089..5ac2abe7b 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -108,7 +108,7 @@ def select_training_points( seed=None, ): """ - Efficient sequential DoE using MSEw or IMSEw with vectorized acquisition. + Efficient sequential DoE using MSEw (equation 10) or IMSEw (equation 9) from Taflanidis et al. (https://doi.org/10.1016/j.ymssp.2025.113014). Parameters ---------- @@ -119,9 +119,9 @@ def select_training_points( mci_samples : array-like Monte Carlo integration samples (original space). use_mse_w : bool - Whether to use MSEw (True) or IMSEw (False). + Whether to use MSE equation (10) (True) or IMSEw equation (9) (False). weights : array-like or None - Optional importance weights. + Optional importance weights for integration points. n_samples : int Number of candidate points to generate (IMSEw only). seed : int or None @@ -139,9 +139,9 @@ def select_training_points( selected_points = [] - # 1. Fixed candidate pool + # 1. Setup candidate pool based on acquisition function if use_mse_w: - # get unique MCI samples to avoid redundant computations + # MSE (equation 10): candidates are the integration points themselves mci_samples, weights = remove_duplicate_inputs( mci_samples, weights @@ -150,62 +150,84 @@ def select_training_points( ) candidate_pool = mci_samples.copy() else: + # IMSEw (equation 9): generate separate candidate pool via LHS bounds = compute_lhs_bounds(x_train, mci_samples, padding=0) candidate_pool = generate_lhs_candidates(n_samples, bounds, seed=seed) + # 2. Sequential selection for _ in range(n_points): - # 2. Scale everything + # Scale all arrays scaled_x_train = self._scale(x_train) scaled_candidates = self._scale(candidate_pool) scaled_mci_samples = self._scale(mci_samples) - # 3. Set GP inputs + # Set GP with current training set self.gp_model_for_doe.set_XY( scaled_x_train, np.zeros((scaled_x_train.shape[0], 1)) ) - # 4. Predict variance at integration points - _, pred_var = self.gp_model_for_doe.predict(scaled_mci_samples) - if weights is not None: - pred_var *= weights.reshape(-1, 1) # shape (n_mci, 1) - - # 5. Vectorized acquisition + # Compute acquisition values based on selected method if use_mse_w: - # Just return the predictive variance + # MSE (equation 10): σ²(θ_new|Θ^(k)) + _, pred_var = self.gp_model_for_doe.predict(scaled_candidates) + + # Apply importance weights if provided + if weights is not None: + pred_var *= weights.reshape(-1, 1) + acquisition_values = pred_var + else: - # Covariance matrix: shape (n_mci, n_candidates) + # IMSEw (equation 9): (1/N_s) * Σ[R(θ^(l), θ_new)^β * σ²(θ^(l)|Θ^(k))] + + # Predict variance at integration points + _, pred_var = self.gp_model_for_doe.predict(scaled_mci_samples) + + # Apply importance weights to integration points + if weights is not None: + pred_var *= weights.reshape(-1, 1) + + # Compute covariance K(θ^(l), θ_new) k_matrix = self.gp_model_for_doe.kern.K( scaled_mci_samples, scaled_candidates ) - # Normalize to get correlation matrix - kernel_var = self.gp_model_for_doe.kern.variance[0] # scalar - corr_matrix = k_matrix / kernel_var # element-wise division + # Get diagonal kernel values for covariance normalization + k_diag_mci = np.diag( + self.gp_model_for_doe.kern.K(scaled_mci_samples) + ) + k_diag_candidates = np.diag( + self.gp_model_for_doe.kern.K(scaled_candidates) + ) + # Correlations : R(θ^(l), θ_new) + corr_matrix = k_matrix / np.sqrt( + np.outer(k_diag_mci, k_diag_candidates) + ) + + # Apply beta exponent beta = 2.0 * scaled_x_train.shape[1] - # beta = 1.0 * scaled_x_train.shape[1] - # Vectorized IMSEw: (corr^beta) * pred_var → mean over MCI samples - weighted_var = ( - corr_matrix**beta - ).T @ pred_var # shape (n_candidates, 1) + # IMSEw computation: (corr^beta).T @ pred_var / N_s + weighted_var = (corr_matrix**beta).T @ pred_var acquisition_values = weighted_var / scaled_mci_samples.shape[0] - # 6. Select best candidate + # Select candidate with maximum acquisition value idx = np.argmax(acquisition_values) next_point = candidate_pool[idx] selected_points.append(next_point) - # 7. Update training set and candidate pool + # Update training set and candidate pool x_train = np.vstack([x_train, next_point]) candidate_pool = np.delete(candidate_pool, idx, axis=0) + if use_mse_w: # For MSE: maintain alignment by deleting from mci_samples too mci_samples = np.delete(mci_samples, idx, axis=0) if weights is not None: weights = np.delete(weights, idx, axis=0) # For IMSEw: keep mci_samples constant (integration domain unchanged) + return np.array(selected_points) def write_gp_for_doe_to_json(self, filepath: str | Path): From e63343f6200ab465bc945aebd6d331a8778eff2d Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 3 Sep 2025 13:10:45 -0700 Subject: [PATCH 6/8] abs- adding output scaling in gp when not using pca, and setting reasonable bounds on lengthscales with or without input scaling --- modules/performUQ/common/gp_model.py | 318 ++++++++++++++++++++------- 1 file changed, 238 insertions(+), 80 deletions(-) diff --git a/modules/performUQ/common/gp_model.py b/modules/performUQ/common/gp_model.py index e08f0f461..cb4ea793a 100644 --- a/modules/performUQ/common/gp_model.py +++ b/modules/performUQ/common/gp_model.py @@ -71,6 +71,7 @@ class GaussianProcessModelSettings(BaseModel): use_pca: bool = False pca_threshold: float = 0.999 scale_inputs: bool = True + scale_outputs: bool = True @model_validator(mode='after') def check_nugget_value(self): # noqa: D102 @@ -105,13 +106,14 @@ def __init__( self.scale_inputs = settings.scale_inputs self.pca = ( - PrincipalComponentAnalysis(self.pca_threshold, perform_scaling=True) + PrincipalComponentAnalysis(self.pca_threshold, perform_scaling=False) if self.use_pca else None ) self.latent_dimension = None self.input_scaler = SafeStandardScaler() if self.scale_inputs else None + self.output_scaler = SafeStandardScaler() if settings.scale_outputs else None self.kernel = self._create_kernel() @@ -122,6 +124,56 @@ def __init__( self.y_train = None self.x_train_scaled = None + def apply_input_scaling(self, x, *, fit=False): + """Scale inputs using the input scaler.""" + if self.input_scaler is None: + return x + + if fit: + return self.input_scaler.fit_transform(x) + return self.input_scaler.transform(x) + + def _preprocess_outputs(self, y, *, fit: bool) -> np.ndarray: + """Scale outputs and apply PCA. If fit=False, use existing transformers.""" + y_processed = y + if self.output_scaler is not None: + y_processed = ( + self.output_scaler.fit_transform(y_processed) + if fit + else self.output_scaler.transform(y_processed) + ) + if self.pca is not None: + if fit: + self.pca.fit(y_processed) + y_processed = self.pca.project_to_latent_space(y_processed) + return y_processed + + def _postprocess_outputs(self, y_latent): + """Single point for all output postprocessing.""" + y_processed = y_latent + + # Step 1: Inverse PCA if used + if self.pca is not None: + y_processed = self.pca.project_back_to_original_space(y_processed) + + # Step 2: Inverse scaling if used + if self.output_scaler is not None: + y_processed = self.output_scaler.inverse_transform(y_processed) + + return y_processed + + # def _create_surrogate(self): + # use_gpy_normalizer = self.output_scaler is None + # model = GPy.models.GPRegression( + # X=np.zeros((1, self.input_dimension), dtype=float), + # Y=np.zeros((1, 1), dtype=float), + # kernel=self.kernel.copy(), + # mean_function=None, + # normalizer=use_gpy_normalizer, + # ) + # self._configure_nugget(model) + # return model + def _create_kernel(self): if self.kernel_type == 'exponential': return GPy.kern.Exponential(input_dim=self.input_dimension, ARD=self.ARD) @@ -134,26 +186,6 @@ def _create_kernel(self): msg = f"Unknown kernel_type '{self.kernel_type}'" raise ValueError(msg) - def apply_input_scaling(self, x, *, fit=False): - """Scale inputs using the input scaler.""" - if self.input_scaler is None: - return x - - if fit: - return self.input_scaler.fit_transform(x) - return self.input_scaler.transform(x) - - def _create_surrogate(self): - model = GPy.models.GPRegression( - X=np.zeros((1, self.input_dimension), dtype=float), - Y=np.zeros((1, 1), dtype=float), - kernel=self.kernel.copy(), - mean_function=None, - normalizer=True, - ) - self._configure_nugget(model) - return model - def _configure_nugget(self, model: GPy.models.GPRegression): """Configure the nugget (Gaussian noise variance) for each GP model.""" if self.settings.fix_nugget: @@ -163,19 +195,92 @@ def _configure_nugget(self, model: GPy.models.GPRegression): model.Gaussian_noise.variance = self.settings.nugget_value model.likelihood.variance.constrain_bounded(1e-8, 1e-3) + # def _set_kernel_hyperparameter_bounds(self, kernel): + # """Set reasonable hyperparameter bounds for scaled inputs.""" + # if self.scale_inputs: + # # For standardized inputs, set reasonable bounds + # if hasattr(kernel, 'lengthscale'): + # kernel.lengthscale.constrain_bounded(0.01, 10.0) + # # if hasattr(kernel, 'variance'): + # # kernel.variance.constrain_bounded(0.01, 100.0) + # # TODO (ABS): For unscaled inputs, use default bounds or wider ranges + def _set_kernel_hyperparameter_bounds(self, kernel): - """Set reasonable hyperparameter bounds for scaled inputs.""" - if self.scale_inputs: - # For standardized inputs, set reasonable bounds - if hasattr(kernel, 'lengthscale'): - kernel.lengthscale.constrain_bounded(0.01, 100.0) - # if hasattr(kernel, 'variance'): - # kernel.variance.constrain_bounded(0.01, 100.0) - # TODO (ABS): For unscaled inputs, use default bounds or wider ranges + """Set reasonable hyperparameter bounds for inputs.""" + if hasattr(kernel, 'lengthscale'): + if self.scale_inputs: + # Standard approach for scaled inputs + kernel.lengthscale.constrain_bounded(0.01, 10.0) + # Engineering approach: bounds based on physical scales + elif hasattr(self, 'x_train') and self.x_train is not None: + input_ranges = np.ptp(self.x_train, axis=0) + + # Handle zero-variance dimensions + input_ranges = np.where(input_ranges == 0, 1.0, input_ranges) + + # Apply bounds to lengthscale parameter + if self.ARD and len(input_ranges) > 1: + # For ARD kernels, set bounds per dimension using GPy's indexing + for nx in range(len(input_ranges)): + lb = input_ranges[nx] * 0.01 # 1% of range + ub = input_ranges[nx] * 10.0 # 10x range + + # Ensure valid bounds + lb = max(lb, 1e-6) + if lb >= ub: + lb = ub * 0.01 + + # Use GPy's double bracket indexing for individual lengthscale parameters + kernel.lengthscale[[nx]].constrain_bounded( + lb, ub, warning=False + ) + else: + # Single lengthscale case + lb = float(np.min(input_ranges) * 0.01) + ub = float(np.max(input_ranges) * 10.0) + lb = max(lb, 1e-6) + kernel.lengthscale.constrain_bounded(lb, ub) + + def _initialize_kernel_lengthscales(self, kernel, *, force_initialization=False): + """Initialize kernel lengthscales to reasonable values based on input data. + + Parameters + ---------- + kernel : GPy kernel + Kernel to initialize + force_initialization : bool + If True, always initialize. If False, only initialize if lengthscales + are at default values (close to 1.0). + """ + if not hasattr(kernel, 'lengthscale'): + return + + # Check if training data is available + if self.x_train is None: + return # Can't initialize without training data + + # Check if hyperparameters appear to be at default values + if not force_initialization: + current_lengthscales = kernel.lengthscale.values.flatten() # noqa: PD011 + if not np.allclose(current_lengthscales, 1.0, rtol=0.01): + return # Skip initialization - hyperparameters appear optimized + + # Simple consistent initialization: use input data ranges + input_ranges = np.ptp(self.x_train, axis=0) + # Handle zero-variance dimensions + input_ranges = np.where(input_ranges == 0, 1.0, input_ranges) + + # Use 20% of range as initial lengthscale + initial_lengthscales = input_ranges * 0.2 + + # Apply to kernel + if self.ARD and len(initial_lengthscales) > 1: + kernel.lengthscale[:] = initial_lengthscales + else: + kernel.lengthscale = float(np.median(initial_lengthscales)) def _fit_pca_and_create_models(self, *, reoptimize=True, num_random_restarts=10): self.x_train_scaled = self.apply_input_scaling(self.x_train, fit=True) - if self.scale_inputs: self.loginfo( f'Input scaling enabled. Scaled training inputs with shape {self.x_train_scaled.shape}' @@ -183,9 +288,15 @@ def _fit_pca_and_create_models(self, *, reoptimize=True, num_random_restarts=10) else: self.loginfo('Input scaling disabled. Using original input scale.') + y_latent = self._preprocess_outputs(self.y_train, fit=True) + if self.output_scaler is not None: + self.loginfo('Output scaling enabled. Scaled training outputs.') + else: + self.loginfo('Output scaling disabled. Using original output scale.') + if self.pca is not None: - self.pca.fit(self.y_train) - y_latent = self.pca.project_to_latent_space(self.y_train) + # self.pca.fit(self.y_train) + # y_latent = self.pca.project_to_latent_space(self.y_train) self.latent_dimension = int(self.pca.n_components) # type: ignore self.output_dimension = self.latent_dimension self.loginfo( @@ -194,7 +305,7 @@ def _fit_pca_and_create_models(self, *, reoptimize=True, num_random_restarts=10) f'{self.pca_threshold:.2%} variance.' ) else: - y_latent = self.y_train + # y_latent = self.y_train self.output_dimension = self.model_output_dimension self.loginfo( f'Using model outputs of dimension {self.model_output_dimension} ' @@ -217,14 +328,20 @@ def _fit_pca_and_create_models(self, *, reoptimize=True, num_random_restarts=10) y_detrended = y kernel_copy = self.kernel.copy() + # Data-adaptive lengthscale initialization (always for new models) + self._initialize_kernel_lengthscales( + kernel_copy, force_initialization=True + ) self._set_kernel_hyperparameter_bounds(kernel_copy) + use_gpy_normalizer = self.output_scaler is None + gp = GPy.models.GPRegression( X=x, Y=y_detrended, kernel=kernel_copy, mean_function=None, - normalizer=True, + normalizer=use_gpy_normalizer, ) self._configure_nugget(gp) if reoptimize: @@ -295,13 +412,14 @@ def update_training_dataset( self.y_train = outputs self.x_train_scaled = self.apply_input_scaling(self.x_train, fit=True) + y_latent = self._preprocess_outputs(self.y_train, fit=True) if self.pca is not None: - self.pca.fit(self.y_train) - y_latent = self.pca.project_to_latent_space(self.y_train) + # self.pca.fit(self.y_train) + # y_latent = self.pca.project_to_latent_space(self.y_train) new_latent_dim = int(self.pca.n_components) # type: ignore else: - y_latent = self.y_train + # y_latent = self.y_train new_latent_dim = self.model_output_dimension if ( @@ -321,6 +439,11 @@ def update_training_dataset( ) else: y_detrended = y + + # Initialize only if hyperparameters appear unoptimized + self._initialize_kernel_lengthscales( + self.model[i].kern, force_initialization=False + ) self.model[i].set_XY(self.x_train_scaled, y_detrended) if reoptimize: self.model[i].optimize_restarts(num_random_restarts) @@ -361,11 +484,13 @@ def predict(self, x_predict): y_mean[:, i] = mean.flatten() y_var[:, i] = var.flatten() - if self.pca is not None: - y_mean = self.pca.project_back_to_original_space(y_mean) - # y_var = self.pca.inverse_transform_variance(y_var) # Not transforming variance for now + # if self.pca is not None: + # y_mean = self.pca.project_back_to_original_space(y_mean) + # # y_var = self.pca.inverse_transform_variance(y_var) # Not transforming variance for now + + y_mean_final = self._postprocess_outputs(y_mean) - return y_mean, y_var # variance is in latent space, if pca is used + return y_mean_final, y_var # variance is in latent space, if pca is used def loo_predictions(self, x_train, y_train, epsilon=1e-8): """ @@ -394,15 +519,16 @@ def loo_predictions(self, x_train, y_train, epsilon=1e-8): raise ValueError(msg) x_train_scaled = self.apply_input_scaling(x_train, fit=False) + y_train_processed = self._preprocess_outputs(y_train, fit=False) - if self.pca is not None: - y_train = self.pca.project_to_latent_space(y_train) + # if self.pca is not None: + # y_train = self.pca.project_to_latent_space(y_train) - loo_pred_latent = np.zeros_like(y_train, dtype=np.float64) + loo_pred_latent = np.zeros_like(y_train_processed, dtype=np.float64) for i in range(self.output_dimension): x = x_train_scaled - y = y_train[:, i].reshape(-1, 1) + y = y_train_processed[:, i].reshape(-1, 1) # Step 1: Detrend if self.linear_models[i] is not None: @@ -429,12 +555,13 @@ def loo_predictions(self, x_train, y_train, epsilon=1e-8): # Step 4: Add back trend to get final prediction loo_pred_latent[:, i] = loo_gp + trend.flatten() - # Step 5: Inverse PCA if used - if self.pca is not None: - loo_pred_latent = self.pca.project_back_to_original_space( - loo_pred_latent - ) - return loo_pred_latent + # Step 5: Postprocess outputs + # if self.pca is not None: + # loo_pred_latent = self.pca.project_back_to_original_space( + # loo_pred_latent + # ) + loo_pred_latent = self._postprocess_outputs(loo_pred_latent) + return loo_pred_latent # noqa: RET504 @property def pca_info(self): @@ -489,6 +616,7 @@ def write_model_parameters_to_json( 'use_pca': self.use_pca, 'pca_threshold': self.pca_threshold, 'scale_inputs': self.scale_inputs, + 'scale_outputs': bool(self.output_scaler is not None), 'pca_info': self.pca_info, 'models': [], } @@ -500,6 +628,13 @@ def write_model_parameters_to_json( 'n_features_in_': int(self.input_scaler.n_features_in_), # type: ignore } + if self.output_scaler is not None: + model_params['output_scaler'] = { + 'mean_': self.output_scaler.mean_.tolist(), # type: ignore + 'scale_': self.output_scaler.scale_.tolist(), # type: ignore + 'n_features_in_': int(self.output_scaler.n_features_in_), + } + if include_training_data: if self.x_train is None or self.y_train is None: self.logger.warning( @@ -547,19 +682,11 @@ def load_model_parameters_from_json( # noqa: C901 *, has_training_data: bool = True, ): - """ - Load GP model settings, kernel, linear regression, and optionally training data from a JSON file. - - Parameters - ---------- - filepath : str or Path - Path to the JSON file. - has_training_data : bool - Whether the JSON file includes training data (default is True). - """ + """Load GP model settings, kernel, linear regression, and optionally training data from a JSON file.""" with Path(filepath).open() as f: model_params = json.load(f) + # --- Settings / flags --- self.input_dimension = model_params['input_dimension'] self.output_dimension = model_params['output_dimension'] self.ARD = model_params['ARD'] @@ -570,60 +697,82 @@ def load_model_parameters_from_json( # noqa: C901 self.use_pca = model_params['use_pca'] self.pca_threshold = model_params['pca_threshold'] self.scale_inputs = model_params.get('scale_inputs', True) + scale_outputs_flag = model_params.get('scale_outputs', True) + # --- Input scaler --- + self.input_scaler = None if self.scale_inputs and 'input_scaler' in model_params: self.input_scaler = SafeStandardScaler() - scaler_params = model_params['input_scaler'] - self.input_scaler.mean_ = np.array(scaler_params['mean_']) - self.input_scaler.scale_ = np.array(scaler_params['scale_']) - self.input_scaler.n_features_in_ = scaler_params['n_features_in_'] # type: ignore + s = model_params['input_scaler'] + self.input_scaler.mean_ = np.array(s['mean_']) + self.input_scaler.scale_ = np.array(s['scale_']) + self.input_scaler.n_features_in_ = s['n_features_in_'] + + # --- Output scaler --- + self.output_scaler = None + if scale_outputs_flag: + self.output_scaler = SafeStandardScaler() + if 'output_scaler' in model_params: + s = model_params['output_scaler'] + self.output_scaler.mean_ = np.array(s['mean_']) + self.output_scaler.scale_ = np.array(s['scale_']) + self.output_scaler.n_features_in_ = s['n_features_in_'] + # If missing, it's fine: we'll fit during preprocessing when fit=True. + + # Keep as local (don't assign to property) + _pca_info_from_file = model_params.get('pca_info') + + # Reset PCA/latent bookkeeping; let training path fit as needed + self.pca = None + self.latent_dimension = None + if self.use_pca: + # Keep perform_scaling=False for consistency with external scalers + self.pca = PrincipalComponentAnalysis( + self.pca_threshold, perform_scaling=False + ) if has_training_data: + # --- Training data path: let _fit_pca_and_create_models do all fitting --- self.x_train = np.array(model_params['x_train']) self.y_train = np.array(model_params['y_train']) + # Transform inputs using restored scaler (no fit) self.x_train_scaled = self.apply_input_scaling(self.x_train, fit=False) - self.pca = None - self.latent_dimension = None - if self.use_pca: - self.pca_info = model_params['pca_info'] # type: ignore - self.pca = PrincipalComponentAnalysis( - self.pca_threshold, perform_scaling=True - ) - self.pca.fit(self.y_train) - + # Build models without re-optimizing (since we'll overwrite params below) self._fit_pca_and_create_models(reoptimize=False) else: + # --- No training data: rebuild empty models and then overwrite params --- models_data = model_params['models'] + # Trust the JSON if it differs from the header output_dimension if len(models_data) != self.output_dimension: - msg = f'Expected {self.output_dimension} models, but got {len(models_data)}' - raise ValueError(msg) + self.output_dimension = len(models_data) self.model = [] self.linear_models = [] base_kernel = self._create_kernel() - for _ in models_data: x_dummy = np.zeros((1, self.input_dimension)) y_dummy = np.zeros((1, 1)) kernel_copy = base_kernel.copy() self._set_kernel_hyperparameter_bounds(kernel_copy) + use_gpy_normalizer = self.output_scaler is None gp = GPy.models.GPRegression( X=x_dummy, Y=y_dummy, kernel=kernel_copy, mean_function=None, - normalizer=True, + normalizer=use_gpy_normalizer, ) self._configure_nugget(gp) self.model.append(gp) self.linear_models.append(None) - # Overwrite kernel and linear model parameters - for i, m in enumerate(model_params['models']): + # --- Overwrite kernel / nugget / linear trend params from file --- + models_data = model_params['models'] + for i, m in enumerate(models_data): gp = self.model[i] for name, value in m['kernel_parameters'].items(): @@ -637,9 +786,12 @@ def load_model_parameters_from_json( # noqa: C901 lr_params = m.get('linear_regression') if lr_params is not None: lr = LinearRegression() + # shape: (1, n_features) to match sklearn internals lr.coef_ = np.array(lr_params['coefficients']).reshape(1, -1) lr.intercept_ = np.array([lr_params['intercept']]) self.linear_models[i] = lr + else: + self.linear_models[i] = None # =============== @@ -838,6 +990,12 @@ def parse_gp_arguments(args=None) -> dict: default=True, help='Whether to scale inputs before training the GP.', ) + optional.add_argument( + '--scale_outputs', + type=lambda x: x.lower() in ('true', '1'), + default=True, + help='Whether to scale outputs before training the GP.', + ) parsed_args = parser.parse_args(args) return vars(parsed_args) From 8b6bbbfa448cafafbba0cbda8092567229aa3aff Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 3 Sep 2025 13:12:47 -0700 Subject: [PATCH 7/8] abs - weighing hyperparameters by variance when not using pca --- modules/performUQ/common/adaptive_doe.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index 5ac2abe7b..3097b5935 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -57,8 +57,18 @@ def _hyperparameters_for_doe(self): else: w = explained_variance_ratio / np.sum(explained_variance_ratio) else: - n_models = len(self.gp_model.model) - w = np.full(n_models, 1.0 / n_models) + if self.gp_model.output_scaler is not None: + Y_proc = self.gp_model.output_scaler.transform(self.gp_model.y_train) + else: + Y_proc = self.gp_model.y_train + + var_per_output = np.nanvar(Y_proc, axis=0, ddof=1) # shape (p,) + var_per_output = np.where( + np.isfinite(var_per_output) & (var_per_output > 0), + var_per_output, + 1.0, + ) + w = var_per_output / np.sum(var_per_output) hyperparameters_matrix = [ np.atleast_2d(model.kern.param_array) for model in self.gp_model.model From 77f23b1d138cd8c279eddbadeea49f895a327728 Mon Sep 17 00:00:00 2001 From: bsaakash <11618528+bsaakash@users.noreply.github.com> Date: Wed, 3 Sep 2025 14:04:12 -0700 Subject: [PATCH 8/8] abs - refactoring gp hyperparameter aggregation --- modules/performUQ/common/adaptive_doe.py | 156 ++++++++++++++--------- 1 file changed, 99 insertions(+), 57 deletions(-) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index 3097b5935..4c5c19a85 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -39,36 +39,95 @@ def __init__(self, gp_model): def _scale(self, x): return self.gp_model.apply_input_scaling(x, fit=False) - def _hyperparameters_for_doe(self): + def _safe_normalize(self, v: np.ndarray) -> np.ndarray: + s = float(np.nansum(v)) + if not np.isfinite(s) or s <= 0: + return np.full_like(v, 1.0 / max(len(v), 1)) + return v / s + + def _compute_doe_weights(self) -> tuple[np.ndarray, dict]: """ - Compute the weighted average of kernel hyperparameters for DoE. + Compute weights for aggregating across GP output components. - If PCA is used, weights are based on explained variance. - Otherwise, uniform averaging is used. + Return (weights, meta) for aggregating across GP output components: + - PCA on -> eigenvalue weights (explained_variance_ratio) + - PCA off -> per-output variance in processed training space. """ + n_models = len(self.gp_model.model) + if self.gp_model.use_pca: - pca_info = self.gp_model.pca_info - n_components = pca_info['n_components'] - explained_variance_ratio = np.asarray( - pca_info['explained_variance_ratio'] - )[:n_components] - if np.sum(explained_variance_ratio) == 0: - w = np.full(n_components, 1.0 / n_components) + pinfo = self.gp_model.pca_info + ncomp = int(pinfo['n_components']) + evr = np.asarray(pinfo['explained_variance_ratio'])[:ncomp].astype(float) + w = self._safe_normalize(evr) + if len(w) != n_models: # defensive + w = np.full(n_models, 1.0 / max(n_models, 1)) + basis = 'uniform_fallback' + basis_vals = np.ones(n_models, dtype=float) else: - w = explained_variance_ratio / np.sum(explained_variance_ratio) + basis = 'pca_explained_variance_ratio' + basis_vals = evr + return w, { + 'mode': 'pca', + 'basis': basis, + 'basis_values': basis_vals.tolist(), + 'sum': float(np.sum(w)), + 'count': int(len(w)), + } + + # no PCA: weight by per-output variance in the space the GP trained on + Y = getattr(self.gp_model, 'y_train', None) # noqa: N806 + if Y is None or Y.size == 0: + w = np.full(n_models, 1.0 / max(n_models, 1)) + return w, { + 'mode': 'no_pca', + 'basis': 'uniform_no_training_data', + 'basis_values': [1.0] * n_models, + 'sum': float(np.sum(w)), + 'count': int(len(w)), + } + + if self.gp_model.output_scaler is not None: + Y_proc = self.gp_model.output_scaler.transform(Y) # noqa: N806 else: - if self.gp_model.output_scaler is not None: - Y_proc = self.gp_model.output_scaler.transform(self.gp_model.y_train) + Y_proc = Y # noqa: N806 + + raw = np.nanvar(Y_proc, axis=0, ddof=1).astype(float) + raw = np.where(np.isfinite(raw) & (raw > 0), raw, 1.0) + w = self._safe_normalize(raw) + + if len(w) != n_models: # defensive + if len(w) > n_models: + basis_vals = raw[:n_models] + w = self._safe_normalize(basis_vals) + basis = 'output_variance_truncated' else: - Y_proc = self.gp_model.y_train + pad = np.full( + n_models - len(w), float(np.mean(raw)) if raw.size else 1.0 + ) + basis_vals = np.concatenate([raw, pad]) + w = self._safe_normalize(basis_vals) + basis = 'output_variance_padded' + else: + basis = 'output_variance' + basis_vals = raw + + return w, { + 'mode': 'no_pca', + 'basis': basis, + 'basis_values': basis_vals.tolist(), + 'sum': float(np.sum(w)), + 'count': int(len(w)), + } - var_per_output = np.nanvar(Y_proc, axis=0, ddof=1) # shape (p,) - var_per_output = np.where( - np.isfinite(var_per_output) & (var_per_output > 0), - var_per_output, - 1.0, - ) - w = var_per_output / np.sum(var_per_output) + def _hyperparameters_for_doe(self): + """ + Compute the weighted average of kernel hyperparameters for DoE. + + If PCA is used, weights are based on explained variance. + Otherwise, uniform averaging is used. + """ + w, _ = self._compute_doe_weights() hyperparameters_matrix = [ np.atleast_2d(model.kern.param_array) for model in self.gp_model.model @@ -241,61 +300,44 @@ def select_training_points( return np.array(selected_points) def write_gp_for_doe_to_json(self, filepath: str | Path): - """ - Write DoE GP kernel hyperparameters and contributing model param_arrays to JSON. - - Parameters - ---------- - filepath : str or Path - Output file path. - """ + """Write DoE GP kernel hyperparameters and contributing model param_arrays to JSON.""" if not hasattr(self, 'gp_model_for_doe'): msg = 'gp_model_for_doe has not been initialized.' raise RuntimeError(msg) - kernel = self.gp_model_for_doe.kern + # Ensure we have the current weighted vector available + if not hasattr(self, 'weighted_hyperparameters'): + _ = self._hyperparameters_for_doe() + + # Reuse the exact same weights/provenance used for aggregation + weights, weights_meta = self._compute_doe_weights() - # Detailed hyperparameters for the DoE model + kernel = self.gp_model_for_doe.kern doe_hyperparams = { - param.name: { - 'value': param.values.tolist() # noqa: PD011 - if param.size > 1 - else float(param.values), - 'shape': param.shape, + p.name: { + 'value': p.values.tolist() if p.size > 1 else float(p.values), # noqa: PD011 + 'shape': p.shape, } - for param in kernel.parameters + for p in kernel.parameters } - # Aggregation weights - if self.gp_model.use_pca: - weights = np.asarray(self.gp_model.pca_info['explained_variance_ratio'])[ - : self.gp_model.pca_info['n_components'] - ] - weights = (weights / np.sum(weights)).tolist() - else: - weights = ( - np.ones(len(self.gp_model.model)) / len(self.gp_model.model) - ).tolist() - - # Contributing models' param_arrays only contributing_param_arrays = [ gp.kern.param_array.tolist() for gp in self.gp_model.model ] - # Output structure output = { 'doe_kernel_type': kernel.name, - 'doe_ARD': kernel.ARD, + 'doe_ARD': getattr(kernel, 'ARD', None), 'doe_hyperparameters': doe_hyperparams, 'weighted_param_array': self.weighted_hyperparameters.tolist(), - 'aggregation_weights': weights, + 'aggregation_weights': weights.tolist(), + 'aggregation_weights_meta': weights_meta, # provenance for reproducibility 'contributing_param_arrays': contributing_param_arrays, } - # Write JSON file - filepath = Path(filepath) - filepath.parent.mkdir(parents=True, exist_ok=True) - with filepath.open('w') as f: + path = Path(filepath) + path.parent.mkdir(parents=True, exist_ok=True) + with path.open('w') as f: json.dump(uq_utilities.make_json_serializable(output), f, indent=4)