From 40a86ff0adf0cde550e4bce6510cc5c4824fe66d Mon Sep 17 00:00:00 2001
From: bsaakash <11618528+bsaakash@users.noreply.github.com>
Date: Tue, 23 Sep 2025 11:05:35 -0700
Subject: [PATCH 1/3] abs - saving exploration candidates to json
---
modules/performUQ/common/gp_ab_algorithm.py | 60 +++++++++++++++++++--
1 file changed, 56 insertions(+), 4 deletions(-)
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
From f2ab32a1de95fe185bd0023921ede2ccd4ba3440 Mon Sep 17 00:00:00 2001
From: bsaakash <11618528+bsaakash@users.noreply.github.com>
Date: Tue, 23 Sep 2025 15:48:38 -0700
Subject: [PATCH 2/3] abs - method for visualizing the variance field and
acquisition function in doe
---
modules/performUQ/common/adaptive_doe.py | 499 ++++++++++++++++++++++-
1 file changed, 495 insertions(+), 4 deletions(-)
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):
"""
From f161cd6cc5fb2531a90bb00b303310ca924c8d17 Mon Sep 17 00:00:00 2001
From: bsaakash <11618528+bsaakash@users.noreply.github.com>
Date: Thu, 25 Sep 2025 18:42:08 -0700
Subject: [PATCH 3/3] abs - more memory efficient calculation of response
spectra
---
.../IntensityMeasureComputer.py | 81 +++++++++----------
1 file changed, 39 insertions(+), 42 deletions(-)
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
)
}
)