Add structural assertions for plot tests

python/tests/unit/aout/test_analyze_error_by_phase.py

@@ -23,6 +23,7 @@ def test_analyze_error_by_phase_pure_cases():
     A = 0.49
     DC = 0.5
     phase_clean = 2 * np.pi * Fin * t
+    rng = np.random.default_rng(2026062210)
 
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}, A={A}")
 
@@ -38,9 +39,9 @@ def test_analyze_error_by_phase_pure_cases():
 
     # Generate signals
     for case in test_cases:
-        am_noise = np.random.randn(N) * case['am'] if case['am'] > 0 else 0
-        pm_noise = np.random.randn(N) * case['pm'] / A if case['pm'] > 0 else 0
-        th_noise = np.random.randn(N) * case['thermal'] if case['thermal'] > 0 else 0
+        am_noise = rng.standard_normal(N) * case['am'] if case['am'] > 0 else 0
+        pm_noise = rng.standard_normal(N) * case['pm'] / A if case['pm'] > 0 else 0
+        th_noise = rng.standard_normal(N) * case['thermal'] if case['thermal'] > 0 else 0
         case['signal'] = (A + am_noise) * np.sin(phase_clean + pm_noise) + DC + th_noise
 
     # Test both baseline modes
@@ -93,6 +94,7 @@ def test_analyze_error_by_phase_mixed_cases():
     A = 0.49
     DC = 0.5
     phase_clean = 2 * np.pi * Fin * t
+    rng = np.random.default_rng(2026062211)
 
     # Define mixed test cases
     test_cases = [
@@ -106,9 +108,9 @@ def test_analyze_error_by_phase_mixed_cases():
 
     # Generate signals
     for case in test_cases:
-        am_noise = np.random.randn(N) * case['am'] if case['am'] > 0 else 0
-        pm_noise = np.random.randn(N) * case['pm'] / A if case['pm'] > 0 else 0
-        th_noise = np.random.randn(N) * case['thermal'] if case['thermal'] > 0 else 0
+        am_noise = rng.standard_normal(N) * case['am'] if case['am'] > 0 else 0
+        pm_noise = rng.standard_normal(N) * case['pm'] / A if case['pm'] > 0 else 0
+        th_noise = rng.standard_normal(N) * case['thermal'] if case['thermal'] > 0 else 0
         case['signal'] = (A + am_noise) * np.sin(phase_clean + pm_noise) + DC + th_noise
 
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))

python/tests/unit/aout/test_analyze_error_by_value.py

@@ -22,15 +22,16 @@ def test_analyze_error_by_value_thermal_vs_nonlinearity():
     A = 0.49
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062212)
 
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = A * np.sin(2 * np.pi * Fin * t) + DC + np.random.randn(N) * base_noise
+    sig_noise = A * np.sin(2 * np.pi * Fin * t) + DC + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with 3rd Order Nonlinearity
     k3 = 0.01
-    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create figure with 3 subplots
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
@@ -65,8 +66,9 @@ def test_analyze_error_by_value_bin_count(n_bins):
     DC = 0.5
     base_noise = 50e-6
     k3 = 0.01
+    rng = np.random.default_rng(2026062213)
 
-    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create a simple plot
     fig, ax = plt.subplots(1, 1, figsize=(6, 5))
@@ -86,17 +88,18 @@ def test_analyze_error_by_value_different_nonlinearities():
     A = 0.49
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062214)
 
     # Case 1: Thermal only
-    sig_thermal = A * np.sin(2 * np.pi * Fin * t) + DC + np.random.randn(N) * base_noise
+    sig_thermal = A * np.sin(2 * np.pi * Fin * t) + DC + rng.standard_normal(N) * base_noise
 
     # Case 2: 2nd order nonlinearity
     k2 = 0.01
-    sig_k2 = A * np.sin(2 * np.pi * Fin * t) + DC + k2 * (A * np.sin(2 * np.pi * Fin * t))**2 + np.random.randn(N) * base_noise
+    sig_k2 = A * np.sin(2 * np.pi * Fin * t) + DC + k2 * (A * np.sin(2 * np.pi * Fin * t))**2 + rng.standard_normal(N) * base_noise
 
     # Case 3: 3rd order nonlinearity
     k3 = 0.01
-    sig_k3 = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_k3 = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create figure with 3 subplots
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))

python/tests/unit/aout/test_decompose_harmonics.py

@@ -22,25 +22,26 @@ def test_decompose_harmonics_basic():
     A = 0.25
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062215)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
     sig_ideal = sig_ac + DC
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = sig_ideal + np.random.randn(N) * base_noise
+    sig_noise = sig_ideal + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with Nonlinearity
     k2 = 0.001
     k3 = 0.005
-    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Case 3: ADC with Glitches
     glitch_prob = 0.01
     glitch_amplitude = 0.1
-    glitch_mask = np.random.rand(N) < glitch_prob
+    glitch_mask = rng.random(N) < glitch_prob
     glitch = glitch_mask * glitch_amplitude
-    sig_glitch = sig_ideal + glitch + np.random.randn(N) * base_noise
+    sig_glitch = sig_ideal + glitch + rng.standard_normal(N) * base_noise
 
     # Create figure
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
@@ -74,9 +75,10 @@ def test_decompose_harmonics_nonlinearity_levels(k3_value):
     A = 0.25
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062216)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
-    sig_nonlin = DC + sig_ac + k3_value * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k3_value * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Just run the analysis without creating plot
     fig, ax = plt.subplots(1, 1, figsize=(6, 5))

python/tests/unit/aout/test_decompose_harmonics_polar.py

@@ -24,25 +24,26 @@ def test_decompose_harmonics_polar_basic():
     DC = 0.5
     base_noise = 50e-6
     adc_range = [0, 1]
+    rng = np.random.default_rng(2026062217)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
     sig_ideal = sig_ac + DC
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.6f} MHz, Bin={Fin_bin}, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = sig_ideal + np.random.randn(N) * base_noise
+    sig_noise = sig_ideal + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with Nonlinearity
     k2 = 0.001
     k3 = 0.005
-    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Case 3: ADC with Glitches
     glitch_prob = 0.01
     glitch_amplitude = 0.1
-    glitch_mask = np.random.rand(N) < glitch_prob
+    glitch_mask = rng.random(N) < glitch_prob
     glitch = glitch_mask * glitch_amplitude
-    sig_glitch = sig_ideal + glitch + np.random.randn(N) * base_noise
+    sig_glitch = sig_ideal + glitch + rng.standard_normal(N) * base_noise
 
     # Create 2x3 subplot grid
     fig = plt.figure(figsize=(18, 12))

python/tests/unit/aout/test_fit_sine_4param.py

@@ -21,11 +21,12 @@ def test_fit_sine_4param_basic():
     A = 0.499
     DC = 0.5
     noise_rms = 20e-3
+    rng = np.random.default_rng(2026062218)
 
     print(f"\n[Sinewave] [Fs={Fs/1e6:.1f} MHz] [Fin={Fin/1e6:.6f} MHz] [Amplitude={A:.3f} V] [DC={DC:.3f} V] [Noise RMS={noise_rms*1e3:.2f} mV]")
 
     sig_ideal = A * np.cos(2 * np.pi * Fin * t) + DC
-    sig_noisy = sig_ideal + np.random.randn(N) * noise_rms
+    sig_noisy = sig_ideal + rng.standard_normal(N) * noise_rms
 
     # Fit 4-parameter sine wave
     result = fit_sine_4param(sig_noisy)
@@ -110,9 +111,10 @@ def test_fit_sine_4param_noise_levels(noise_level):
     t = np.arange(N) / Fs
     A = 0.499
     DC = 0.5
+    rng = np.random.default_rng(2026062219)
 
     sig_ideal = A * np.cos(2 * np.pi * Fin * t) + DC
-    sig_noisy = sig_ideal + np.random.randn(N) * noise_level
+    sig_noisy = sig_ideal + rng.standard_normal(N) * noise_level
 
     result = fit_sine_4param(sig_noisy)
     sig_fit = result['fitted_signal']

python/tests/unit/aout/test_inl_from_sine_sweep_length.py

@@ -25,6 +25,7 @@ def test_inl_from_sine_sweep_length():
     DC = 0.5
     base_noise = 50e-6
     hd2_dB, hd3_dB = -80, -66
+    rng = np.random.default_rng(2026062220)
 
     # Compute HD coefficients
     hd2_amp = 10**(hd2_dB/20)
@@ -43,7 +44,7 @@ def test_inl_from_sine_sweep_length():
         fin, J = find_coherent_frequency(fs, fin_target, N)
         t = np.arange(N) / fs
         sinewave = A * np.sin(2 * np.pi * fin * t)
-        signal_distorted = sinewave + k2 * sinewave**2 + k3 * sinewave**3 + DC + np.random.randn(N) * base_noise
+        signal_distorted = sinewave + k2 * sinewave**2 + k3 * sinewave**3 + DC + rng.standard_normal(N) * base_noise
 
         result = analyze_spectrum(signal_distorted, fs=fs, create_plot=False)
 

python/tests/unit/aout/test_jitter_calculation.py

@@ -21,6 +21,7 @@ def test_jitter_recovery_at_1ghz():
     A = 0.49
     DC = 0.0
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062221)
 
     # Test 3 frequencies
     fin_targets = [100e6, 1000e6, 2000e6]  # 100 MHz, 1 GHz, 2 GHz
@@ -52,9 +53,9 @@ def test_jitter_recovery_at_1ghz():
 
                 # Phase jitter model
                 phase_noise_rms = 2 * np.pi * Fin * jitter_rms
-                phase_jitter = np.random.randn(N) * phase_noise_rms
+                phase_jitter = rng.standard_normal(N) * phase_noise_rms
 
-                signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC + np.random.randn(N) * noise_level
+                signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC + rng.standard_normal(N) * noise_level
 
                 # Measure jitter using analyze_error_by_phase
                 results = analyze_error_by_phase(signal, norm_freq=Fin/Fs, n_bins=100,
@@ -164,6 +165,7 @@ def test_jitter_recovery_frequency_sweep(fin_target, expected_correlation):
     Fs = 7e9
     A = 0.49
     DC = 0.0
+    rng = np.random.default_rng(2026062222)
 
     # Find coherent frequency
     Fin, Fin_bin = find_coherent_frequency(fs=Fs, fin_target=fin_target, n_fft=N)
@@ -179,7 +181,7 @@ def test_jitter_recovery_frequency_sweep(fin_target, expected_correlation):
         # Generate signal
         t = np.arange(N) / Fs
         phase_noise_rms = 2 * np.pi * Fin * jitter_rms
-        phase_jitter = np.random.randn(N) * phase_noise_rms
+        phase_jitter = rng.standard_normal(N) * phase_noise_rms
         signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC
 
         # Measure jitter

python/tests/unit/aout/test_verify_spec_plot_phase.py

@@ -21,8 +21,7 @@ def generate_test_signal():
         signal: Test signal with 4 harmonics + noise
         params: Dictionary of expected values
     """
-    # Set random seed for reproducibility
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
 
     N = 8192           # Number of samples
     Fs = 1e9           # Sampling frequency
@@ -49,7 +48,7 @@ def generate_test_signal():
               A_HD4 * np.sin(2*np.pi*4*Fin*t + phi_HD4))
 
     # Add DC offset and small noise
-    signal = signal + 0.5 + np.random.randn(N) * 1e-5
+    signal = signal + 0.5 + rng.standard_normal(N) * 1e-5
 
     # Calculate maxSignal for normalization (needed for expected values)
     maxSignal = np.max(signal) - np.min(signal)

python/tests/unit/calibration/test_verify_calibration_full.py

@@ -62,7 +62,7 @@ def test_calibration_single_dataset_shuffled():
             shift_amounts_order = np.arange(bit_width - 1, -1, -1)
 
 
-            # shuffled_indices = np.random.permutation(bit_width)
+            # shuffled_indices = np.random.default_rng(2026062200).permutation(bit_width)
             shuffled_indices = np.arange(bit_width)
             shuffled_weights = true_weights[shuffled_indices]
             current_shifts = shift_amounts_order[shuffled_indices]
@@ -170,7 +170,7 @@ def test_calibration_single_dataset_shuffled_search_freq():
             shift_amounts_order = np.arange(bit_width - 1, -1, -1)
 
 
-            # shuffled_indices = np.random.permutation(bit_width)
+            # shuffled_indices = np.random.default_rng(2026062200).permutation(bit_width)
             shuffled_indices = np.arange(bit_width)
             shuffled_weights = true_weights[shuffled_indices]
             current_shifts = shift_amounts_order[shuffled_indices]

python/tests/unit/calibration/test_verify_lstsq_solver.py

@@ -142,12 +142,13 @@ def test_weight_recovery_absolute():
     4. Force the Solver's 'b' vector to match our Target Signal.
     """
     n_samples = 1000
+    rng = np.random.default_rng(2026062205)
     # 1. Define Ground Truth
     true_weights = np.array([0.25, 0.5, 1.0])
     bit_width = len(true_weights)
 
     # 2. Generate random bit matrix (n_samples, 3)
-    bits = np.random.randint(0, 2, (n_samples, bit_width)).astype(float)
+    bits = rng.integers(0, 2, (n_samples, bit_width)).astype(float)
 
     # 3. Construct the Target Signal directly from bits and weights
     # This signal is exactly what the solver should aim to reconstruct
@@ -238,15 +239,16 @@ def test_solve_weights_shared_recovery():
     Test 6: Verify that weights are shared across datasets but DCs are independent.
     """
     n_samples = 500
+    rng = np.random.default_rng(2026062206)
     true_weights = np.array([10.0, 20.0]) # Shared
     dc1, dc2 = 5.0, -3.0 # Different DCs
 
     # Dataset 1: bits @ weights + dc1
-    bits1 = np.random.randint(0, 2, (n_samples, 2)).astype(float)
+    bits1 = rng.integers(0, 2, (n_samples, 2)).astype(float)
     sig1 = bits1 @ true_weights + dc1
 
     # Dataset 2: bits @ weights + dc2
-    bits2 = np.random.randint(0, 2, (n_samples, 2)).astype(float)
+    bits2 = rng.integers(0, 2, (n_samples, 2)).astype(float)
     sig2 = bits2 @ true_weights + dc2
 
     # Mocking the basis to be the signals themselves to force direct recovery

python/tests/unit/calibration/test_verify_patch_rank_deficiency.py

@@ -10,7 +10,8 @@ def test_patch_rank_deficiency_logic():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062207)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -54,7 +55,8 @@ def test_patch_rank_deficiency_logic_reverse():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1][::-1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062208)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -99,7 +101,8 @@ def test_patch_rank_deficiency_logic_recover():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1][::-1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062209)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -124,4 +127,4 @@ def test_patch_rank_deficiency_logic_recover():
 
     # --- Assertions ---
     assert weights_recovered[1] == 0.0
-    assert weights_recovered[2] == 0.0
+    assert weights_recovered[2] == 0.0