From 90d3e8a43679e3ed99ea9b9974166f391a3e491b Mon Sep 17 00:00:00 2001 From: John Date: Thu, 16 Jul 2026 13:06:17 -0400 Subject: [PATCH 1/7] Fixup for time scale in synthetic melt pool --- .../run_synthetic_rve.py | 160 +++++++++++++----- src/raptor/api.py | 52 ++++-- src/raptor/utilities.py | 22 ++- tests/test_api.py | 29 ++++ tests/test_utilities.py | 52 ++++++ 5 files changed, 254 insertions(+), 61 deletions(-) create mode 100644 tests/test_utilities.py diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index 78ce13a..6f9dfcd 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -9,6 +9,7 @@ # https://github.com/ORNL-MDF/Raptor/LICENSE # ============================================================================= import numpy as np +import matplotlib.pyplot as plt from raptor.api import ( compute_morphology, @@ -16,85 +17,159 @@ create_grid, create_melt_pool, create_path_vectors, + visualize, write_morphology, write_vtk, ) from raptor.utilities import MeltPoolFilter -# ============================================================================= -# 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. +N_SPECTRAL_MODES = 50 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]) +RVE_MAX_POINT = np.array([10.0e-4, 10.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" -def build_melt_pool(): - """Create a stochastic melt pool from the top-level configuration.""" - np.random.seed(RANDOM_SEED) +def plot_width_data(width_data): + time_ms = width_data[:, 0] * 1.0e3 + width_um = width_data[:, 1] * 1.0e6 + mean_um = MELT_POOL_WIDTH * 1.0e6 + std_dev_um = MELT_POOL_WIDTH_STD_DEV * 1.0e6 - # Internal signal-generation settings. These control the numerical - # representation of the stochastic melt-pool history. - sampling_frequency = 250000.0 - time_series_duration = 0.08 + plot_limits = ( + min(width_um.min(), mean_um - 4.0 * std_dev_um), + max(width_um.max(), mean_um + 4.0 * std_dev_um), + ) + gaussian_width = np.linspace(*plot_limits, 500) + gaussian_density = np.exp(-0.5 * ((gaussian_width - mean_um) / std_dev_um) ** 2) / ( + std_dev_um * np.sqrt(2.0 * np.pi) + ) - melt_pool_filter = MeltPoolFilter( + style = { + "font.size": 8, + "axes.labelsize": 8, + "legend.fontsize": 7, + "xtick.labelsize": 7, + "ytick.labelsize": 7, + } + with plt.rc_context(style): + fig, (time_axis, distribution_axis) = plt.subplots( + 1, + 2, + figsize=(6.5, 3), + constrained_layout=True, + ) + + time_axis.plot(time_ms, width_um, color="#0072B2", linewidth=0.8) + time_axis.axhline( + mean_um, + color="black", + linestyle="--", + linewidth=1.0, + label="Specified mean", + ) + time_axis.set_xlabel("Time (ms)") + time_axis.set_ylabel("Melt-pool width (µm)") + time_axis.legend(frameon=False, loc="upper right") + + distribution_axis.hist( + width_um, + bins=40, + density=True, + color="#56B4E9", + edgecolor="white", + linewidth=0.4, + alpha=0.75, + label="Filtered samples", + ) + distribution_axis.plot( + gaussian_width, + gaussian_density, + color="#D55E00", + linewidth=2.0, + label=f"Gaussian (µ={mean_um:.0f}, σ={std_dev_um:.0f} µm)", + ) + distribution_axis.set_xlim(plot_limits) + distribution_axis.set_xlabel("Melt-pool width (µm)") + distribution_axis.set_ylabel("Probability density (µm⁻¹)") + distribution_axis.legend(frameon=False, loc="upper right") + + for label, axis in zip(("(a)", "(b)"), (time_axis, distribution_axis)): + axis.text( + 0.02, + 0.96, + label, + transform=axis.transAxes, + fontweight="bold", + va="top", + ) + axis.spines["top"].set_visible(False) + axis.spines["right"].set_visible(False) + axis.tick_params(direction="out") + + fig.savefig(WIDTH_DATA_PLOT, dpi=300, facecolor="white") + plt.close(fig) + + +def build_melt_pool(): + # Create melt pools from convolution filter + fluctuation_time_scale = MELT_POOL_LENGTH / SCAN_SPEED + sampling_frequency = SCAN_SPEED / VOXEL_RESOLUTION + duration = 200.0 * fluctuation_time_scale + + # Instantiate object + mp_filter = MeltPoolFilter( MELT_POOL_WIDTH, MELT_POOL_WIDTH_STD_DEV, SCAN_SPEED, - [sampling_frequency, time_series_duration], + [sampling_frequency, duration], ) - # 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) + # Define physical scales + mp_filter.add_effect("melt_pool", [MELT_POOL_LENGTH, None, 1]) - # Depth and height reuse the width history, scaled to their requested means. + # Generate stochastic melt pool + mp_filter.initialize() + width_data = mp_filter.generate_fluctuations(1) + plot_width_data(width_data) + + # scale melt pool data by constant factor + depth_scale = MELT_POOL_DEPTH / MELT_POOL_WIDTH + height_scale = MELT_POOL_HEIGHT / MELT_POOL_WIDTH + + # 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, N_SPECTRAL_MODES, 1.0, 2.0), "depth": ( width_data, N_SPECTRAL_MODES, - MELT_POOL_DEPTH / MELT_POOL_WIDTH, + depth_scale, DEPTH_SHAPE_FACTOR, ), "height": ( width_data, N_SPECTRAL_MODES, - MELT_POOL_HEIGHT / MELT_POOL_WIDTH, + height_scale, HEIGHT_SHAPE_FACTOR, ), } @@ -103,12 +178,11 @@ def build_melt_pool(): 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, @@ -120,11 +194,16 @@ def main(): extra_layers=0, ) + # 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 +211,9 @@ def main(): ) write_morphology(morphology, MORPHOLOGY_OUTPUT) + # 7. Visualize using PyVista + visualize(VTK_OUTPUT) + if __name__ == "__main__": main() diff --git a/src/raptor/api.py b/src/raptor/api.py index 2fb6518..0218f8b 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -68,28 +68,54 @@ def create_path_vectors( def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: int) -> np.ndarray: 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) + values = melt_pool_data[:, 1] + mode0 = values.mean() + n_fft = len(values) 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] + fft_resolution = np.fft.rfft(values - mode0) + frequencies = np.fft.rfftfreq(n_fft, dt) + available_indices = np.arange(1, len(fft_resolution)) + n_fluctuation_modes = min(n_modes - 1, len(available_indices)) + + dominant = available_indices[ + np.argsort(np.abs(fft_resolution[available_indices]))[ + -n_fluctuation_modes: + ] + ] + dominant.sort() + + selected_fft = fft_resolution[dominant] + selected_frequencies = frequencies[dominant] + amplitudes = 2.0 * np.abs(selected_fft) / n_fft + + # The Nyquist term has no negative-frequency partner. + if n_fft % 2 == 0: + amplitudes[dominant == n_fft // 2] *= 0.5 + + # Preserve the variance of the input after truncating the spectrum. + represented_variance = 0.5 * np.sum(amplitudes**2) + if n_fft % 2 == 0: + represented_variance += 0.5 * np.sum( + amplitudes[dominant == n_fft // 2] ** 2 + ) + target_std = np.std(values) + if represented_variance > 0.0: + amplitudes *= target_std / np.sqrt(represented_variance) - frequencies = np.float64(1 / (dt * n_fft)) * np.arange( - n_modes, dtype=np.float64 - ) - phases = np.float64(np.angle(F[:n_modes])) - amplitudes = np.float64(np.abs(F[:n_modes]) / n_fft) + phases = np.angle(selected_fft) spectral_array = np.vstack( [ np.array([mode0, 0, 0]), - np.vstack([amplitudes[1:], frequencies[1:], phases[1:]]).transpose(), + np.vstack([amplitudes, selected_frequencies, phases]).transpose(), ] ) + + if spectral_array.shape[0] < n_modes: + spectral_array = np.vstack( + [spectral_array, np.zeros((n_modes - spectral_array.shape[0], 3))] + ) return np.float64(spectral_array) diff --git a/src/raptor/utilities.py b/src/raptor/utilities.py index 345efa8..d61d694 100644 --- a/src/raptor/utilities.py +++ b/src/raptor/utilities.py @@ -11,7 +11,7 @@ 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 class ScanPathBuilder: @@ -239,8 +239,8 @@ def __init__( # 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) + self.n_points = int(np.floor(self.duration * self.fs)) + 1 + self.t = np.arange(self.n_points) * self.dt # parametric representations of fluctuation scales self.physical_effects = {} # contains scale description and parameters @@ -287,12 +287,13 @@ def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4): 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) + 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): base_white_noise = np.random.normal( @@ -310,7 +311,10 @@ 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 ) diff --git a/tests/test_api.py b/tests/test_api.py index df22e9c..06ff1ff 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -332,6 +332,35 @@ 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 that the truncated cosine representation preserves variance.""" + spectral_array = compute_spectral_components(sample_time_series_data, 3) + 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_compute_spectral_components_selects_dominant_modes(self): + """Test that high-frequency signal content is not discarded.""" + 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(2000.0) + def test_compute_spectral_components_invalid_input(self): """Test spectral component computation with invalid input.""" with pytest.raises(IndexError): diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 0000000..3fadec7 --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,52 @@ +# ============================================================================= +# 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.utilities import MeltPoolFilter + + +def test_melt_pool_filter_recovers_requested_statistics(): + mean = 148.0e-6 + standard_deviation = 40.0e-6 + scan_speed = 1.7 + melt_pool_length = 300.0e-6 + time_scale = melt_pool_length / scan_speed + sampling_frequency = 340_000.0 + + melt_pool_filter = MeltPoolFilter( + mean, + standard_deviation, + scan_speed, + [sampling_frequency, 30.0 * time_scale], + ) + melt_pool_filter.add_effect("melt_pool", [melt_pool_length, None, 1.0]) + melt_pool_filter.initialize() + width_data = melt_pool_filter.generate_fluctuations(1.0) + + assert np.isfinite(width_data).all() + np.testing.assert_allclose(width_data[:, 1].mean(), mean, atol=1.0e-15) + np.testing.assert_allclose( + width_data[:, 1].std(), standard_deviation, atol=1.0e-15 + ) + + +def test_melt_pool_filter_remains_finite_at_high_sampling_frequency(): + melt_pool_filter = MeltPoolFilter( + 148.0e-6, + 40.0e-6, + 1.7, + [2_500_000.0, 0.002], + ) + melt_pool_filter.add_effect("melt_pool", [300.0e-6, None, 1.0]) + melt_pool_filter.initialize() + width_data = melt_pool_filter.generate_fluctuations(1.0) + + assert np.isfinite(width_data).all() From 562e0692d15417a01d2ebb37340ec122b926f08c Mon Sep 17 00:00:00 2001 From: John Date: Thu, 16 Jul 2026 13:06:36 -0400 Subject: [PATCH 2/7] Fixup for time scale in synthetic melt pool --- src/raptor/api.py | 4 +--- tests/test_api.py | 8 ++------ tests/test_utilities.py | 4 +--- 3 files changed, 4 insertions(+), 12 deletions(-) diff --git a/src/raptor/api.py b/src/raptor/api.py index 0218f8b..12ae4f9 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -80,9 +80,7 @@ def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: int) -> np. n_fluctuation_modes = min(n_modes - 1, len(available_indices)) dominant = available_indices[ - np.argsort(np.abs(fft_resolution[available_indices]))[ - -n_fluctuation_modes: - ] + np.argsort(np.abs(fft_resolution[available_indices]))[-n_fluctuation_modes:] ] dominant.sort() diff --git a/tests/test_api.py b/tests/test_api.py index 06ff1ff..1ff5a3a 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -346,18 +346,14 @@ def test_compute_spectral_components_preserves_standard_deviation( assert reconstructed.mean() == pytest.approx( sample_time_series_data[:, 1].mean() ) - assert reconstructed.std() == pytest.approx( - sample_time_series_data[:, 1].std() - ) + assert reconstructed.std() == pytest.approx(sample_time_series_data[:, 1].std()) def test_compute_spectral_components_selects_dominant_modes(self): """Test that high-frequency signal content is not discarded.""" 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 - ) + spectral_array = compute_spectral_components(np.column_stack([time, values]), 2) assert spectral_array[1, 1] == pytest.approx(2000.0) diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 3fadec7..f72baf8 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -33,9 +33,7 @@ def test_melt_pool_filter_recovers_requested_statistics(): assert np.isfinite(width_data).all() np.testing.assert_allclose(width_data[:, 1].mean(), mean, atol=1.0e-15) - np.testing.assert_allclose( - width_data[:, 1].std(), standard_deviation, atol=1.0e-15 - ) + np.testing.assert_allclose(width_data[:, 1].std(), standard_deviation, atol=1.0e-15) def test_melt_pool_filter_remains_finite_at_high_sampling_frequency(): From bd0f4f53ff2ffd3199fdbad10cf703aae97f76e1 Mon Sep 17 00:00:00 2001 From: Vamsi Subraveti Date: Fri, 17 Jul 2026 13:16:53 -0400 Subject: [PATCH 3/7] mpfilter fix for proper timeseries generation --- src/raptor/utilities.py | 86 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 77 insertions(+), 9 deletions(-) diff --git a/src/raptor/utilities.py b/src/raptor/utilities.py index d61d694..76a1173 100644 --- a/src/raptor/utilities.py +++ b/src/raptor/utilities.py @@ -12,6 +12,7 @@ from typing import List, Tuple from .structures import PathVector from scipy.signal import butter, sosfilt +from scipy.stats import chi2 class ScanPathBuilder: @@ -221,7 +222,7 @@ 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, confidence: float, ci_relative_width: float, voxel_resolution: float ): """ Filtration of disparate fluctuation scales to infer a melt pool oscillations sequence. @@ -237,10 +238,13 @@ def __init__( # 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(np.floor(self.duration * self.fs)) + 1 - self.t = np.arange(self.n_points) * self.dt + self.voxel_resolution = voxel_resolution + self.fs = self.scan_speed / self.voxel_resolution + + # confidence interval properties + self.confidence = confidence + self.ci_relative_width = ci_relative_width + # parametric representations of fluctuation scales self.physical_effects = {} # contains scale description and parameters @@ -280,6 +284,30 @@ def initialize(self): params["sigma_weight"] / self.normalization_factor ) * self.sigma + max_timescale = max( + p["length_scale_m"] / self.scan_speed + for p in self.physical_effects.values() + if p["length_scale_m"] is not None + ) + duration_guess = max_timescale + npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 + t_guess = np.arange(npoints_guess) * (1 / self.fs) + timeseries_guess = self.generate_fluctuations(self.sigma, npoints_guess, t_guess) + ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) + while True: + if ci_result["target_within_ci"] and ci_result["precision_satisfied"]: + break + duration_guess *=2 + npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 + t_guess = np.arange(npoints_guess) * (1 / self.fs) + timeseries_guess = self.generate_fluctuations(self.sigma, npoints_guess, t_guess) + ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) + print("Confidence interval evaluation: ", ci_result) + + self.duration = duration_guess + self.n_points = int(np.floor(self.duration * self.fs)) + 1 + self.t = t_guess + 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) @@ -295,11 +323,11 @@ def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4): sos = butter(order, [low, high], btype="band", output="sos") return sosfilt(sos, data) - def generate_fluctuations(self, noise_scale): + def generate_fluctuations(self, noise_scale, n_points, t): base_white_noise = np.random.normal( - loc=0, scale=noise_scale, size=self.n_points + loc=0, scale=noise_scale, size=n_points ) - final_series = np.zeros(self.n_points) + final_series = np.zeros(n_points) self.component_series = {} # Create each component series, scale it, and add to the final series @@ -324,4 +352,44 @@ def generate_fluctuations(self, noise_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))) # Ensure effective_n is at least 1 and not greater than n + + # chi-squared confidence interval for variance + dof = effective_n - 1 + 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 + + allowed_width = 2 * self.ci_relative_width * self.sigma**2 + target_within_ci = (lower_bound <= self.sigma**2 <= upper_bound) + precision_satisfied = (upper_bound - lower_bound) <= allowed_width + + 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, + } \ No newline at end of file From 1d45e52f129496ff50268039163c464073e50f33 Mon Sep 17 00:00:00 2001 From: Vamsi Subraveti Date: Fri, 17 Jul 2026 14:22:44 -0400 Subject: [PATCH 4/7] fixed spectral component determination with rmse construction --- .gitignore | 1 + .../run_synthetic_rve.py | 22 ++-- src/raptor/api.py | 107 ++++++++++-------- 3 files changed, 72 insertions(+), 58 deletions(-) diff --git a/.gitignore b/.gitignore index 63342b5..49720fa 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ build # other typical case file names *.vti *morphology.csv +*.png # vscode settings cache .vscode diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index 6f9dfcd..36e39c7 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -37,6 +37,9 @@ MELT_POOL_WIDTH_STD_DEV = 18.0e-6 MELT_POOL_LENGTH = 300.0e-6 +VARIANCE_SIGNIFICANCE = 0.025 # \alpha value for significance of variance test +CI_TOLERANCE_WINDOW = 0.01 # relative tolerance window around the user supplied std + N_SPECTRAL_MODES = 50 HEIGHT_SHAPE_FACTOR = 1.0 DEPTH_SHAPE_FACTOR = 1.0 @@ -133,16 +136,15 @@ def plot_width_data(width_data): def build_melt_pool(): # Create melt pools from convolution filter - fluctuation_time_scale = MELT_POOL_LENGTH / SCAN_SPEED - sampling_frequency = SCAN_SPEED / VOXEL_RESOLUTION - duration = 200.0 * fluctuation_time_scale - + # Instantiate object mp_filter = MeltPoolFilter( MELT_POOL_WIDTH, MELT_POOL_WIDTH_STD_DEV, SCAN_SPEED, - [sampling_frequency, duration], + VARIANCE_SIGNIFICANCE, + CI_TOLERANCE_WINDOW, + VOXEL_RESOLUTION ) # Define physical scales @@ -150,7 +152,7 @@ def build_melt_pool(): # Generate stochastic melt pool mp_filter.initialize() - width_data = mp_filter.generate_fluctuations(1) + width_data = mp_filter.generate_fluctuations(1, mp_filter.n_points, mp_filter.t) plot_width_data(width_data) # scale melt pool data by constant factor @@ -159,22 +161,22 @@ def build_melt_pool(): # assign shape to melt pool and cap (1 = parabola, 2 = ellipse) melt_pool_dict = { - "width": (width_data, N_SPECTRAL_MODES, 1.0, 2.0), + "width": (width_data, None, 1.0, 2.0), "depth": ( width_data, - N_SPECTRAL_MODES, + None, depth_scale, DEPTH_SHAPE_FACTOR, ), "height": ( width_data, - N_SPECTRAL_MODES, + 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(): diff --git a/src/raptor/api.py b/src/raptor/api.py index 12ae4f9..9c77db7 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -65,74 +65,85 @@ def create_path_vectors( return scan_path_builder.process_vectors() -def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: int) -> np.ndarray: - - dt = melt_pool_data[1, 0] - melt_pool_data[0, 0] - values = melt_pool_data[:, 1] - mode0 = values.mean() - n_fft = len(values) - if n_modes == 1: - spectral_array = np.array([[mode0, 0, 0]]) - else: - fft_resolution = np.fft.rfft(values - mode0) - frequencies = np.fft.rfftfreq(n_fft, dt) - available_indices = np.arange(1, len(fft_resolution)) - n_fluctuation_modes = min(n_modes - 1, len(available_indices)) - - dominant = available_indices[ - np.argsort(np.abs(fft_resolution[available_indices]))[-n_fluctuation_modes:] - ] - dominant.sort() +def compute_required_modes(data: np.ndarray, reconstruction_rmse: float) -> int: + """ + Computes the number of modes required to reconstruct the melt pool data within a specified RMSE threshold. + """ + fft_resolution = np.fft.fft(data[:, 1]) + F = np.zeros_like(fft_resolution) + n_fft = len(fft_resolution) + - selected_fft = fft_resolution[dominant] - selected_frequencies = frequencies[dominant] - amplitudes = 2.0 * np.abs(selected_fft) / n_fft +def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: Optional[int]=None, tolerance: Optional[float]=None) -> np.ndarray: - # The Nyquist term has no negative-frequency partner. - if n_fft % 2 == 0: - amplitudes[dominant == n_fft // 2] *= 0.5 + 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) - # Preserve the variance of the input after truncating the spectrum. - represented_variance = 0.5 * np.sum(amplitudes**2) - if n_fft % 2 == 0: - represented_variance += 0.5 * np.sum( - amplitudes[dominant == n_fft // 2] ** 2 + if n_modes is None: + if tolerance is None: + raise ValueError( + "Either n_modes or tolerance must be provided for spectral component computation." ) - target_std = np.std(values) - if represented_variance > 0.0: - amplitudes *= target_std / np.sqrt(represented_variance) + for k in range(1, n_fft // 2 + 1): + # Keep only the n_modes with largest magnitude + fft_truncated = fft_resolution.copy() + indices = np.argsort(np.abs(fft_resolution)) + fft_truncated[indices[:-k]] = 0 + + # Reconstruct signal + reconstructed_signal = np.fft.ifft(fft_truncated).real + + # Compute RMSE + rmse = np.sqrt(np.mean((melt_pool_data[:, 1] - reconstructed_signal) ** 2)) + + # Check if threshold is met + if rmse <= tolerance: + n_modes = k + break + print(f"Computed {n_modes} modes to achieve reconstruction RMSE of {rmse:.6f} within tolerance {tolerance}.") + + elif n_modes == 1: + spectral_array = np.array([[mode0, 0, 0]]) - phases = np.angle(selected_fft) - spectral_array = np.vstack( - [ - np.array([mode0, 0, 0]), - np.vstack([amplitudes, selected_frequencies, phases]).transpose(), - ] - ) + for i in range(1, n_modes): + F[i] = fft_resolution[i] + F[n_fft - i] = fft_resolution[n_fft - i] - if spectral_array.shape[0] < n_modes: - spectral_array = np.vstack( - [spectral_array, np.zeros((n_modes - spectral_array.shape[0], 3))] - ) + frequencies = np.float64(1 / (dt * n_fft)) * np.arange( + n_modes, dtype=np.float64 + ) + 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(), + ] + ) return np.float64(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) + # for _, nmodes, _, _ in melt_pool_dict.values(): + # max_modes = max(max_modes, n_modes if nmodes is not None else 0) # 2. Process each component into its spectral format + max_modes = 0 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) + spectral_array = compute_spectral_components(data, n_modes,tolerance) + max_modes = max(max_modes, spectral_array.shape[0]) spectral_array[:, 0] *= scale # Option B: Input data is a spectral array [amplitude, frequency, phase] From 9bb32e200930fd407d26cc59d2718a3d743b9ed6 Mon Sep 17 00:00:00 2001 From: Vamsi Subraveti Date: Fri, 17 Jul 2026 14:24:42 -0400 Subject: [PATCH 5/7] precommit lol --- .../run_synthetic_rve.py | 12 ++--- src/raptor/api.py | 31 +++++++------ src/raptor/utilities.py | 44 ++++++++++++------- 3 files changed, 52 insertions(+), 35 deletions(-) diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index 36e39c7..2fdde1f 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -37,8 +37,8 @@ MELT_POOL_WIDTH_STD_DEV = 18.0e-6 MELT_POOL_LENGTH = 300.0e-6 -VARIANCE_SIGNIFICANCE = 0.025 # \alpha value for significance of variance test -CI_TOLERANCE_WINDOW = 0.01 # relative tolerance window around the user supplied std +VARIANCE_SIGNIFICANCE = 0.025 # \alpha value for significance of variance test +CI_TOLERANCE_WINDOW = 0.01 # relative tolerance window around the user supplied std N_SPECTRAL_MODES = 50 HEIGHT_SHAPE_FACTOR = 1.0 @@ -136,7 +136,7 @@ def plot_width_data(width_data): def build_melt_pool(): # Create melt pools from convolution filter - + # Instantiate object mp_filter = MeltPoolFilter( MELT_POOL_WIDTH, @@ -144,7 +144,7 @@ def build_melt_pool(): SCAN_SPEED, VARIANCE_SIGNIFICANCE, CI_TOLERANCE_WINDOW, - VOXEL_RESOLUTION + VOXEL_RESOLUTION, ) # Define physical scales @@ -176,7 +176,9 @@ def build_melt_pool(): ), } - return create_melt_pool(melt_pool_dict, enable_random_phases=True, tolerance = VOXEL_RESOLUTION) + return create_melt_pool( + melt_pool_dict, enable_random_phases=True, tolerance=VOXEL_RESOLUTION + ) def main(): diff --git a/src/raptor/api.py b/src/raptor/api.py index 9c77db7..fe6402c 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -72,9 +72,13 @@ def compute_required_modes(data: np.ndarray, reconstruction_rmse: float) -> int: fft_resolution = np.fft.fft(data[:, 1]) F = np.zeros_like(fft_resolution) n_fft = len(fft_resolution) - -def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: Optional[int]=None, tolerance: Optional[float]=None) -> np.ndarray: + +def compute_spectral_components( + melt_pool_data: np.ndarray, + n_modes: Optional[int] = None, + tolerance: Optional[float] = None, +) -> np.ndarray: dt = melt_pool_data[1, 0] - melt_pool_data[0, 0] mode0 = melt_pool_data[:, 1].mean() @@ -92,19 +96,21 @@ def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: Optional[in fft_truncated = fft_resolution.copy() indices = np.argsort(np.abs(fft_resolution)) fft_truncated[indices[:-k]] = 0 - + # Reconstruct signal reconstructed_signal = np.fft.ifft(fft_truncated).real - + # Compute RMSE rmse = np.sqrt(np.mean((melt_pool_data[:, 1] - reconstructed_signal) ** 2)) - + # Check if threshold is met if rmse <= tolerance: n_modes = k break - print(f"Computed {n_modes} modes to achieve reconstruction RMSE of {rmse:.6f} within tolerance {tolerance}.") - + print( + f"Computed {n_modes} modes to achieve reconstruction RMSE of {rmse:.6f} within tolerance {tolerance}." + ) + elif n_modes == 1: spectral_array = np.array([[mode0, 0, 0]]) @@ -112,9 +118,7 @@ def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: Optional[in 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 - ) + frequencies = np.float64(1 / (dt * n_fft)) * np.arange(n_modes, dtype=np.float64) phases = np.float64(np.angle(F[:n_modes])) amplitudes = np.float64(np.abs(F[:n_modes]) / n_fft) spectral_array = np.vstack( @@ -127,11 +131,12 @@ def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: Optional[in def create_melt_pool( - melt_pool_dict: Dict[str, Any], enable_random_phases: bool, tolerance: Optional[float] = None + melt_pool_dict: Dict[str, Any], + enable_random_phases: bool, + tolerance: Optional[float] = None, ) -> MeltPool: processed_components: Dict[str, Tuple[np.ndarray, float]] = {} - # 1. Determine the maximum number of modes required. # for _, nmodes, _, _ in melt_pool_dict.values(): @@ -142,7 +147,7 @@ def create_melt_pool( 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,tolerance) + spectral_array = compute_spectral_components(data, n_modes, tolerance) max_modes = max(max_modes, spectral_array.shape[0]) spectral_array[:, 0] *= scale diff --git a/src/raptor/utilities.py b/src/raptor/utilities.py index 76a1173..6ff18b9 100644 --- a/src/raptor/utilities.py +++ b/src/raptor/utilities.py @@ -222,7 +222,13 @@ def write_layers(self, output_name, mode="layers"): class MeltPoolFilter: def __init__( - self, mu: float, sigma: float, scan_speed: float, confidence: float, ci_relative_width: float, voxel_resolution: float + self, + mu: float, + sigma: float, + scan_speed: float, + confidence: float, + ci_relative_width: float, + voxel_resolution: float, ): """ Filtration of disparate fluctuation scales to infer a melt pool oscillations sequence. @@ -292,18 +298,22 @@ def initialize(self): duration_guess = max_timescale npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 t_guess = np.arange(npoints_guess) * (1 / self.fs) - timeseries_guess = self.generate_fluctuations(self.sigma, npoints_guess, t_guess) + timeseries_guess = self.generate_fluctuations( + self.sigma, npoints_guess, t_guess + ) ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) while True: if ci_result["target_within_ci"] and ci_result["precision_satisfied"]: break - duration_guess *=2 + duration_guess *= 2 npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 t_guess = np.arange(npoints_guess) * (1 / self.fs) - timeseries_guess = self.generate_fluctuations(self.sigma, npoints_guess, t_guess) + timeseries_guess = self.generate_fluctuations( + self.sigma, npoints_guess, t_guess + ) ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) print("Confidence interval evaluation: ", ci_result) - + self.duration = duration_guess self.n_points = int(np.floor(self.duration * self.fs)) + 1 self.t = t_guess @@ -324,9 +334,7 @@ def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4): return sosfilt(sos, data) def generate_fluctuations(self, noise_scale, n_points, t): - base_white_noise = np.random.normal( - loc=0, scale=noise_scale, size=n_points - ) + base_white_noise = np.random.normal(loc=0, scale=noise_scale, size=n_points) final_series = np.zeros(n_points) self.component_series = {} @@ -353,7 +361,7 @@ def generate_fluctuations(self, noise_scale, n_points, t): final_series += self.mu return np.column_stack([t, final_series]) - + def evaluate_variance_ci(self, data): """ Evaluates the confidence interval for the variance of the data. @@ -361,28 +369,30 @@ def evaluate_variance_ci(self, data): """ 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] + 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))) # Ensure effective_n is at least 1 and not greater than n + effective_n = int( + max(1, min(effective_n, n)) + ) # Ensure effective_n is at least 1 and not greater than n # chi-squared confidence interval for variance dof = effective_n - 1 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 - + allowed_width = 2 * self.ci_relative_width * self.sigma**2 - target_within_ci = (lower_bound <= self.sigma**2 <= upper_bound) + target_within_ci = lower_bound <= self.sigma**2 <= upper_bound precision_satisfied = (upper_bound - lower_bound) <= allowed_width return { @@ -392,4 +402,4 @@ def evaluate_variance_ci(self, data): "upper_bound": upper_bound, "target_within_ci": target_within_ci, "precision_satisfied": precision_satisfied, - } \ No newline at end of file + } From 828f065ffc6c3e66cddfed4106901b7719a3957c Mon Sep 17 00:00:00 2001 From: John Date: Fri, 17 Jul 2026 17:00:13 -0400 Subject: [PATCH 6/7] changes to autoselection of melt pool filter --- .gitignore | 3 + examples/.DS_Store | Bin 6148 -> 0 bytes .../run_synthetic_rve.py | 123 ++----- pyproject.toml | 2 + src/raptor/api.py | 158 +++++---- src/raptor/utilities.py | 329 +++++++++++++++--- tests/test_api.py | 74 +++- tests/test_utilities.py | 138 ++++++-- 8 files changed, 573 insertions(+), 254 deletions(-) delete mode 100644 examples/.DS_Store diff --git a/.gitignore b/.gitignore index 49720fa..8ad58ca 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,6 @@ build # vscode settings cache .vscode + +# macOS Finder metadata +.DS_Store diff --git a/examples/.DS_Store b/examples/.DS_Store deleted file mode 100644 index 55d96b26e7c7bf14ec34313ded502889cda7b2f7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHK%}T>S5T0$TO(;SS3VI88E!aOQ;w8lT0!H+pQWFw1m}X0VtWXNM>kIiJK94iI zTd~w0yoi(;nEhsFXR_?KVP_Wrgg1#A09636Pzeh*Ha`fBlded?cnF31Mhq_yKpzH> zu0*rrKQcgjSHV3SKmtAZzJGBz9Sy=%h8Xc44$@>auGc?Ap;+45E?Z@*V%_dN|}ZWI|#3mc-pJ(p2=(+Bw0LE2}u-V$nABKMRMxNNfwP%uCE6y+p>Gr zgW0Un?9^+{ajQMAIkS@{+U-_nKDVvC!=uxS-c#Do=>MjSyAIu%i;3iHGeIvxFviE|C+8g)7dH8Xz4%q+|cMX1@)@2GGPu10Q| z0cK#CfuinKss5j@fB!EhagP~b1~!TTQR@0#50_+W>&oJ&);g&7s3a7ZYy3( int: - """ - Computes the number of modes required to reconstruct the melt pool data within a specified RMSE threshold. - """ - fft_resolution = np.fft.fft(data[:, 1]) - F = np.zeros_like(fft_resolution) - n_fft = len(fft_resolution) + """Return the minimum number of real Fourier modes for an RMSE target.""" + return compute_spectral_components(data, tolerance=reconstruction_rmse).shape[0] def compute_spectral_components( @@ -79,55 +75,82 @@ def compute_spectral_components( n_modes: Optional[int] = None, tolerance: Optional[float] = None, ) -> np.ndarray: + """Convert a uniformly sampled real signal to a sparse cosine expansion. - 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 is None: - if tolerance is None: - raise ValueError( - "Either n_modes or tolerance must be provided for spectral component computation." + ``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 + ) + 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" ) - for k in range(1, n_fft // 2 + 1): - # Keep only the n_modes with largest magnitude - fft_truncated = fft_resolution.copy() - indices = np.argsort(np.abs(fft_resolution)) - fft_truncated[indices[:-k]] = 0 - - # Reconstruct signal - reconstructed_signal = np.fft.ifft(fft_truncated).real - - # Compute RMSE - rmse = np.sqrt(np.mean((melt_pool_data[:, 1] - reconstructed_signal) ** 2)) - - # Check if threshold is met - if rmse <= tolerance: - n_modes = k - break - print( - f"Computed {n_modes} modes to achieve reconstruction RMSE of {rmse:.6f} within tolerance {tolerance}." + + 1 ) + 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] - elif n_modes == 1: - spectral_array = np.array([[mode0, 0, 0]]) - - 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) - 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(), - ] - ) - return np.float64(spectral_array) + # 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( @@ -136,19 +159,21 @@ def create_melt_pool( tolerance: Optional[float] = None, ) -> MeltPool: - processed_components: Dict[str, Tuple[np.ndarray, float]] = {} + processed_components: Dict[str, np.ndarray] = {} - # 1. Determine the maximum number of modes required. - # for _, nmodes, _, _ in melt_pool_dict.values(): - # max_modes = max(max_modes, n_modes if nmodes is not None else 0) - - # 2. Process each component into its spectral format - max_modes = 0 + # 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, tolerance) - max_modes = max(max_modes, spectral_array.shape[0]) + 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] @@ -161,16 +186,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/utilities.py b/src/raptor/utilities.py index 6ff18b9..1e400e5 100644 --- a/src/raptor/utilities.py +++ b/src/raptor/utilities.py @@ -8,8 +8,11 @@ # 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 butter, sosfilt from scipy.stats import chi2 @@ -226,33 +229,30 @@ def __init__( mu: float, sigma: float, scan_speed: float, - confidence: float, - ci_relative_width: 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.voxel_resolution = voxel_resolution self.fs = self.scan_speed / self.voxel_resolution - - # confidence interval properties self.confidence = confidence self.ci_relative_width = ci_relative_width - - # parametric representations of fluctuation scales - self.physical_effects = {} # contains scale description and parameters + 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): """ @@ -267,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 @@ -285,43 +295,101 @@ 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( - p["length_scale_m"] / self.scan_speed - for p in self.physical_effects.values() - if p["length_scale_m"] is not None + 1.0 / p["frequency_hz"] for p in self.physical_effects.values() ) - duration_guess = max_timescale - npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 - t_guess = np.arange(npoints_guess) * (1 / self.fs) - timeseries_guess = self.generate_fluctuations( - self.sigma, npoints_guess, t_guess - ) - ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) - while True: - if ci_result["target_within_ci"] and ci_result["precision_satisfied"]: + 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 - duration_guess *= 2 - npoints_guess = int(np.floor(duration_guess * self.fs)) + 1 - t_guess = np.arange(npoints_guess) * (1 / self.fs) - timeseries_guess = self.generate_fluctuations( - self.sigma, npoints_guess, t_guess + upper_points = 2 * upper_points - 1 + else: + raise RuntimeError( + "Unable to satisfy the variance precision within " + "16 duration refinements." ) - ci_result = self.evaluate_variance_ci(timeseries_guess[:, 1]) - print("Confidence interval evaluation: ", ci_result) - self.duration = duration_guess - self.n_points = int(np.floor(self.duration * self.fs)) + 1 - self.t = t_guess + # 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 @@ -334,7 +402,7 @@ def bandpass_filter(self, data, f0, bandwidth_fraction, fs, order=4): return sosfilt(sos, data) def generate_fluctuations(self, noise_scale, n_points, t): - base_white_noise = np.random.normal(loc=0, scale=noise_scale, size=n_points) + base_white_noise = self.rng.normal(loc=0.0, scale=noise_scale, size=n_points) final_series = np.zeros(n_points) self.component_series = {} @@ -357,6 +425,17 @@ def generate_fluctuations(self, noise_scale, n_points, t): 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 @@ -378,12 +457,13 @@ def evaluate_variance_ci(self, data): 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)) - ) # Ensure effective_n is at least 1 and not greater than n + effective_n = int(max(1, min(effective_n, n))) + + return self._variance_ci(sample_variance, effective_n) - # chi-squared confidence interval for variance - dof = effective_n - 1 + 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) @@ -391,9 +471,13 @@ def evaluate_variance_ci(self, data): lower_bound = dof * sample_variance / chi2_upper upper_bound = dof * sample_variance / chi2_lower - allowed_width = 2 * self.ci_relative_width * self.sigma**2 target_within_ci = lower_bound <= self.sigma**2 <= upper_bound - precision_satisfied = (upper_bound - lower_bound) <= allowed_width + 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, @@ -403,3 +487,132 @@ def evaluate_variance_ci(self, data): "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 1ff5a3a..092c294 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]) @@ -335,8 +336,8 @@ def test_compute_spectral_components_mean_value(self, sample_time_series_data): def test_compute_spectral_components_preserves_standard_deviation( self, sample_time_series_data ): - """Test that the truncated cosine representation preserves variance.""" - spectral_array = compute_spectral_components(sample_time_series_data, 3) + """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) @@ -348,21 +349,57 @@ def test_compute_spectral_components_preserves_standard_deviation( ) assert reconstructed.std() == pytest.approx(sample_time_series_data[:, 1].std()) - def test_compute_spectral_components_selects_dominant_modes(self): - """Test that high-frequency signal content is not discarded.""" + 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(2000.0) + 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))) # ============================================================================= @@ -406,8 +443,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, @@ -415,7 +452,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.""" @@ -435,6 +472,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_utilities.py b/tests/test_utilities.py index f72baf8..36be2b2 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -8,43 +8,129 @@ # 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 +from raptor.utilities import ( + MeltPoolFilter, + plot_melt_pool_signal, + reconstruct_spectral_signal, +) -def test_melt_pool_filter_recovers_requested_statistics(): - mean = 148.0e-6 - standard_deviation = 40.0e-6 - scan_speed = 1.7 - melt_pool_length = 300.0e-6 - time_scale = melt_pool_length / scan_speed - sampling_frequency = 340_000.0 - +def make_filter(*, voxel_resolution=5.0e-6, random_seed=None): melt_pool_filter = MeltPoolFilter( - mean, - standard_deviation, - scan_speed, - [sampling_frequency, 30.0 * time_scale], + 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", [melt_pool_length, None, 1.0]) + 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) + 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(), mean, atol=1.0e-15) - np.testing.assert_allclose(width_data[:, 1].std(), standard_deviation, atol=1.0e-15) + 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_remains_finite_at_high_sampling_frequency(): - melt_pool_filter = MeltPoolFilter( - 148.0e-6, - 40.0e-6, - 1.7, - [2_500_000.0, 0.002], +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 ) - melt_pool_filter.add_effect("melt_pool", [300.0e-6, None, 1.0]) + + 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() - width_data = melt_pool_filter.generate_fluctuations(1.0) + model_correlation = melt_pool_filter._model_autocorrelation() - assert np.isfinite(width_data).all() + 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) From a8a036d241970e47abd6d9fa72b40897f7b41192 Mon Sep 17 00:00:00 2001 From: John Date: Tue, 21 Jul 2026 10:45:07 -0400 Subject: [PATCH 7/7] Fix CI, API, and CLI for melt pool filter --- .../run_synthetic_rve.py | 6 +++-- src/raptor/api.py | 1 + src/raptor/cli.py | 13 ++++++--- src/raptor/io.py | 5 +++- tests/test_api.py | 18 +++++++++++++ tests/test_io.py | 27 +++++++++++++++++++ 6 files changed, 63 insertions(+), 7 deletions(-) create mode 100644 tests/test_io.py diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index dca018e..086c18b 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -51,6 +51,7 @@ VTK_OUTPUT = "rve.vti" MORPHOLOGY_OUTPUT = "rve_morphology.csv" WIDTH_DATA_PLOT = "melt_pool_width_timeseries.png" +ENABLE_VISUALIZATION = False def build_melt_pool(): @@ -152,8 +153,9 @@ def main(): ) write_morphology(morphology, MORPHOLOGY_OUTPUT) - # 7. Visualize using PyVista - visualize(VTK_OUTPUT) + # 7. Optionally visualize using PyVista (requires a graphical display) + if ENABLE_VISUALIZATION: + visualize(VTK_OUTPUT) if __name__ == "__main__": diff --git a/src/raptor/api.py b/src/raptor/api.py index b49db92..c6d2ead 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -179,6 +179,7 @@ def create_melt_pool( # 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( 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/tests/test_api.py b/tests/test_api.py index 092c294..4f8aed0 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -435,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]]) 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)