diff --git a/.gitignore b/.gitignore index 63342b5..8ad58ca 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,10 @@ build # other typical case file names *.vti *morphology.csv +*.png # vscode settings cache .vscode + +# macOS Finder metadata +.DS_Store diff --git a/examples/.DS_Store b/examples/.DS_Store deleted file mode 100644 index 55d96b2..0000000 Binary files a/examples/.DS_Store and /dev/null differ diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index 78ce13a..086c18b 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -13,102 +13,118 @@ from raptor.api import ( compute_morphology, compute_porosity, + compute_spectral_components, create_grid, create_melt_pool, create_path_vectors, + visualize, write_morphology, write_vtk, ) -from raptor.utilities import MeltPoolFilter +from raptor.utilities import ( + MeltPoolFilter, + plot_melt_pool_signal, + reconstruct_spectral_signal, +) -# ============================================================================= -# User configuration -# ============================================================================= -# These are the main parameters to vary in a standard AM optimization study. -# All lengths are in meters, speed is in m/s, power is in W, and angles are in -# degrees. - -# Process parameters -# Accepted by the API, but does not affect the result (see note below). +# 0. User defined process parameters LASER_POWER = 370.0 SCAN_SPEED = 1.7 HATCH_SPACING = 140.0e-6 LAYER_HEIGHT = 30.0e-6 SCAN_ROTATION = 67.0 -# Mean melt-pool geometry MELT_POOL_WIDTH = 148.0e-6 MELT_POOL_DEPTH = 118.0e-6 -MELT_POOL_HEIGHT = 59.2e-6 - -# Stochastic melt-pool model +MELT_POOL_HEIGHT = 30.0e-6 MELT_POOL_WIDTH_STD_DEV = 18.0e-6 MELT_POOL_LENGTH = 300.0e-6 -N_SPECTRAL_MODES = 50 -RANDOM_SEED = 42 -# Lamé-curve exponents: 1 = parabola, 2 = ellipse. HEIGHT_SHAPE_FACTOR = 1.0 DEPTH_SHAPE_FACTOR = 1.0 -# RVE and output controls RVE_MIN_POINT = np.array([0.0, 0.0, 0.0]) RVE_MAX_POINT = np.array([5.0e-4, 5.0e-4, 5.0e-4]) VOXEL_RESOLUTION = 5.0e-6 + VTK_OUTPUT = "rve.vti" MORPHOLOGY_OUTPUT = "rve_morphology.csv" +WIDTH_DATA_PLOT = "melt_pool_width_timeseries.png" +ENABLE_VISUALIZATION = False def build_melt_pool(): - """Create a stochastic melt pool from the top-level configuration.""" - np.random.seed(RANDOM_SEED) - - # Internal signal-generation settings. These control the numerical - # representation of the stochastic melt-pool history. - sampling_frequency = 250000.0 - time_series_duration = 0.08 + # Create melt pools from convolution filter - melt_pool_filter = MeltPoolFilter( + # Instantiate object + mp_filter = MeltPoolFilter( MELT_POOL_WIDTH, MELT_POOL_WIDTH_STD_DEV, SCAN_SPEED, - [sampling_frequency, time_series_duration], + VOXEL_RESOLUTION, + random_seed=42, + ) + + # Define physical scales + mp_filter.add_effect("melt_pool", [MELT_POOL_LENGTH, None, 1]) + + # Generate stochastic melt pool + mp_filter.initialize() + width_data = mp_filter.generate_fluctuations(1.0, mp_filter.n_points, mp_filter.t) + width_spectral = compute_spectral_components(width_data, tolerance=VOXEL_RESOLUTION) + reconstructed_width = reconstruct_spectral_signal(width_data[:, 0], width_spectral) + reconstruction_rmse = np.sqrt( + np.mean((width_data[:, 1] - reconstructed_width) ** 2) + ) + print( + "Melt-pool signal: " + f"duration={mp_filter.duration:.6e} s, " + f"samples={mp_filter.n_points}, " + f"n_modes={width_spectral.shape[0]}, " + f"reconstruction_rmse={reconstruction_rmse:.6e} m" + ) + plot_melt_pool_signal( + width_data, + width_spectral, + MELT_POOL_WIDTH, + MELT_POOL_WIDTH_STD_DEV, + WIDTH_DATA_PLOT, + value_label="Melt-pool width (µm)", ) - # MeltPoolFilter converts this length scale to a temporal frequency - # using SCAN_SPEED / MELT_POOL_LENGTH. - melt_pool_filter.add_effect("melt_pool", [MELT_POOL_LENGTH, None, 1.0]) - melt_pool_filter.initialize() - width_data = melt_pool_filter.generate_fluctuations(noise_scale=1.0) + # scale melt pool data by constant factor + depth_scale = MELT_POOL_DEPTH / MELT_POOL_WIDTH + height_scale = MELT_POOL_HEIGHT / MELT_POOL_WIDTH - # Depth and height reuse the width history, scaled to their requested means. + # assign shape to melt pool and cap (1 = parabola, 2 = ellipse) melt_pool_dict = { - "width": (width_data, N_SPECTRAL_MODES, 1.0, 1.0), + "width": (width_data, None, 1.0, 2.0), "depth": ( width_data, - N_SPECTRAL_MODES, - MELT_POOL_DEPTH / MELT_POOL_WIDTH, + None, + depth_scale, DEPTH_SHAPE_FACTOR, ), "height": ( width_data, - N_SPECTRAL_MODES, - MELT_POOL_HEIGHT / MELT_POOL_WIDTH, + None, + height_scale, HEIGHT_SHAPE_FACTOR, ), } - return create_melt_pool(melt_pool_dict, enable_random_phases=True) + return create_melt_pool( + melt_pool_dict, enable_random_phases=True, tolerance=VOXEL_RESOLUTION + ) def main(): + # 1. Create voxel grid for the representative volume element (RVE) bound_box = np.array([RVE_MIN_POINT, RVE_MAX_POINT]) grid = create_grid(VOXEL_RESOLUTION, bound_box=bound_box) - # LASER_POWER is accepted by the current scan-path API but is not used by - # the porosity calculation. An optimization must supply its relationship to - # melt-pool geometry through the parameters. + # 2. Create path vectors through the representative volume element (RVE) path_vectors = create_path_vectors( bound_box, LASER_POWER, @@ -117,14 +133,19 @@ def main(): LAYER_HEIGHT, SCAN_ROTATION, scan_extension=max(RVE_MAX_POINT - RVE_MIN_POINT), - extra_layers=0, + extra_layers=10, ) + # 3. Create melt pools from convolution filter melt_pool = build_melt_pool() + + # 4. Compute porosity using conic section / superellipse curves for melt pool mask porosity = compute_porosity(grid, path_vectors, melt_pool, jit_warmup=True) + # 5. Write porosity field to .VTI write_vtk(grid.origin, grid.resolution, porosity, VTK_OUTPUT) + # 6. Compute morphology morphology = compute_morphology( porosity, VOXEL_RESOLUTION, @@ -132,6 +153,10 @@ def main(): ) write_morphology(morphology, MORPHOLOGY_OUTPUT) + # 7. Optionally visualize using PyVista (requires a graphical display) + if ENABLE_VISUALIZATION: + visualize(VTK_OUTPUT) + if __name__ == "__main__": main() diff --git a/pyproject.toml b/pyproject.toml index e7b3eb9..bb9e0d4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,8 @@ classifiers = [ dependencies = [ "numpy", "numba", + "scipy", + "matplotlib", "PyYAML", "vtk", "scikit-image", diff --git a/src/raptor/api.py b/src/raptor/api.py index 2fb6518..c6d2ead 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -9,7 +9,7 @@ # https://github.com/ORNL-MDF/Raptor/LICENSE # ============================================================================= import time -from typing import List, Tuple, Optional, Dict, Any +from typing import Any, Dict, List, Optional import numpy as np import pandas as pd import vtk @@ -65,55 +65,121 @@ def create_path_vectors( return scan_path_builder.process_vectors() -def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: int) -> np.ndarray: +def compute_required_modes(data: np.ndarray, reconstruction_rmse: float) -> int: + """Return the minimum number of real Fourier modes for an RMSE target.""" + return compute_spectral_components(data, tolerance=reconstruction_rmse).shape[0] - dt = melt_pool_data[1, 0] - melt_pool_data[0, 0] - mode0 = melt_pool_data[:, 1].mean() - fft_resolution = np.fft.fft(melt_pool_data[:, 1]) - F = np.zeros_like(fft_resolution) - n_fft = len(fft_resolution) - if n_modes == 1: - spectral_array = np.array([[mode0, 0, 0]]) - else: - for i in range(1, n_modes): - F[i] = fft_resolution[i] - F[n_fft - i] = fft_resolution[n_fft - i] - frequencies = np.float64(1 / (dt * n_fft)) * np.arange( - n_modes, dtype=np.float64 +def compute_spectral_components( + melt_pool_data: np.ndarray, + n_modes: Optional[int] = None, + tolerance: Optional[float] = None, +) -> np.ndarray: + """Convert a uniformly sampled real signal to a sparse cosine expansion. + + ``n_modes`` includes the mean (DC) mode and retains the corresponding + low-frequency prefix. When ``tolerance`` is used, the smallest set of + Fourier bins whose discarded energy satisfies the requested reconstruction + RMSE is retained. + """ + data = np.asarray(melt_pool_data, dtype=np.float64) + if data.ndim != 2 or data.shape[1] != 2 or data.shape[0] < 2: + raise ValueError("melt_pool_data must have shape (n, 2) with n >= 2.") + if not np.isfinite(data).all(): + raise ValueError("melt_pool_data must contain only finite values.") + if (n_modes is None) == (tolerance is None): + raise ValueError("Provide exactly one of n_modes or tolerance.") + + time_values = data[:, 0] + signal = data[:, 1] + time_steps = np.diff(time_values) + dt = time_steps[0] + if dt <= 0.0 or not np.allclose(time_steps, dt, rtol=1.0e-7, atol=0.0): + raise ValueError("Time samples must be strictly increasing and uniform.") + + n_samples = signal.size + fft_values = np.fft.rfft(signal) + frequencies = np.fft.rfftfreq(n_samples, d=dt) + candidate_bins = np.arange(1, fft_values.size) + + # Parseval energy represented by each positive-frequency cosine. Interior + # rFFT bins represent a conjugate pair; the Nyquist bin does not. + energy_weights = np.full(fft_values.size, 2.0) + energy_weights[0] = 1.0 + if n_samples % 2 == 0: + energy_weights[-1] = 1.0 + if tolerance is not None: + if not np.isfinite(tolerance) or tolerance < 0.0: + raise ValueError("tolerance must be a finite, non-negative value.") + candidate_energy = ( + energy_weights[candidate_bins] * np.abs(fft_values[candidate_bins]) ** 2 ) - phases = np.float64(np.angle(F[:n_modes])) - amplitudes = np.float64(np.abs(F[:n_modes]) / n_fft) - spectral_array = np.vstack( - [ - np.array([mode0, 0, 0]), - np.vstack([amplitudes[1:], frequencies[1:], phases[1:]]).transpose(), - ] + energy_indices = np.argsort(candidate_energy)[::-1] + energy_order = candidate_bins[energy_indices] + allowed_energy = (tolerance * n_samples) ** 2 + ordered_energy = candidate_energy[energy_indices] + required_energy = ordered_energy.sum() - allowed_energy + retained_count = ( + 0 + if required_energy <= 0.0 + else np.searchsorted( + np.cumsum(ordered_energy), required_energy, side="left" + ) + + 1 ) - return np.float64(spectral_array) + selected_bins = energy_order[:retained_count] + else: + if not isinstance(n_modes, (int, np.integer)) or n_modes < 1: + raise ValueError("n_modes must be a positive integer.") + if n_modes > fft_values.size: + raise ValueError( + f"n_modes cannot exceed {fft_values.size} for this time series." + ) + selected_bins = candidate_bins[: n_modes - 1] + + # Frequency order is convenient for evaluation and deterministic output. + selected_bins = np.sort(selected_bins) + amplitudes = energy_weights[selected_bins] * np.abs(fft_values[selected_bins]) + amplitudes /= n_samples + phases = np.angle(fft_values[selected_bins]) + # FFT phases are relative to sample zero; account for an absolute time axis. + phases -= 2.0 * np.pi * frequencies[selected_bins] * time_values[0] + + spectral_array = np.empty((selected_bins.size + 1, 3), dtype=np.float64) + spectral_array[0] = (signal.mean(), 0.0, 0.0) + spectral_array[1:, 0] = amplitudes + spectral_array[1:, 1] = frequencies[selected_bins] + spectral_array[1:, 2] = phases + return spectral_array def create_melt_pool( - melt_pool_dict: Dict[str, Any], enable_random_phases: bool + melt_pool_dict: Dict[str, Any], + enable_random_phases: bool, + tolerance: Optional[float] = None, ) -> MeltPool: - processed_components: Dict[str, Tuple[np.ndarray, float]] = {} - max_modes = 0 - - # 1. Determine the maximum number of modes required. - for _, nmodes, _, _ in melt_pool_dict.values(): - max_modes = max(max_modes, nmodes) + processed_components: Dict[str, np.ndarray] = {} - # 2. Process each component into its spectral format + # Process every component before padding so results do not depend on + # dictionary insertion order. for key, (data, n_modes, scale, shape_factor) in melt_pool_dict.items(): # Option A: Input data is a raw time-series [time, value] if data.shape[1] == 2: - spectral_array = compute_spectral_components(data, n_modes) + component_tolerance = tolerance + if tolerance is not None and scale != 0.0: + component_tolerance = tolerance / abs(scale) + spectral_array = compute_spectral_components( + data, + n_modes=n_modes, + tolerance=component_tolerance if n_modes is None else None, + ) spectral_array[:, 0] *= scale # Option B: Input data is a spectral array [amplitude, frequency, phase] elif data.shape[1] == 3: spectral_array = data.copy() + spectral_array[:, 0] *= scale else: raise ValueError( @@ -121,16 +187,17 @@ def create_melt_pool( f"Must be [time, value] or [amplitude, frequency, phase]" ) - # Pad the array with zeros if it has fewer modes than the max. + processed_components[key] = np.asarray(spectral_array, dtype=np.float64) + + max_modes = max(array.shape[0] for array in processed_components.values()) + for key, spectral_array in processed_components.items(): current_modes = spectral_array.shape[0] - if current_modes < max_modes: + if current_modes != max_modes: pad_array = np.zeros( shape=(max_modes - current_modes, spectral_array.shape[1]), dtype=np.float64, ) - spectral_array = np.vstack([spectral_array, pad_array]) - - processed_components[key] = spectral_array + processed_components[key] = np.vstack([spectral_array, pad_array]) # 3. Create the MeltPool object width_oscillations = processed_components["width"] diff --git a/src/raptor/cli.py b/src/raptor/cli.py index 43642b8..2a0f1a0 100644 --- a/src/raptor/cli.py +++ b/src/raptor/cli.py @@ -95,9 +95,14 @@ def main() -> int: filepath = melt_pool_dict[key]["file_name"] scale = melt_pool_dict[key]["scale"] nmodes = int(melt_pool_dict[key]["nmodes"]) - print("Scaling the zeroth mode (mean) of {} spectral array.") + shape_factor = 2 if key == "width" else melt_pool_dict[key]["shape"] spec_array = read_data(filepath) - melt_pool_data[key] = (spec_array, scale, nmodes) + melt_pool_data[key] = ( + spec_array, + nmodes, + scale, + shape_factor, + ) except: print( "Error reading the specified {} {} data format.".format( @@ -165,8 +170,8 @@ def main() -> int: for key in melt_pool_data: print(f" {key} datatype: " + melt_pool_dict[key]["type"]) print(f" {key} path : " + melt_pool_dict[key]["file_name"]) - print(f" {key} scaling : {melt_pool_data[key][1]}") - print(f" {key} modes : {melt_pool_data[key][2]}") + print(f" {key} scaling : {melt_pool_data[key][2]}") + print(f" {key} modes : {melt_pool_data[key][1]}") if bounding_box is not None: print(" RVE data:") diff --git a/src/raptor/io.py b/src/raptor/io.py index 8077b59..fa272d5 100644 --- a/src/raptor/io.py +++ b/src/raptor/io.py @@ -25,7 +25,10 @@ def read_data(fname: str) -> np.ndarray: """ if not os.path.exists(fname): raise FileNotFoundError(f"Melt pool measurement file not found: {fname}") - return np.loadtxt(fname, delimiter=",", dtype=np.float32) + # Preserve timestamp precision. Casting finely spaced absolute times to + # float32 can make an otherwise uniform series appear non-uniform to FFT + # consumers such as ``compute_spectral_components``. + return np.loadtxt(fname, delimiter=",", dtype=np.float64) def read_scan_path(fname: str) -> List[PathVector]: diff --git a/src/raptor/utilities.py b/src/raptor/utilities.py index 345efa8..1e400e5 100644 --- a/src/raptor/utilities.py +++ b/src/raptor/utilities.py @@ -8,10 +8,14 @@ # For details, see the top-level LICENSE file at: # https://github.com/ORNL-MDF/Raptor/LICENSE # ============================================================================= +from pathlib import Path +from typing import List, Optional, Tuple, Union + +import matplotlib.pyplot as plt import numpy as np -from typing import List, Tuple from .structures import PathVector -from scipy.signal import lfilter, butter +from scipy.signal import butter, sosfilt +from scipy.stats import chi2 class ScanPathBuilder: @@ -221,28 +225,34 @@ def write_layers(self, output_name, mode="layers"): class MeltPoolFilter: def __init__( - self, mu: float, sigma: float, scan_speed: float, timeseries_params: list + self, + mu: float, + sigma: float, + scan_speed: float, + voxel_resolution: float, + *, + confidence: float = 0.95, + ci_relative_width: float = 0.10, + correlation_tolerance: float = 0.10, + random_seed: Optional[int] = None, ): + """Generate a Gaussian melt-pool history from physical length scales. + + Sampling is tied to the spatial grid through ``fs = scan_speed / + voxel_resolution``. ``initialize`` chooses the shortest duration that + satisfies the requested variance precision while resolving each + effect's characteristic wavelength. ``correlation_tolerance`` controls + the relative wavelength-resolution check; it is not the FFT RMSE. """ - Filtration of disparate fluctuation scales to infer a melt pool oscillations sequence. - Uses scipy.butter to convolve scales of fluctuations together. - Initialization parameters: - mu: mean melt pool dimension - sigma: target standard deviation of the fluctuations - scan_speed: speed in m/s - timeseries_params: list of [fs,duration] - """ - # statistical properties self.mu, self.sigma = mu, sigma - # process parameters self.scan_speed = scan_speed - # timeseries related properties - self.fs, self.duration = timeseries_params - self.dt = 1 / self.fs - self.n_points = int(self.duration // self.dt + 1) - self.t = np.linspace(0, self.duration, self.n_points) - # parametric representations of fluctuation scales - self.physical_effects = {} # contains scale description and parameters + self.voxel_resolution = voxel_resolution + self.fs = self.scan_speed / self.voxel_resolution + self.confidence = confidence + self.ci_relative_width = ci_relative_width + self.correlation_tolerance = correlation_tolerance + self.rng = np.random.default_rng(random_seed) + self.physical_effects = {} def add_effect(self, effect_name: str, effect_params: list): """ @@ -257,17 +267,27 @@ def add_effect(self, effect_name: str, effect_params: list): "sigma_weight": sigma_weight, } - def initialize(self): + def initialize(self) -> None: + """Choose the minimum planned duration. + + Duration selection is deterministic. It uses the autocorrelation of + the configured linear filters to approximate the effective degrees of + freedom of a Gaussian-process variance estimate. + """ # Calculate frequencies from length scales - for name, params in self.physical_effects.items(): + for params in self.physical_effects.values(): if params["length_scale_m"] is not None: params["frequency_hz"] = self.scan_speed / params["length_scale_m"] - # Check for Nyquist limit violations - max_freq = max(p["frequency_hz"] for p in self.physical_effects.values()) - if max_freq > self.fs / 2: + # The complete passband, rather than only its center, must satisfy + # Nyquist. The same limits are used by bandpass_filter below. + highest_passband_frequency = max( + 1.5 * p["frequency_hz"] for p in self.physical_effects.values() + ) + if highest_passband_frequency >= self.fs / 2.0: raise ValueError( - f"Error: Maximum frequency ({max_freq/1000:.1f} kHz) exceeds Nyquist limit ({self.fs/2000:.1f} kHz). Increase sampling rate 'fs'." + "The filter passband exceeds the Nyquist frequency. " + "Decrease voxel_resolution or increase the physical length scale." ) # Normalize sigma weights so the variances sum correctly @@ -275,30 +295,115 @@ def initialize(self): sum_of_sq_weights = np.sum(weights**2) self.normalization_factor = np.sqrt(sum_of_sq_weights) - for name, params in self.physical_effects.items(): + for params in self.physical_effects.values(): params["sigma_contribution"] = ( params["sigma_weight"] / self.normalization_factor ) * self.sigma + max_timescale = max( + 1.0 / p["frequency_hz"] for p in self.physical_effects.values() + ) + minimum_periods = max(4, int(np.ceil(1.0 / self.correlation_tolerance))) + minimum_duration = minimum_periods * max_timescale + lower_points = max(3, int(np.ceil(minimum_duration * self.fs)) + 1) + + # The IIR response decays over a small number of characteristic periods; + # fifty periods provides a conservative autocorrelation horizon without + # making it scale with the final time-series duration. + autocorrelation = self._model_autocorrelation() + + upper_points = lower_points + for _ in range(16): + effective_n = self._effective_sample_size(upper_points, autocorrelation) + planned_ci = self._variance_ci(self.sigma**2, effective_n) + if planned_ci["precision_satisfied"]: + break + upper_points = 2 * upper_points - 1 + else: + raise RuntimeError( + "Unable to satisfy the variance precision within " + "16 duration refinements." + ) + + # Precision improves monotonically with sample count for the fixed model + # autocorrelation, so binary search gives the minimum accepted count. + left, right = lower_points, upper_points + while left < right: + midpoint = (left + right) // 2 + effective_n = self._effective_sample_size(midpoint, autocorrelation) + if self._variance_ci(self.sigma**2, effective_n)["precision_satisfied"]: + right = midpoint + else: + left = midpoint + 1 + + self.n_points = left + self.duration = (self.n_points - 1) / self.fs + self.t = np.arange(self.n_points, dtype=np.float64) / self.fs + self.planned_effective_n = self._effective_sample_size( + self.n_points, autocorrelation + ) + self.planned_variance_ci = self._variance_ci( + self.sigma**2, self.planned_effective_n + ) + + def _model_autocorrelation(self) -> np.ndarray: + """Return the normalized autocorrelation implied by the shared driver.""" + minimum_frequency = min( + params["frequency_hz"] for params in self.physical_effects.values() + ) + response_points = max(256, int(np.ceil(50.0 * self.fs / minimum_frequency))) + impulse = np.zeros(response_points, dtype=np.float64) + impulse[0] = 1.0 + combined_response = np.zeros(response_points, dtype=np.float64) + + for params in self.physical_effects.values(): + response = self.bandpass_filter( + impulse, + params["frequency_hz"], + 1.0, + self.fs, + ) + response_norm = np.sqrt(np.sum(response**2)) + combined_response += params["sigma_contribution"] * response / response_norm + + n_fft = 1 << (2 * response_points - 1).bit_length() + response_fft = np.fft.rfft(combined_response, n=n_fft) + autocovariance = np.fft.irfft( + response_fft * np.conjugate(response_fft), n=n_fft + )[:response_points] + return autocovariance / autocovariance[0] + + @staticmethod + def _effective_sample_size(n_samples: int, autocorrelation: np.ndarray) -> float: + """Approximate effective sample count for a Gaussian variance estimate.""" + max_lag = min(n_samples - 1, autocorrelation.size - 1) + if max_lag < 1: + return float(n_samples) + lags = np.arange(1, max_lag + 1, dtype=np.float64) + finite_sample_weights = 1.0 - lags / n_samples + correlation_penalty = 1.0 + 2.0 * np.sum( + finite_sample_weights * autocorrelation[1 : max_lag + 1] ** 2 + ) + return float(np.clip(n_samples / correlation_penalty, 2.0, n_samples)) + def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4): """Applies a bandpass filter around a center frequency f0.""" - lowcut = f0 * (1 - bandwidth_fraction / 2) - highcut = f0 * (1 + bandwidth_fraction / 2) + lowcut = f0 * (1.0 - bandwidth_fraction / 2.0) + highcut = f0 * (1.0 + bandwidth_fraction / 2.0) nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq - if high >= 1.0: - high = 0.999 - if low <= 0.0001: - low = 0.0001 - b, a = butter(order, [low, high], btype="band") - return lfilter(b, a, data) - - def generate_fluctuations(self, noise_scale): - base_white_noise = np.random.normal( - loc=0, scale=noise_scale, size=self.n_points - ) - final_series = np.zeros(self.n_points) + if not 0.0 < low < high < 1.0: + raise ValueError( + "Bandpass range must lie strictly between 0 Hz and the " + "Nyquist frequency. Increase the sampling frequency." + ) + sos = butter(order, [low, high], btype="band", output="sos") + return sosfilt(sos, data) + + def generate_fluctuations(self, noise_scale, n_points, t): + base_white_noise = self.rng.normal(loc=0.0, scale=noise_scale, size=n_points) + final_series = np.zeros(n_points) self.component_series = {} # Create each component series, scale it, and add to the final series @@ -310,14 +415,204 @@ def generate_fluctuations(self, noise_scale): fs=self.fs, ) + component_noise -= np.mean(component_noise) std_dev = np.std(component_noise) + if not np.isfinite(std_dev) or std_dev == 0.0: + raise ValueError("Unable to generate finite melt-pool fluctuations.") scaled_component = component_noise * ( params["sigma_contribution"] / std_dev ) self.component_series[name] = scaled_component final_series += scaled_component + # Shared noise can correlate filtered effects, so normalizing their + # individual variances does not generally normalize the variance of + # their sum. Apply one common correction to preserve relative weights. + aggregate_std = np.std(final_series) + if not np.isfinite(aggregate_std) or aggregate_std == 0.0: + raise ValueError("Unable to generate finite melt-pool fluctuations.") + covariance_scale = self.sigma / aggregate_std + final_series *= covariance_scale + for name in self.component_series: + self.component_series[name] *= covariance_scale + # Adding the mean final_series += self.mu - return np.column_stack([self.t, final_series]) + return np.column_stack([t, final_series]) + + def evaluate_variance_ci(self, data): + """ + Evaluates the confidence interval for the variance of the data. + Returns (lower_bound, upper_bound) for the variance. + """ + n = len(data) + sample_variance = np.var(data, ddof=1) + + # Compute autocorrelation to estimate effective sample size + autocorr = np.correlate(data - np.mean(data), data - np.mean(data), mode="full") + autocorr = autocorr[autocorr.size // 2 :] / autocorr[autocorr.size // 2] + + # Effective sample size + zero_crossings = np.where(autocorr < 0)[0] + cutoff = zero_crossings[0] if zero_crossings.size > 0 else len(autocorr) + effective_n = n / (1 + 2 * np.sum(autocorr[1:cutoff])) + effective_n = int(max(1, min(effective_n, n))) + + return self._variance_ci(sample_variance, effective_n) + + def _variance_ci(self, sample_variance: float, effective_n: float) -> dict: + """Construct the approximate Gaussian-process variance interval.""" + dof = effective_n - 1.0 + alpha = 1 - self.confidence + chi2_lower = chi2.ppf(alpha / 2, dof) + chi2_upper = chi2.ppf(1 - alpha / 2, dof) + + lower_bound = dof * sample_variance / chi2_upper + upper_bound = dof * sample_variance / chi2_lower + + target_within_ci = lower_bound <= self.sigma**2 <= upper_bound + target_std_lower = self.sigma * (1.0 - self.ci_relative_width) + target_std_upper = self.sigma * (1.0 + self.ci_relative_width) + precision_satisfied = ( + np.sqrt(lower_bound) >= target_std_lower + and np.sqrt(upper_bound) <= target_std_upper + ) + + return { + "sample_variance": sample_variance, + "effective_n": effective_n, + "lower_bound": lower_bound, + "upper_bound": upper_bound, + "target_within_ci": target_within_ci, + "precision_satisfied": precision_satisfied, + } + + +def reconstruct_spectral_signal( + time_values: np.ndarray, spectral_components: np.ndarray +) -> np.ndarray: + """Evaluate a ``[amplitude, frequency, phase]`` cosine expansion.""" + time_values = np.asarray(time_values, dtype=np.float64) + spectral_components = np.asarray(spectral_components, dtype=np.float64) + if spectral_components.ndim != 2 or spectral_components.shape[1] != 3: + raise ValueError("spectral_components must have shape (n_modes, 3).") + reconstructed = np.zeros_like(time_values) + for amplitude, frequency, phase in spectral_components: + reconstructed += amplitude * np.cos( + 2.0 * np.pi * frequency * time_values + phase + ) + return reconstructed + + +def plot_melt_pool_signal( + time_series: np.ndarray, + spectral_components: np.ndarray, + mean: float, + standard_deviation: float, + output_path: Union[str, Path], + *, + time_scale: float = 1.0e3, + value_scale: float = 1.0e6, + time_label: str = "Time (ms)", + value_label: str = "Melt-pool dimension (µm)", + bins: int = 40, +) -> Path: + """Write a 6.5-by-3 inch, 300 dpi signal and distribution figure.""" + data = np.asarray(time_series, dtype=np.float64) + if data.ndim != 2 or data.shape[1] != 2: + raise ValueError("time_series must have shape (n, 2).") + if standard_deviation <= 0.0: + raise ValueError("standard_deviation must be positive.") + + time_values = data[:, 0] + signal = data[:, 1] + reconstructed = reconstruct_spectral_signal(time_values, spectral_components) + plot_time = time_values * time_scale + plot_signal = signal * value_scale + plot_reconstructed = reconstructed * value_scale + plot_mean = mean * value_scale + plot_std = standard_deviation * value_scale + + distribution_limits = ( + min(plot_signal.min(), plot_mean - 4.0 * plot_std), + max(plot_signal.max(), plot_mean + 4.0 * plot_std), + ) + gaussian_values = np.linspace(*distribution_limits, 500) + gaussian_density = np.exp( + -0.5 * ((gaussian_values - plot_mean) / plot_std) ** 2 + ) / (plot_std * np.sqrt(2.0 * np.pi)) + + style = { + "font.size": 8, + "axes.labelsize": 8, + "legend.fontsize": 7, + "xtick.labelsize": 7, + "ytick.labelsize": 7, + "lines.linewidth": 1.0, + } + with plt.rc_context(style): + figure, (signal_axis, distribution_axis) = plt.subplots( + 1, 2, figsize=(6.5, 3.0), constrained_layout=True + ) + signal_axis.plot( + plot_time, plot_signal, color="#0072B2", linewidth=0.9, label="Generated" + ) + signal_axis.plot( + plot_time, + plot_reconstructed, + color="#D55E00", + linewidth=0.9, + linestyle="--", + label=f"Truncated FFT ({len(spectral_components)} modes)", + ) + signal_axis.set_xlabel(time_label) + signal_axis.set_ylabel(value_label) + signal_axis.legend(frameon=False, loc="upper right") + + distribution_axis.hist( + plot_signal, + bins=bins, + density=True, + color="#56B4E9", + edgecolor="white", + linewidth=0.4, + alpha=0.8, + label="Generated samples", + ) + distribution_axis.plot( + gaussian_values, + gaussian_density, + color="#D55E00", + linewidth=1.8, + label=f"Gaussian (µ={plot_mean:.0f}, σ={plot_std:.0f})", + ) + distribution_axis.axvline( + plot_mean, + color="black", + linewidth=0.8, + linestyle=":", + label="Specified mean", + ) + distribution_axis.set_xlim(distribution_limits) + distribution_axis.set_xlabel(value_label) + distribution_axis.set_ylabel("Probability density") + distribution_axis.legend(frameon=False, loc="upper right") + + for panel, axis in zip(("(a)", "(b)"), (signal_axis, distribution_axis)): + axis.text( + 0.02, + 0.96, + panel, + transform=axis.transAxes, + fontweight="bold", + va="top", + ) + axis.spines["top"].set_visible(False) + axis.spines["right"].set_visible(False) + axis.tick_params(direction="out") + + output = Path(output_path) + figure.savefig(output, dpi=300, facecolor="white") + plt.close(figure) + return output diff --git a/tests/test_api.py b/tests/test_api.py index df22e9c..4f8aed0 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -28,6 +28,7 @@ from raptor.api import ( create_grid, create_path_vectors, + compute_required_modes, compute_spectral_components, create_melt_pool, compute_porosity, @@ -88,7 +89,7 @@ def sample_process_parameters(): @pytest.fixture def sample_time_series_data(): """Fixture providing sample time series data for melt pool.""" - t = np.linspace(0, 1, 100) + t = np.arange(100, dtype=np.float64) / 100.0 values = 0.0001 + 0.00002 * np.sin(2 * np.pi * 5 * t) return np.column_stack([t, values]) @@ -332,12 +333,73 @@ def test_compute_spectral_components_mean_value(self, sample_time_series_data): np.testing.assert_allclose(spectral_array[0, 0], expected_mean) + def test_compute_spectral_components_preserves_standard_deviation( + self, sample_time_series_data + ): + """Test correct amplitude when the signal frequency is retained.""" + spectral_array = compute_spectral_components(sample_time_series_data, 6) + time = sample_time_series_data[:, 0] + reconstructed = np.zeros_like(time) + + for amplitude, frequency, phase in spectral_array: + reconstructed += amplitude * np.cos(2.0 * np.pi * frequency * time + phase) + + assert reconstructed.mean() == pytest.approx( + sample_time_series_data[:, 1].mean() + ) + assert reconstructed.std() == pytest.approx(sample_time_series_data[:, 1].std()) + + def test_explicit_mode_count_retains_low_frequency_prefix(self): + """Explicit mode counts preserve the original truncation semantics.""" + sampling_frequency = 10_000.0 + time = np.arange(1000) / sampling_frequency + values = 1.0 + 0.2 * np.sin(2.0 * np.pi * 2000.0 * time) + spectral_array = compute_spectral_components(np.column_stack([time, values]), 2) + + assert spectral_array[1, 1] == pytest.approx(10.0) + + def test_tolerance_selects_minimum_modes_and_meets_rmse(self): + sampling_frequency = 1000.0 + time = np.arange(1000) / sampling_frequency + values = ( + 2.0 + + 0.30 * np.cos(2.0 * np.pi * 20.0 * time + 0.2) + + 0.04 * np.cos(2.0 * np.pi * 170.0 * time - 0.4) + ) + data = np.column_stack([time, values]) + + spectral_array = compute_spectral_components(data, tolerance=0.03) + reconstructed = sum( + amplitude * np.cos(2.0 * np.pi * frequency * time + phase) + for amplitude, frequency, phase in spectral_array + ) + + assert spectral_array.shape == (2, 3) + assert spectral_array[1, 1] == pytest.approx(20.0) + assert np.sqrt(np.mean((values - reconstructed) ** 2)) <= 0.03 + assert compute_required_modes(data, 0.03) == 2 + + def test_tolerance_zero_reconstructs_even_length_signal(self): + time = 0.25 + np.arange(100) / 100.0 + values = 3.0 + 0.2 * np.cos(2.0 * np.pi * 5.0 * time) + values += 0.1 * np.cos(2.0 * np.pi * 50.0 * time) + spectral_array = compute_spectral_components( + np.column_stack([time, values]), tolerance=0.0 + ) + reconstructed = sum( + amplitude * np.cos(2.0 * np.pi * frequency * time + phase) + for amplitude, frequency, phase in spectral_array + ) + np.testing.assert_allclose(reconstructed, values, atol=1.0e-13) + def test_compute_spectral_components_invalid_input(self): """Test spectral component computation with invalid input.""" - with pytest.raises(IndexError): + with pytest.raises(ValueError, match="shape"): compute_spectral_components(np.array([[0.0, 1.0]]), 1) - with pytest.raises(IndexError): + with pytest.raises(ValueError, match="shape"): compute_spectral_components(np.ones((4, 1)), 2) + with pytest.raises(ValueError, match="exactly one"): + compute_spectral_components(np.ones((4, 2))) # ============================================================================= @@ -373,6 +435,24 @@ def test_create_melt_pool_spectral_input(self, sample_spectral_components): assert isinstance(melt_pool, MeltPool) + def test_create_melt_pool_scales_spectral_amplitudes( + self, sample_spectral_components + ): + scale_factor = 2.0 + melt_pool = create_melt_pool( + { + "width": (sample_spectral_components, 3, scale_factor, 2.0), + "depth": (sample_spectral_components, 3, 1.0, 2.0), + "height": (sample_spectral_components, 3, 1.0, 2.0), + }, + enable_random_phases=False, + ) + + np.testing.assert_allclose( + melt_pool.width_oscillations[:, 0], + scale_factor * sample_spectral_components[:, 0], + ) + def test_create_melt_pool_mode_padding(self): """Test that melt pool correctly pads modes to match maximum.""" one_mode = np.array([[1.0e-4, 0.0, 0.0]]) @@ -381,8 +461,8 @@ def test_create_melt_pool_mode_padding(self): ) melt_pool = create_melt_pool( { - "width": (three_modes, 3, 1.0, 2.0), - "depth": (one_mode, 1, 1.0, 2.0), + "width": (one_mode, 1, 1.0, 2.0), + "depth": (three_modes, 3, 1.0, 2.0), "height": (one_mode, 1, 1.0, 2.0), }, enable_random_phases=False, @@ -390,7 +470,7 @@ def test_create_melt_pool_mode_padding(self): assert melt_pool.width_oscillations.shape == (3, 3) assert melt_pool.depth_oscillations.shape == (3, 3) assert melt_pool.height_oscillations.shape == (3, 3) - np.testing.assert_array_equal(melt_pool.depth_oscillations[1:], 0.0) + np.testing.assert_array_equal(melt_pool.width_oscillations[1:], 0.0) def test_create_melt_pool_scaling(self, sample_time_series_data): """Test that scaling is correctly applied.""" @@ -410,6 +490,21 @@ def test_create_melt_pool_scaling(self, sample_time_series_data): sample_time_series_data[:, 1].mean() ) + def test_create_melt_pool_allows_fixed_and_tolerance_mode_counts( + self, sample_time_series_data + ): + melt_pool = create_melt_pool( + { + "width": (sample_time_series_data, 2, 1.0, 2.0), + "depth": (sample_time_series_data, None, 1.0, 2.0), + "height": (sample_time_series_data, None, 1.0, 2.0), + }, + enable_random_phases=False, + tolerance=1.0e-7, + ) + assert melt_pool.width_oscillations.shape == melt_pool.depth_oscillations.shape + assert melt_pool.depth_oscillations.shape == melt_pool.height_oscillations.shape + def test_create_melt_pool_shape_factors(self, sample_melt_pool_dict): """Test that shape factors are correctly set.""" melt_pool = create_melt_pool(sample_melt_pool_dict, enable_random_phases=False) diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..b2d668d --- /dev/null +++ b/tests/test_io.py @@ -0,0 +1,27 @@ +# ============================================================================= +# Copyright (c) 2025 Oak Ridge National Laboratory +# +# All rights reserved. +# +# This file is part of Raptor. +# +# For details, see the top-level LICENSE file at: +# https://github.com/ORNL-MDF/Raptor/LICENSE +# ============================================================================= +import numpy as np + +from raptor.api import compute_spectral_components +from raptor.io import read_data + + +def test_read_data_preserves_uniform_timestamp_precision(tmp_path): + time = np.arange(30_000, dtype=np.float64) * (1.0 / 6_137_000.0) + values = 1.0e-4 + 1.0e-5 * np.cos(2.0 * np.pi * 100.0 * time) + input_path = tmp_path / "melt_pool_data.csv" + np.savetxt(input_path, np.column_stack([time, values]), delimiter=",") + + data = read_data(input_path) + spectral_components = compute_spectral_components(data, n_modes=2) + + assert data.dtype == np.float64 + assert spectral_components.shape == (2, 3) diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 0000000..36be2b2 --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,136 @@ +# ============================================================================= +# Copyright (c) 2025 Oak Ridge National Laboratory +# +# All rights reserved. +# +# This file is part of Raptor. +# +# For details, see the top-level LICENSE file at: +# https://github.com/ORNL-MDF/Raptor/LICENSE +# ============================================================================= +import matplotlib.image as mpimg +import numpy as np + +from raptor.utilities import ( + MeltPoolFilter, + plot_melt_pool_signal, + reconstruct_spectral_signal, +) + + +def make_filter(*, voxel_resolution=5.0e-6, random_seed=None): + melt_pool_filter = MeltPoolFilter( + mu=148.0e-6, + sigma=40.0e-6, + scan_speed=1.7, + confidence=0.95, + ci_relative_width=0.30, + voxel_resolution=voxel_resolution, + correlation_tolerance=0.20, + random_seed=random_seed, + ) + melt_pool_filter.add_effect("melt_pool", [300.0e-6, None, 1.0]) + return melt_pool_filter + + +def test_melt_pool_filter_recovers_requested_statistics(): + melt_pool_filter = make_filter() + melt_pool_filter.initialize() + width_data = melt_pool_filter.generate_fluctuations( + 1.0, melt_pool_filter.n_points, melt_pool_filter.t + ) + + assert np.isfinite(width_data).all() + np.testing.assert_allclose(width_data[:, 1].mean(), melt_pool_filter.mu, atol=1e-15) + np.testing.assert_allclose( + width_data[:, 1].std(), melt_pool_filter.sigma, atol=1e-15 + ) + assert melt_pool_filter.planned_variance_ci["precision_satisfied"] + + +def test_melt_pool_filter_normalizes_covariance_between_effects(): + melt_pool_filter = make_filter() + melt_pool_filter.add_effect("overlapping", [300.0e-6, None, 1.0]) + melt_pool_filter.initialize() + width_data = melt_pool_filter.generate_fluctuations( + 1.0, melt_pool_filter.n_points, melt_pool_filter.t + ) + + np.testing.assert_allclose(width_data[:, 1].mean(), melt_pool_filter.mu, atol=1e-15) + np.testing.assert_allclose( + width_data[:, 1].std(), melt_pool_filter.sigma, atol=1e-15 + ) + + +def test_melt_pool_filter_seed_reproduces_filtered_time_series(): + first_filter = make_filter(random_seed=42) + second_filter = make_filter(random_seed=42) + first_filter.initialize() + second_filter.initialize() + + first = first_filter.generate_fluctuations( + 1.0, first_filter.n_points, first_filter.t + ) + second = second_filter.generate_fluctuations( + 1.0, second_filter.n_points, second_filter.t + ) + + np.testing.assert_array_equal(first, second) + + +def test_duration_is_the_minimum_accepted_sample_count(): + melt_pool_filter = make_filter() + melt_pool_filter.initialize() + model_correlation = melt_pool_filter._model_autocorrelation() + + previous_n = melt_pool_filter.n_points - 1 + previous_effective_n = melt_pool_filter._effective_sample_size( + previous_n, model_correlation + ) + previous_ci = melt_pool_filter._variance_ci( + melt_pool_filter.sigma**2, previous_effective_n + ) + + minimum_periods = max(4, int(np.ceil(1.0 / melt_pool_filter.correlation_tolerance))) + minimum_n = ( + int( + np.ceil( + minimum_periods + * (300.0e-6 / melt_pool_filter.scan_speed) + * melt_pool_filter.fs + ) + ) + + 1 + ) + if previous_n >= minimum_n: + assert not previous_ci["precision_satisfied"] + + +def test_melt_pool_filter_uses_statistical_defaults(): + melt_pool_filter = MeltPoolFilter(148.0e-6, 40.0e-6, 1.7, 5.0e-6) + assert melt_pool_filter.confidence == 0.95 + assert melt_pool_filter.ci_relative_width == 0.10 + assert melt_pool_filter.correlation_tolerance == 0.10 + + +def test_reconstruct_spectral_signal(): + time = np.arange(100) / 100.0 + components = np.array([[2.0, 0.0, 0.0], [0.5, 4.0, np.pi / 3.0]]) + expected = 2.0 + 0.5 * np.cos(2.0 * np.pi * 4.0 * time + np.pi / 3.0) + np.testing.assert_allclose( + reconstruct_spectral_signal(time, components), expected, atol=1e-14 + ) + + +def test_publication_plot_has_requested_pixel_dimensions(tmp_path): + time = np.arange(500) / 10_000.0 + values = 148.0e-6 + 18.0e-6 * np.cos(2.0 * np.pi * 200.0 * time) + data = np.column_stack([time, values]) + components = np.array([[148.0e-6, 0.0, 0.0], [18.0e-6, 200.0, 0.0]]) + output = tmp_path / "melt_pool_signal.png" + + returned_path = plot_melt_pool_signal(data, components, 148.0e-6, 18.0e-6, output) + + assert returned_path == output + image = mpimg.imread(output) + assert image.shape[:2] == (900, 1950)