diff --git a/examples/synthetic_api_example/run_synthetic_rve.py b/examples/synthetic_api_example/run_synthetic_rve.py index 086c18b..035e39a 100644 --- a/examples/synthetic_api_example/run_synthetic_rve.py +++ b/examples/synthetic_api_example/run_synthetic_rve.py @@ -40,6 +40,7 @@ MELT_POOL_HEIGHT = 30.0e-6 MELT_POOL_WIDTH_STD_DEV = 18.0e-6 MELT_POOL_LENGTH = 300.0e-6 +N_MODES = 50 HEIGHT_SHAPE_FACTOR = 1.0 DEPTH_SHAPE_FACTOR = 1.0 @@ -99,16 +100,16 @@ def build_melt_pool(): # assign shape to melt pool and cap (1 = parabola, 2 = ellipse) melt_pool_dict = { - "width": (width_data, None, 1.0, 2.0), + "width": (width_data, N_MODES, 1.0, 2.0), "depth": ( width_data, - None, + N_MODES, depth_scale, DEPTH_SHAPE_FACTOR, ), "height": ( width_data, - None, + N_MODES, height_scale, HEIGHT_SHAPE_FACTOR, ), diff --git a/src/raptor/api.py b/src/raptor/api.py index c6d2ead..bfc365f 100644 --- a/src/raptor/api.py +++ b/src/raptor/api.py @@ -9,6 +9,7 @@ # https://github.com/ORNL-MDF/Raptor/LICENSE # ============================================================================= import time +import warnings from typing import Any, Dict, List, Optional import numpy as np import pandas as pd @@ -77,18 +78,24 @@ def compute_spectral_components( ) -> 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. + ``n_modes`` includes the mean (DC) mode. By itself, it 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. If both are provided, ``n_modes`` caps + that set while preserving its highest-energy bins; the requested tolerance + cannot be met when the cap is smaller than the required set. """ 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.") + if n_modes is None and tolerance is None: + raise ValueError("Provide either n_modes, tolerance, or both.") + if n_modes is not None and ( + not isinstance(n_modes, (int, np.integer)) or n_modes < 1 + ): + raise ValueError("n_modes must be a positive integer.") time_values = data[:, 0] signal = data[:, 1] @@ -101,6 +108,10 @@ def compute_spectral_components( fft_values = np.fft.rfft(signal) frequencies = np.fft.rfftfreq(n_samples, d=dt) candidate_bins = np.arange(1, fft_values.size) + if n_modes is not None and n_modes > fft_values.size: + raise ValueError( + f"n_modes cannot exceed {fft_values.size} for this time series." + ) # Parseval energy represented by each positive-frequency cosine. Interior # rFFT bins represent a conjugate pair; the Nyquist bin does not. @@ -128,17 +139,21 @@ def compute_spectral_components( + 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." + if n_modes is not None and selected_bins.size > n_modes - 1: + warnings.warn( + f"Requested tolerance requires {selected_bins.size + 1} modes; " + f"limiting the result to {n_modes} modes, so the tolerance " + "will not be met.", + RuntimeWarning, + stacklevel=2, ) + selected_bins = selected_bins[: n_modes - 1] + else: selected_bins = candidate_bins[: n_modes - 1] - # Frequency order is convenient for evaluation and deterministic output. + # Return modes in frequency order after any energy-ranked selection and cap. 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]) @@ -172,7 +187,7 @@ def create_melt_pool( spectral_array = compute_spectral_components( data, n_modes=n_modes, - tolerance=component_tolerance if n_modes is None else None, + tolerance=component_tolerance, ) spectral_array[:, 0] *= scale diff --git a/tests/test_api.py b/tests/test_api.py index 4f8aed0..8d29432 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -379,6 +379,39 @@ def test_tolerance_selects_minimum_modes_and_meets_rmse(self): assert np.sqrt(np.mean((values - reconstructed) ** 2)) <= 0.03 assert compute_required_modes(data, 0.03) == 2 + def test_mode_cap_retains_highest_energy_tolerance_modes(self): + sampling_frequency = 1000.0 + time = np.arange(1000) / sampling_frequency + values = ( + 1.0 + + 0.1 * np.cos(2.0 * np.pi * 20.0 * time) + + 0.9 * np.cos(2.0 * np.pi * 170.0 * time) + ) + data = np.column_stack([time, values]) + + with pytest.warns(RuntimeWarning, match="tolerance will not be met"): + spectral_array = compute_spectral_components( + data, n_modes=2, tolerance=0.01 + ) + + assert spectral_array.shape == (2, 3) + assert spectral_array[1, 1] == pytest.approx(170.0) + + def test_tolerance_modes_are_returned_in_frequency_order(self): + sampling_frequency = 1000.0 + time = np.arange(1000) / sampling_frequency + values = ( + 1.0 + + 0.1 * np.cos(2.0 * np.pi * 20.0 * time) + + 0.9 * np.cos(2.0 * np.pi * 170.0 * time) + ) + + spectral_array = compute_spectral_components( + np.column_stack([time, values]), tolerance=0.01 + ) + + assert spectral_array[:, 1].tolist() == pytest.approx([0.0, 20.0, 170.0]) + 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) @@ -398,9 +431,26 @@ def test_compute_spectral_components_invalid_input(self): compute_spectral_components(np.array([[0.0, 1.0]]), 1) with pytest.raises(ValueError, match="shape"): compute_spectral_components(np.ones((4, 1)), 2) - with pytest.raises(ValueError, match="exactly one"): + with pytest.raises(ValueError, match="either n_modes"): compute_spectral_components(np.ones((4, 2))) + @pytest.mark.parametrize("n_modes", [0, -1, 1.5]) + def test_compute_spectral_components_rejects_invalid_mode_counts( + self, sample_time_series_data, n_modes + ): + with pytest.raises(ValueError, match="positive integer"): + compute_spectral_components(sample_time_series_data, n_modes=n_modes) + + def test_tolerance_rejects_mode_count_above_available_bins( + self, sample_time_series_data + ): + with pytest.raises(ValueError, match="cannot exceed"): + compute_spectral_components( + sample_time_series_data, + n_modes=sample_time_series_data.shape[0], + tolerance=0.0, + ) + # ============================================================================= # Tests for create_melt_pool