diff --git a/modules/common/groundMotionIM/IntensityMeasureComputer.py b/modules/common/groundMotionIM/IntensityMeasureComputer.py index 09d5911c3..2cf924861 100644 --- a/modules/common/groundMotionIM/IntensityMeasureComputer.py +++ b/modules/common/groundMotionIM/IntensityMeasureComputer.py @@ -343,13 +343,12 @@ def compute_response_spectrum(self, periods=[], damping=0.05, im_units=dict()): return elif type(periods) == list: # noqa: RET505, E721 periods = np.array(periods) - num_periods = len(periods) for cur_hist_name, cur_hist in self.time_hist_dict.items(): dt = cur_hist[1] ground_acc = cur_hist[2] num_steps = len(ground_acc) - # discritize + # discretize dt_disc = 0.005 num_steps_disc = int(np.floor(num_steps * dt / dt_disc)) f = interp1d( @@ -360,75 +359,73 @@ def compute_response_spectrum(self, periods=[], damping=0.05, im_units=dict()): ) tmp_time = [dt_disc * x for x in range(num_steps_disc)] ground_acc = f(tmp_time) + # circular frequency, damping, and stiffness terms - omega = (2 * np.pi) / periods - cval = damping * 2 * omega - kval = ((2 * np.pi) / periods) ** 2 + omega = (2.0 * np.pi) / periods + cval = damping * 2.0 * omega + kval = ((2.0 * np.pi) / periods) ** 2 + denom = 1.0 + 0.5 * dt_disc * cval + 0.25 * (dt_disc**2) * kval + + # state vectors (length = num_periods) + accel = np.zeros_like(omega) + vel = np.zeros_like(omega) + disp = np.zeros_like(omega) + + # rolling maxima for spectra + max_disp = np.zeros_like(omega) + max_vel = np.zeros_like(omega) + max_at = np.zeros_like(omega) + # Newmark-Beta - accel = np.zeros([num_steps_disc, num_periods]) - vel = np.zeros([num_steps_disc, num_periods]) - disp = np.zeros([num_steps_disc, num_periods]) - a_t = np.zeros([num_steps_disc, num_periods]) - accel[0, :] = (-ground_acc[0] - (cval * vel[0, :])) - (kval * disp[0, :]) + accel[:] = -ground_acc[0] - (cval * vel) - (kval * disp) for j in range(1, num_steps_disc): delta_acc = ground_acc[j] - ground_acc[j - 1] delta_d2u = ( -delta_acc - - dt_disc * cval * accel[j - 1, :] - - dt_disc - * kval - * (vel[j - 1, :] + 0.5 * dt_disc * accel[j - 1, :]) - ) / (1.0 + 0.5 * dt_disc * cval + 0.25 * dt_disc**2 * kval) - delta_du = dt_disc * accel[j - 1, :] + 0.5 * dt_disc * delta_d2u + - dt_disc * cval * accel + - dt_disc * kval * (vel + 0.5 * dt_disc * accel) + ) / denom + delta_du = dt_disc * accel + 0.5 * dt_disc * delta_d2u delta_u = ( - dt_disc * vel[j - 1, :] - + 0.5 * dt_disc**2 * accel[j - 1, :] + dt_disc * vel + + 0.5 * dt_disc**2 * accel + 0.25 * dt_disc**2 * delta_d2u ) - accel[j, :] = delta_d2u + accel[j - 1, :] - vel[j, :] = delta_du + vel[j - 1, :] - disp[j, :] = delta_u + disp[j - 1, :] - a_t[j, :] = ground_acc[j] + accel[j, :] + accel += delta_d2u + vel += delta_du + disp += delta_u + a_t = ground_acc[j] + accel + + # update rolling maxima in-place + np.maximum(np.fabs(disp), max_disp, out=max_disp) + np.maximum(np.fabs(vel), max_vel, out=max_vel) + np.maximum(np.fabs(a_t), max_at, out=max_at) + # collect data self.disp_spectrum.update( - { - cur_hist_name: np.ndarray.tolist( - unit_factor_psd * np.max(np.fabs(disp), axis=0) - ) - } + {cur_hist_name: np.ndarray.tolist(unit_factor_psd * max_disp)} ) self.vel_spectrum.update( - { - cur_hist_name: np.ndarray.tolist( - unit_factor_vspec * np.max(np.fabs(vel), axis=0) - ) - } + {cur_hist_name: np.ndarray.tolist(unit_factor_vspec * max_vel)} ) self.acc_spectrum.update( { cur_hist_name: np.ndarray.tolist( - unit_factor_aspec - * np.max(np.fabs(a_t), axis=0) - / 100.0 - / self.g + unit_factor_aspec * max_at / 100.0 / self.g ) } ) self.psv.update( { cur_hist_name: np.ndarray.tolist( - unit_factor_psv * omega * np.max(np.fabs(disp), axis=0) + unit_factor_psv * omega * max_disp ) } ) self.psa.update( { cur_hist_name: np.ndarray.tolist( - unit_factor_psa - * omega**2 - * np.max(np.fabs(disp), axis=0) - / 100.0 - / self.g + unit_factor_psa * omega**2 * max_disp / 100.0 / self.g ) } ) diff --git a/modules/performUQ/common/adaptive_doe.py b/modules/performUQ/common/adaptive_doe.py index 4c5c19a85..bc288cd93 100644 --- a/modules/performUQ/common/adaptive_doe.py +++ b/modules/performUQ/common/adaptive_doe.py @@ -6,6 +6,8 @@ from pathlib import Path import numpy as np +import plotly.graph_objects as go +import plotly.subplots as sp import uq_utilities from gp_model import remove_duplicate_inputs from scipy.stats import qmc @@ -72,7 +74,7 @@ def _compute_doe_weights(self) -> tuple[np.ndarray, dict]: 'basis': basis, 'basis_values': basis_vals.tolist(), 'sum': float(np.sum(w)), - 'count': int(len(w)), + 'count': len(w), } # no PCA: weight by per-output variance in the space the GP trained on @@ -84,7 +86,7 @@ def _compute_doe_weights(self) -> tuple[np.ndarray, dict]: 'basis': 'uniform_no_training_data', 'basis_values': [1.0] * n_models, 'sum': float(np.sum(w)), - 'count': int(len(w)), + 'count': len(w), } if self.gp_model.output_scaler is not None: @@ -117,7 +119,7 @@ def _compute_doe_weights(self) -> tuple[np.ndarray, dict]: 'basis': basis, 'basis_values': basis_vals.tolist(), 'sum': float(np.sum(w)), - 'count': int(len(w)), + 'count': len(w), } def _hyperparameters_for_doe(self): @@ -315,7 +317,7 @@ def write_gp_for_doe_to_json(self, filepath: str | Path): kernel = self.gp_model_for_doe.kern doe_hyperparams = { p.name: { - 'value': p.values.tolist() if p.size > 1 else float(p.values), # noqa: PD011 + 'value': p.values.tolist() if p.size > 1 else float(p.values), 'shape': p.shape, } for p in kernel.parameters @@ -340,6 +342,495 @@ def write_gp_for_doe_to_json(self, filepath: str | Path): with path.open('w') as f: json.dump(uq_utilities.make_json_serializable(output), f, indent=4) + def visualize_acquisition_analysis( # noqa: C901 + self, + x_train, + candidate_pool, + mci_samples, + *, + use_mse_w=True, + weights=None, + plot_dims=(0, 1), + n_grid=100, + show_before_after=True, + save_path=None, + ): + """ + Visualize acquisition function and other relevant fields. + + Visualize GP variance field, weight function, and acquisition function + for 2D projections of the design space using Plotly. + + Parameters + ---------- + x_train : array-like + Current training data (original space) + candidate_pool : array-like + Candidate points for selection + mci_samples : array-like + Monte Carlo integration samples + use_mse_w : bool + Whether to use MSEw (True) or IMSEw (False) + weights : array-like or None + Optional importance weights + plot_dims : tuple + Dimensions to plot (default: (0,1)) + n_grid : int + Grid resolution for plotting + show_before_after : bool + Whether to show before/after adding one point + save_path : str or None + Path to save figure (HTML format) + """ + # Ensure we have the DoE GP model + if not hasattr(self, 'gp_model_for_doe'): + msg = 'DoE GP model not initialized. Call select_training_points first.' + raise RuntimeError(msg) + + x_train = np.atleast_2d(x_train) + x_train, _ = remove_duplicate_inputs( + x_train, np.zeros((x_train.shape[0], 1)) + ) + + dim1, dim2 = plot_dims + + # Create plotting grid + x1_min = min(np.min(x_train[:, dim1]), np.min(candidate_pool[:, dim1])) + x1_max = max(np.max(x_train[:, dim1]), np.max(candidate_pool[:, dim1])) + x2_min = min(np.min(x_train[:, dim2]), np.min(candidate_pool[:, dim2])) + x2_max = max(np.max(x_train[:, dim2]), np.max(candidate_pool[:, dim2])) + + # Add padding + x1_range = x1_max - x1_min + x2_range = x2_max - x2_min + padding = 0.1 + x1_min -= padding * x1_range + x1_max += padding * x1_range + x2_min -= padding * x2_range + x2_max += padding * x2_range + + x1_grid = np.linspace(x1_min, x1_max, n_grid) + x2_grid = np.linspace(x2_min, x2_max, n_grid) + X1, X2 = np.meshgrid(x1_grid, x2_grid) + + # Create grid points (assuming other dimensions at mean of training data) + n_dims = x_train.shape[1] + grid_points = np.zeros((n_grid * n_grid, n_dims)) + grid_points[:, dim1] = X1.flatten() + grid_points[:, dim2] = X2.flatten() + + # Set other dimensions to mean of training data + for i in range(n_dims): + if i not in [dim1, dim2]: + grid_points[:, i] = np.mean(x_train[:, i]) + + def compute_fields(current_x_train): + """Compute variance and acquisition fields for given training set.""" + # Scale training data and grid + scaled_x_train = self._scale(current_x_train) + scaled_grid = self._scale(grid_points) + scaled_candidates = self._scale(candidate_pool) + scaled_mci = self._scale(mci_samples) + + # Set GP with current training set + self.gp_model_for_doe.set_XY( + scaled_x_train, np.zeros((scaled_x_train.shape[0], 1)) + ) + + # Compute variance field over grid + _, var_grid = self.gp_model_for_doe.predict(scaled_grid) + variance_field = var_grid.reshape(n_grid, n_grid) + + # Compute acquisition values for candidates + if use_mse_w: + # MSEw: direct variance at candidates + _, pred_var = self.gp_model_for_doe.predict(scaled_candidates) + if weights is not None: + pred_var *= weights.reshape(-1, 1) + acquisition_values = pred_var.flatten() + + # For plotting: acquisition field is just variance field + acquisition_field = variance_field.copy() + + else: + # IMSEw: integrated variance weighted by correlations + _, pred_var_mci = self.gp_model_for_doe.predict(scaled_mci) + if weights is not None: + pred_var_mci *= weights.reshape(-1, 1) + + # Compute acquisition for candidates + k_matrix = self.gp_model_for_doe.kern.K( + scaled_mci, scaled_candidates + ) + k_diag_mci = np.diag(self.gp_model_for_doe.kern.K(scaled_mci)) + k_diag_candidates = np.diag( + self.gp_model_for_doe.kern.K(scaled_candidates) + ) + + corr_matrix = k_matrix / np.sqrt( + np.outer(k_diag_mci, k_diag_candidates) + ) + beta = 2.0 * scaled_x_train.shape[1] + weighted_var = (corr_matrix**beta).T @ pred_var_mci + acquisition_values = (weighted_var / scaled_mci.shape[0]).flatten() + + # Compute acquisition field over grid + k_matrix_grid = self.gp_model_for_doe.kern.K(scaled_mci, scaled_grid) + k_diag_grid = np.diag(self.gp_model_for_doe.kern.K(scaled_grid)) + corr_matrix_grid = k_matrix_grid / np.sqrt( + np.outer(k_diag_mci, k_diag_grid) + ) + + weighted_var_grid = (corr_matrix_grid**beta).T @ pred_var_mci + acquisition_field = ( + weighted_var_grid / scaled_mci.shape[0] + ).reshape(n_grid, n_grid) + + return variance_field, acquisition_field, acquisition_values + + # Compute initial fields + var_field_before, acq_field_before, acq_values = compute_fields(x_train) + + # Select next point + best_idx = np.argmax(acq_values) + next_point = candidate_pool[best_idx] + + # Setup figure + if show_before_after: + fig = sp.make_subplots( + rows=2, + cols=3, + subplot_titles=[ + 'Variance Field (Before)', + 'Acquisition Field (Before)', + 'Combined (Before)', + 'Variance Field (After)', + 'Acquisition Field (After)', + 'Combined (After)', + ], + specs=[ + [{'type': 'xy'}, {'type': 'xy'}, {'type': 'xy'}], + [{'type': 'xy'}, {'type': 'xy'}, {'type': 'xy'}], + ], + ) + else: + fig = sp.make_subplots( + rows=1, + cols=3, + subplot_titles=['Variance Field', 'Acquisition Field', 'Combined'], + specs=[[{'type': 'xy'}, {'type': 'xy'}, {'type': 'xy'}]], + ) + + def add_field_plots(row, var_field, acq_field, title_suffix=''): + # Variance field + fig.add_trace( + go.Contour( + x=x1_grid, + y=x2_grid, + z=var_field, + colorscale='Viridis', + showscale=True, + colorbar=dict(title='Variance', x=0.32, len=0.4), # noqa: C408 + ), + row=row, + col=1, + ) + + # Training points + fig.add_trace( + go.Scatter( + x=x_train[:, dim1], + y=x_train[:, dim2], + mode='markers', + marker=dict(size=10, color='red', symbol='square'), # noqa: C408 + name='Training Points', + showlegend=(row == 1), + ), + row=row, + col=1, + ) + + # Candidates + fig.add_trace( + go.Scatter( + x=candidate_pool[:, dim1], + y=candidate_pool[:, dim2], + mode='markers', + marker=dict(size=4, color='orange', opacity=0.3), # noqa: C408 + name='Candidates', + showlegend=(row == 1), + ), + row=row, + col=1, + ) + + # Selected point (only for "before" plots) + if title_suffix == '': + fig.add_trace( + go.Scatter( + x=[next_point[dim1]], + y=[next_point[dim2]], + mode='markers', + marker=dict( # noqa: C408 + size=15, + color='white', + symbol='star', + line=dict(color='black', width=2), # noqa: C408 + ), + name='Selected Point', + showlegend=(row == 1), + ), + row=row, + col=1, + ) + + # Acquisition field + acq_name = 'MSEw' if use_mse_w else 'IMSEw' + fig.add_trace( + go.Contour( + x=x1_grid, + y=x2_grid, + z=acq_field, + colorscale='Plasma', + showscale=True, + colorbar=dict(title=f'{acq_name} Value', x=0.66, len=0.4), # noqa: C408 + ), + row=row, + col=2, + ) + + # Training points + fig.add_trace( + go.Scatter( + x=x_train[:, dim1], + y=x_train[:, dim2], + mode='markers', + marker=dict(size=10, color='red', symbol='square'), # noqa: C408 + showlegend=False, + ), + row=row, + col=2, + ) + + # Candidates + fig.add_trace( + go.Scatter( + x=candidate_pool[:, dim1], + y=candidate_pool[:, dim2], + mode='markers', + marker=dict(size=4, color='orange', opacity=0.3), # noqa: C408 + showlegend=False, + ), + row=row, + col=2, + ) + + # Selected point + if title_suffix == '': + fig.add_trace( + go.Scatter( + x=[next_point[dim1]], + y=[next_point[dim2]], + mode='markers', + marker=dict( # noqa: C408 + size=15, + color='white', + symbol='star', + line=dict(color='black', width=2), # noqa: C408 + ), + showlegend=False, + ), + row=row, + col=2, + ) + + # Combined view - acquisition field with variance contours + fig.add_trace( + go.Contour( + x=x1_grid, + y=x2_grid, + z=acq_field, + colorscale='Plasma', + opacity=0.8, + showscale=True, + colorbar=dict(title=f'{acq_name} Value', x=1.0, len=0.4), # noqa: C408 + ), + row=row, + col=3, + ) + + # Variance contours + fig.add_trace( + go.Contour( + x=x1_grid, + y=x2_grid, + z=var_field, + showscale=False, + line=dict(color='white', width=1), # noqa: C408 + opacity=0.5, + contours=dict( # noqa: C408 + showlabels=True, + labelfont=dict(size=8, color='white'), # noqa: C408 + ), + ), + row=row, + col=3, + ) + + # Training points + fig.add_trace( + go.Scatter( + x=x_train[:, dim1], + y=x_train[:, dim2], + mode='markers', + marker=dict(size=10, color='red', symbol='square'), # noqa: C408 + showlegend=False, + ), + row=row, + col=3, + ) + + # Selected point + if title_suffix == '': + fig.add_trace( + go.Scatter( + x=[next_point[dim1]], + y=[next_point[dim2]], + mode='markers', + marker=dict( # noqa: C408 + size=15, + color='white', + symbol='star', + line=dict(color='black', width=2), # noqa: C408 + ), + showlegend=False, + ), + row=row, + col=3, + ) + + # Plot before selection + add_field_plots(1, var_field_before, acq_field_before, '') + + if show_before_after: + # Compute fields after adding the selected point + x_train_after = np.vstack([x_train, next_point]) + var_field_after, acq_field_after, _ = compute_fields(x_train_after) + + # Plot after selection + add_field_plots(2, var_field_after, acq_field_after, ' (After)') + + # Update training points for after plots + for col in range(1, 4): + fig.add_trace( + go.Scatter( + x=x_train_after[:, dim1], + y=x_train_after[:, dim2], + mode='markers', + marker=dict(size=10, color='red', symbol='square'), # noqa: C408 + showlegend=False, + ), + row=2, + col=col, + ) + + # Update layout + fig.update_layout( + height=800 if show_before_after else 400, + width=1200, + title_text=f'DoE Acquisition Analysis - {("MSEw" if use_mse_w else "IMSEw")}', + showlegend=True, + ) + + # Update axes labels + for row in range(1, (3 if show_before_after else 2)): + for col in range(1, 4): + fig.update_xaxes( + title_text=f'Dimension {dim1 + 1}', row=row, col=col + ) + fig.update_yaxes( + title_text=f'Dimension {dim2 + 1}', row=row, col=col + ) + + # Add kernel information as annotation + kernel_info = f'Kernel: {self.gp_model_for_doe.kern.name}
' + if hasattr(self.gp_model_for_doe.kern, 'lengthscale'): + ls = self.gp_model_for_doe.kern.lengthscale.values + if self.ARD and len(ls) > 1: + kernel_info += f'Lengthscales: {[f"{l:.3f}" for l in ls]}
' + else: + kernel_info += f'Lengthscale: {ls[0]:.3f}
' + + if hasattr(self.gp_model_for_doe.kern, 'variance'): + kernel_info += ( + f'Variance: {self.gp_model_for_doe.kern.variance.values[0]:.3f}
' + ) + + kernel_info += ( + f'Nugget: {self.gp_model_for_doe.likelihood.variance.values[0]:.3e}
' + ) + kernel_info += f'Training points: {len(x_train)}
' + kernel_info += ( + f'Next point: [{next_point[dim1]:.3f}, {next_point[dim2]:.3f}]' + ) + + fig.add_annotation( + text=kernel_info, + xref='paper', + yref='paper', + x=0.02, + y=0.98, + showarrow=False, + align='left', + bgcolor='wheat', + bordercolor='black', + borderwidth=1, + font=dict(size=10), # noqa: C408 + ) + + if save_path: + fig.write_html(save_path) + print(f'Figure saved to {save_path}') + + fig.show() + + # Print diagnostic information + print('\nDIAGNOSTIC INFORMATION:') # noqa: T201 + print(f'Acquisition method: {"MSEw" if use_mse_w else "IMSEw"}') # noqa: T201 + print(f'GP Kernel: {self.gp_model_for_doe.kern.name}') # noqa: T201 + + if hasattr(self.gp_model_for_doe.kern, 'lengthscale'): + ls = self.gp_model_for_doe.kern.lengthscale.values + print(f'Lengthscales: {ls}') # noqa: T201 + print(f'Min lengthscale: {np.min(ls):.4f}') # noqa: T201 + print(f'Max lengthscale: {np.max(ls):.4f}') # noqa: T201 + + # Check if lengthscales are very small (potential clustering cause) + ls_thresh = 0.1 + if np.any(ls < ls_thresh): + print( # noqa: T201 + 'WARNING: Very small lengthscales detected - may cause clustering!' + ) + + print( # noqa: T201 + f'Variance range in field: [{np.min(var_field_before):.4f}, {np.max(var_field_before):.4f}]' + ) + print( # noqa: T201 + f'Acquisition range: [{np.min(acq_values):.4f}, {np.max(acq_values):.4f}]' + ) + print(f'Selected point acquisition value: {acq_values[best_idx]:.4f}') # noqa: T201 + print(f'Number of candidates: {len(candidate_pool)}') # noqa: T201 + + return { + 'variance_field_before': var_field_before, + 'acquisition_field_before': acq_field_before, + 'selected_point': next_point, + 'selected_point_value': acq_values[best_idx], + 'kernel_lengthscales': self.gp_model_for_doe.kern.lengthscale.values + if hasattr(self.gp_model_for_doe.kern, 'lengthscale') + else None, + } + def generate_lhs_candidates(n_samples, input_bounds, seed=None): """ diff --git a/modules/performUQ/common/gp_ab_algorithm.py b/modules/performUQ/common/gp_ab_algorithm.py index e4e16c88e..f3453613e 100644 --- a/modules/performUQ/common/gp_ab_algorithm.py +++ b/modules/performUQ/common/gp_ab_algorithm.py @@ -1055,7 +1055,7 @@ def _step_4_adaptive_training_point_selection(self): self._get_exploitation_candidates() ) save_exploitation_candidates_by_stage_json( - Path('results') + res_dir / f'exploitation_candidates_{self.iteration_number}.json', candidates_exploit, self.stage_sample_counts, @@ -1099,6 +1099,12 @@ def _step_4_adaptive_training_point_selection(self): self.num_candidate_training_points ) ) + save_exploration_candidates_json( + res_dir / f'exploration_candidates_{self.iteration_number}.json', + candidates_explore, + self.iteration_number, + ) + self.exploration_training_points = np.empty( (0, self.input_dimension) ) @@ -1108,7 +1114,8 @@ def _step_4_adaptive_training_point_selection(self): current_inputs, self.n_explore, candidates_explore, - use_mse_w=False, + # use_mse_w=False, + use_mse_w=True, weights=None, ) ) @@ -1308,7 +1315,7 @@ def save_tabular_results( if length == 1: pred_headers.append(name) else: - pred_headers.extend([f'{name}_{i+1}' for i in range(length)]) + pred_headers.extend([f'{name}_{i + 1}' for i in range(length)]) # Step 2: Create base dataframes df_samples = pd.DataFrame(samples, columns=rv_names_list) @@ -1433,6 +1440,51 @@ def write_results(self, results_dir='results'): ) +def save_exploration_candidates_json( + out_file: Path, + samples: np.ndarray, + iteration: int, + generation_method: str = 'latin_hypercube_sampling', +): + """ + Save exploration candidate samples to a JSON file. + + Parameters + ---------- + out_file : Path + Path to the output JSON file. + samples : np.ndarray + Array of exploration candidate samples of shape (n_candidates, n_parameters). + iteration : int + The iteration number when these candidates were generated. + generation_method : str, optional + Method used to generate candidates (default: "latin_hypercube_sampling"). + + Notes + ----- + The JSON file will contain: + - 'exploration_candidates': list of candidate parameter vectors + - 'num_candidates': total number of candidates + - 'generation_method': method used to generate candidates + - 'iteration': iteration number + """ + output = { + 'exploration_candidates': samples.tolist(), + 'num_candidates': len(samples), + 'generation_method': generation_method, + 'iteration': iteration, + } + + output_json_serializable = uq_utilities.make_json_serializable(output) + + # Ensure parent directory exists + out_file.parent.mkdir(parents=True, exist_ok=True) + + # Save to JSON + with out_file.open('w') as f: + json.dump(output_json_serializable, f, indent=2) + + def save_exploitation_candidates_by_stage_json( out_file: Path, samples: np.ndarray, stage_sample_counts: dict[int, int] ): @@ -1671,7 +1723,7 @@ def preprocess(input_arguments): if count == 1: headings += f'{name} ' else: - headings += ' '.join(f'{name}_{i+1}' for i in range(count)) + ' ' + headings += ' '.join(f'{name}_{i + 1}' for i in range(count)) + ' ' f_out.write(headings.strip() + '\n') linenum = 0