abs - skipping duplicate training points during gp calibration

Author: bsaakash

modules/performUQ/common/adaptive_doe.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 import uq_utilities
+from gp_model import remove_duplicate_inputs
 from scipy.stats import qmc
 
 
@@ -131,6 +132,11 @@ def select_training_points(
         selected_points : np.ndarray of shape (n_points, d)
             Selected new training points.
         """
+        x_train = np.atleast_2d(x_train)
+        x_train, _ = remove_duplicate_inputs(
+            x_train, np.zeros((x_train.shape[0], 1))
+        )
+
         selected_points = []
 
         # 1. Fixed candidate pool
@@ -173,6 +179,7 @@ def select_training_points(
                 corr_matrix = k_matrix / kernel_var  # element-wise division
 
                 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 = (

modules/performUQ/common/gp_ab_algorithm.py

@@ -505,12 +505,16 @@ def _calculate_gcv(self, weights=None):
         """
         # start_time = time.time()
         self.loginfo('Calculating leave-one-out cross validation error measure gCV.')
+        x_train_unique, y_train_unique = (
+            self.current_gp_model.x_train,  # type: ignore
+            self.current_gp_model.y_train,  # type: ignore
+        )
         loo_predictions = self.current_gp_model.loo_predictions(  # type: ignore
-            self.inputs, self.outputs
+            x_train_unique, y_train_unique
         )
         loocv_measure = convergence_metrics.calculate_gcv(
             loo_predictions,
-            self.outputs,
+            y_train_unique,
             self.output_length_list,
             weights=weights,
             weight_combination=(2 / 3, 1 / 3),
@@ -521,21 +525,20 @@ def _calculate_gcv(self, weights=None):
 
     def _get_exploitation_candidates(self):
         """
-        Assemble candidate training points for exploitation by sampling from TMCMC stages
-        after the warm start, proportionally to the provided stage weights.
+        Assemble candidate training points for exploitation.
+
+        This is done by sampling from TMCMC stages after the warm start,
+        proportionally to the provided stage weights.
 
         Returns
         -------
         np.ndarray:
             The exploitation candidate training points of shape (n_candidates, n_parameters).
         dict:
             Dictionary mapping stage number to the number of points sampled from that stage.
-        """  # noqa: D205
-        # self.loginfo('Assembling candidate training points for exploitation.')
-        stage_weights = np.asarray(
-            self.stage_weights
-        )  # aligned with stages_after_warm_start
-        stage_weights = stage_weights / np.sum(stage_weights)
+        """
+        stage_weights = np.asarray(self.stage_weights)
+        stage_weights /= np.sum(stage_weights)  # Normalize weights
 
         n_candidates = self.num_candidate_training_points
         candidate_training_points_exploitation = []
@@ -544,26 +547,29 @@ def _get_exploitation_candidates(self):
         for stage_idx, stage_num in enumerate(self.stages_after_warm_start):
             samples_stage = self.results['model_parameters_dict'][stage_num]
             n_available = samples_stage.shape[0]
+
             n_stage_samples = int(np.round(stage_weights[stage_idx] * n_candidates))
+            if n_stage_samples == 0:
+                continue
 
-            if n_stage_samples >= n_available:
-                selected = samples_stage
-                n_selected = n_available
-            else:
+            if n_stage_samples <= n_available:
                 selected_indices = np.random.choice(
                     n_available, size=n_stage_samples, replace=False
                 )
                 selected = samples_stage[selected_indices]
-                n_selected = n_stage_samples
+            else:
+                # Use all available samples
+                selected = samples_stage
 
             candidate_training_points_exploitation.append(selected)
-            stage_sample_counts[stage_num] = n_selected
+            stage_sample_counts[stage_num] = selected.shape[0]
 
+        # Combine all stage samples
         candidate_training_points_exploitation = np.vstack(
             candidate_training_points_exploitation
         )
 
-        # Trim if needed
+        # Optional: globally truncate if total exceeds n_candidates (e.g., due to rounding)
         if candidate_training_points_exploitation.shape[0] > n_candidates:
             trim_indices = np.random.choice(
                 candidate_training_points_exploitation.shape[0],
@@ -573,9 +579,7 @@ def _get_exploitation_candidates(self):
             candidate_training_points_exploitation = (
                 candidate_training_points_exploitation[trim_indices]
             )
-        # self.loginfo(
-        #     'Completed candidate training point selection for exploitation.'
-        # )
+
         return candidate_training_points_exploitation, stage_sample_counts
 
     def _summarize_iteration(self):
@@ -801,8 +805,13 @@ def _step_1_posterior_approximation(
                 self.gp_prediction_mean, _ = self.current_gp_model.predict(  # type: ignore
                     model_parameters
                 )
+                # self.loo_predictions = self.current_gp_model.loo_predictions(  # type: ignore
+                #     model_parameters, model_outputs
+                # )
+
                 self.loo_predictions = self.current_gp_model.loo_predictions(  # type: ignore
-                    model_parameters, model_outputs
+                    self.current_gp_model.x_train,  # type: ignore
+                    self.current_gp_model.y_train,  # type: ignore
                 )
 
                 self.current_response_approximation = partial(
@@ -856,7 +865,9 @@ def _step_2_bayesian_updating(self, log_prior_fn, log_like_fn):
 
             with self.log_step('Evaluating warm-start of TMCMC.'):
                 if self.iteration_number > 0:
-                    weights = self.kde.evaluate(self.inputs)  # type: ignore
+                    weights = self.kde.evaluate(  # type: ignore
+                        self.current_gp_model.x_train  # type: ignore
+                    )  # inputs deduplicated in the gp model
                     self.gcv = self._calculate_gcv(weights)
                     self.warm_start_possible = self.gcv <= self.gcv_threshold
 
@@ -1058,7 +1069,7 @@ def _step_4_adaptive_training_point_selection(self):
                     (0, self.input_dimension)
                 )
 
-                current_inputs = self.inputs.copy()
+                current_inputs = self.current_gp_model.x_train.copy()  # type: ignore
                 if self.n_exploit > 0:
                     self.exploitation_training_points = (
                         current_doe.select_training_points(

modules/performUQ/common/gp_model.py

@@ -31,6 +31,31 @@
 # =========================================================
 
 
+def remove_duplicate_inputs(
+    x: np.ndarray, y: np.ndarray
+) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Remove duplicate rows in X and corresponding rows in Y.
+
+    Keeps only the first occurrence.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        Input array of shape (n_samples, n_features).
+    y : np.ndarray
+        Output array of shape (n_samples, n_outputs).
+
+    Returns
+    -------
+    x_unique, y_unique : tuple[np.ndarray, np.ndarray]
+        Deduplicated input-output arrays.
+    """
+    _, unique_indices = np.unique(x, axis=0, return_index=True)
+    sorted_indices = np.sort(unique_indices)
+    return x[sorted_indices], y[sorted_indices]
+
+
 class GaussianProcessModelSettings(BaseModel):
     """Settings for creating a Gaussian Process Model."""
 
@@ -227,8 +252,14 @@ def initialize(
         num_random_restarts : int, optional
             Number of random restarts for optimization (default is 10).
         """
-        self.x_train = x_train
-        self.y_train = y_train
+        orig_n = x_train.shape[0]
+        inputs, outputs = remove_duplicate_inputs(x_train, y_train)
+        deduped_n = inputs.shape[0]
+        if self.logger and deduped_n < orig_n:
+            self.logger.info(f'Removed {orig_n - deduped_n} duplicate input points.')
+
+        self.x_train = inputs
+        self.y_train = outputs
         self._fit_pca_and_create_models(
             reoptimize=reoptimize, num_random_restarts=num_random_restarts
         )
@@ -254,8 +285,14 @@ def update_training_dataset(
             msg = 'GP model not initialized. Call `initialize` first.'
             raise ValueError(msg)
 
-        self.x_train = x_train
-        self.y_train = y_train
+        orig_n = x_train.shape[0]
+        inputs, outputs = remove_duplicate_inputs(x_train, y_train)
+        deduped_n = inputs.shape[0]
+        if self.logger and deduped_n < orig_n:
+            self.logger.info(f'Removed {orig_n - deduped_n} duplicate input points.')
+
+        self.x_train = inputs
+        self.y_train = outputs
 
         self.x_train_scaled = self.apply_input_scaling(self.x_train, fit=True)