From fe4582513a943e8339d9ecb230fa074ff9e05c2e Mon Sep 17 00:00:00 2001 From: Luisa Date: Mon, 4 May 2026 14:47:48 -0600 Subject: [PATCH 1/8] updates --- imap_processing/ultra/constants.py | 2 +- imap_processing/ultra/l1b/extendedspin.py | 21 ++++++++------- .../ultra/l1b/ultra_l1b_culling.py | 27 +++++++++++++++---- 3 files changed, 35 insertions(+), 15 deletions(-) diff --git a/imap_processing/ultra/constants.py b/imap_processing/ultra/constants.py index a4a0e0eb79..5112204267 100644 --- a/imap_processing/ultra/constants.py +++ b/imap_processing/ultra/constants.py @@ -206,7 +206,7 @@ class UltraConstants: MAX_ENERGY_THRESHOLD = 116.0 # Angle threshold in radians for ULTRA 45 degree culling. # This is only needed for ULTRA 45 since Earth may be in the FOV. - EARTH_ANGLE_45_THRESHOLD = np.radians(20) + EARTH_ANGLE_45_THRESHOLD = np.radians(15) # An array of energy thresholds to use for culling. Each one corresponds to # the number of energy bins used. # n_bins=len(PSET_ENERGY_BIN_EDGES)[BASE_CULL_EBIN:] // N_CULL_EBINS diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index 147b542afc..d32e84c143 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -75,6 +75,9 @@ def calculate_extendedspin( inst_qf = flag_imap_instruments(de_dataset["spin"].values) spin_bin_size = UltraConstants.SPIN_BIN_SIZE + print("SDC spin[:5]:", spin[:5]) + print("SDC spin_starttime[:5]:", spin_starttime[:5]) + print("SDC len(spin):", len(spin)) spin_tbin_edges = get_binned_spins_edges( spin, spin_period, spin_starttime, spin_bin_size ) @@ -91,6 +94,7 @@ def calculate_extendedspin( intervals, _, _ = build_energy_bins() # Get the energy ranges energy_ranges = get_binned_energy_ranges(intervals) + print("SDC energy_ranges:", energy_ranges) energy_bin_flags = get_energy_range_flags(energy_ranges) # Calculate the high energy quality flags energy_thresholds = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS @@ -102,9 +106,10 @@ def calculate_extendedspin( energy_thresholds, instrument_id, ) - # Combine high energy and voltage flags to use for statistical outlier flagging. - mask = ( - voltage_qf[np.newaxis, :] | high_energy_qf + # For the following culls, mask the spins that have already been flagged for + # low voltage + mask = np.repeat( + voltage_qf[np.newaxis, :], len(energy_ranges) - 1, axis=0 ) # Shape (n_energy_bins, n_spins_bins) upstream_ion_qf_1 = flag_upstream_ion( de_dataset, @@ -114,9 +119,6 @@ def calculate_extendedspin( UltraConstants.UPSTREAM_ION_ENERGY_CHANNELS_1, instrument_id, ) - # Update mask to include upstream ion flags from the first set of energy channels - # before flagging with the second set of energy channels - mask = mask | upstream_ion_qf_1 upstream_ion_qf_2 = flag_upstream_ion( de_dataset, spin_tbin_edges, @@ -129,12 +131,13 @@ def calculate_extendedspin( de_dataset, spin_tbin_edges, energy_ranges, + mask, UltraConstants.SPECTRAL_ENERGY_CHANNELS, instrument_id, ) - # Update mask to include upstream ion flags #2 and spectral flags before flagging - # statistical outliers - mask = mask | upstream_ion_qf_2 | spectral_qf + # Update mask to include high energy, upstream ion flags and spectral flags + # before flagging statistical outliers + mask = mask | upstream_ion_qf_1 | upstream_ion_qf_2 | spectral_qf | high_energy_qf stat_outliers_qf, _, _, _ = flag_statistical_outliers( de_dataset, spin_tbin_edges, diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index 2261cc9039..ab2f63f73c 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -1015,6 +1015,11 @@ def flag_upstream_ion( f" independent flags and will be applied across all {mask.shape[0]} energy" f" channels." ) + # SDC - add inside flag_upstream_ion + # SDC - inside flag_upstream_ion, first call only + print("SDC upstream_1 total_scaled[:5]:", total_scaled[:5]) + print("SDC upstream_1 thresh:", thresh) + print("SDC upstream_1 flagged bins:", np.where(flagged)[0]) return flagged @@ -1022,6 +1027,7 @@ def flag_spectral_events( de_dataset: xr.Dataset, spin_tbin_edges: NDArray, energy_ranges: NDArray, + mask: NDArray, channels: list, sensor_id: int = 90, ) -> NDArray: @@ -1036,6 +1042,8 @@ def flag_spectral_events( Edges of the spin time bins. energy_ranges : NDArray Array of energy range edges. + mask : NDArray + Mask indicating which events to consider for spectral flagging. channels : list List of energy channel indices to use for spectral flagging. sensor_id : int @@ -1058,6 +1066,9 @@ def flag_spectral_events( counts_sum = get_valid_de_count_summary( de_dataset, energy_ranges, spin_tbin_edges, sensor_id=sensor_id )[channels, :] # shape (num_channels, n_spin_bins) + valid_bins = np.all( + ~mask[channels, :], axis=0 + ) # Get valid spin bins across the selected channels; shape (n_spin_bins,) # Flag spin bins where the signed count difference between adjacent selected # energy channels exceeds a Poisson-based threshold. For each pair of # adjacent channels, compute np.diff(counts_sum, axis=0) and compare that @@ -1065,10 +1076,12 @@ def flag_spectral_events( # uncertainty (sqrt(N1 + N2)) of those two channels for each spin bin. # If any adjacent channel pair exceeds the threshold for a spin bin, that # spin bin is flagged across all energy ranges. - diff = np.diff(counts_sum, axis=0) - UltraConstants.SPECTRAL_SIG_THRESHOLD * ( - np.sqrt(counts_sum[:-1] + counts_sum[1:]) + valid_counts_sum = counts_sum[:, valid_bins] + diff = np.diff(valid_counts_sum, axis=0) - UltraConstants.SPECTRAL_SIG_THRESHOLD * ( + np.sqrt(valid_counts_sum[:-1] + valid_counts_sum[1:]) ) # shape (num_channels - 1, n_spin_bins) - flagged = np.any(diff > 0, axis=0) # shape (n_spin_bins,) + flagged = np.zeros(counts_sum.shape[1], dtype=bool) + flagged[valid_bins] = np.any(diff > 0, axis=0) # shape (n_spin_bins,) num_culled: int = np.sum(flagged) logger.info( f"Spectral culling removed {num_culled} spin bins using channels" @@ -1162,12 +1175,14 @@ def get_valid_events_per_energy_range( ) valid_outliers = de_dataset["quality_outliers"].values == 0 valid_scattering = de_dataset["quality_scattering"].values == 0 + print("SDC valid_outliers:", valid_outliers) + print("SDC valid_scattering:", valid_scattering) # TODO what about species non-proton? For those psets dont cull based on # High energy? ebin = de_dataset["ebin"].values valid_ebin = np.isin(ebin, UltraConstants.TOFXPH_SPECIES_GROUPS["proton"]) for i in range(len(energy_ranges) - 1): - energy_mask = (event_energies >= energy_ranges[i]) & ( + energy_mask = (event_energies > energy_ranges[i]) & ( event_energies < energy_ranges[i + 1] ) if not np.any(energy_mask): @@ -1194,7 +1209,9 @@ def get_valid_events_per_energy_range( valid_earth_angle, ] ) - + # SDC - inside get_valid_events_per_energy_range, after computing valid_events + for i in range(len(energy_ranges) - 1): + print(f"SDC energy bin {i} n_valid_events={np.sum(valid_events[i])}") return valid_events From b7fc767abba8fb3eec8d952e893bd0f765f9d313 Mon Sep 17 00:00:00 2001 From: Luisa Date: Mon, 4 May 2026 15:03:15 -0600 Subject: [PATCH 2/8] remove mask --- imap_processing/ultra/l1b/extendedspin.py | 1 - imap_processing/ultra/l1b/ultra_l1b_culling.py | 14 +++----------- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index d32e84c143..247a83e297 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -131,7 +131,6 @@ def calculate_extendedspin( de_dataset, spin_tbin_edges, energy_ranges, - mask, UltraConstants.SPECTRAL_ENERGY_CHANNELS, instrument_id, ) diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index ab2f63f73c..6d9bd34307 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -1027,7 +1027,6 @@ def flag_spectral_events( de_dataset: xr.Dataset, spin_tbin_edges: NDArray, energy_ranges: NDArray, - mask: NDArray, channels: list, sensor_id: int = 90, ) -> NDArray: @@ -1042,8 +1041,6 @@ def flag_spectral_events( Edges of the spin time bins. energy_ranges : NDArray Array of energy range edges. - mask : NDArray - Mask indicating which events to consider for spectral flagging. channels : list List of energy channel indices to use for spectral flagging. sensor_id : int @@ -1066,9 +1063,6 @@ def flag_spectral_events( counts_sum = get_valid_de_count_summary( de_dataset, energy_ranges, spin_tbin_edges, sensor_id=sensor_id )[channels, :] # shape (num_channels, n_spin_bins) - valid_bins = np.all( - ~mask[channels, :], axis=0 - ) # Get valid spin bins across the selected channels; shape (n_spin_bins,) # Flag spin bins where the signed count difference between adjacent selected # energy channels exceeds a Poisson-based threshold. For each pair of # adjacent channels, compute np.diff(counts_sum, axis=0) and compare that @@ -1076,12 +1070,10 @@ def flag_spectral_events( # uncertainty (sqrt(N1 + N2)) of those two channels for each spin bin. # If any adjacent channel pair exceeds the threshold for a spin bin, that # spin bin is flagged across all energy ranges. - valid_counts_sum = counts_sum[:, valid_bins] - diff = np.diff(valid_counts_sum, axis=0) - UltraConstants.SPECTRAL_SIG_THRESHOLD * ( - np.sqrt(valid_counts_sum[:-1] + valid_counts_sum[1:]) + diff = np.diff(counts_sum, axis=0) - UltraConstants.SPECTRAL_SIG_THRESHOLD * ( + np.sqrt(counts_sum[:-1] + counts_sum[1:]) ) # shape (num_channels - 1, n_spin_bins) - flagged = np.zeros(counts_sum.shape[1], dtype=bool) - flagged[valid_bins] = np.any(diff > 0, axis=0) # shape (n_spin_bins,) + flagged = np.any(diff > 0, axis=0) # shape (n_spin_bins,) num_culled: int = np.sum(flagged) logger.info( f"Spectral culling removed {num_culled} spin bins using channels" From ea49130db5202ab53a95292c3f6e7cb2daed11fc Mon Sep 17 00:00:00 2001 From: Luisa Date: Mon, 4 May 2026 15:51:13 -0600 Subject: [PATCH 3/8] remove print statements --- .../tests/ultra/unit/test_ultra_l1b_culling.py | 5 ++--- imap_processing/ultra/l1b/extendedspin.py | 4 ---- imap_processing/ultra/l1b/ultra_l1b_culling.py | 12 +----------- 3 files changed, 3 insertions(+), 18 deletions(-) diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py index 1df4845720..4945ab87bd 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py @@ -811,9 +811,8 @@ def test_flag_statistical_outliers(): # iterated twice. assert np.all(iterations[:-1] == 1) assert iterations[-1] == 2 - # Check that all std_diff values were set (not zero) except the last one - assert np.all(std_diff[:-1] != 0) - assert std_diff[-1] == 0 + # Check that all std_diff values are zero + assert np.all(std_diff[:-1] == 0) def test_flag_statistical_outliers_invalid_events(): diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index 247a83e297..67ea169325 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -75,9 +75,6 @@ def calculate_extendedspin( inst_qf = flag_imap_instruments(de_dataset["spin"].values) spin_bin_size = UltraConstants.SPIN_BIN_SIZE - print("SDC spin[:5]:", spin[:5]) - print("SDC spin_starttime[:5]:", spin_starttime[:5]) - print("SDC len(spin):", len(spin)) spin_tbin_edges = get_binned_spins_edges( spin, spin_period, spin_starttime, spin_bin_size ) @@ -94,7 +91,6 @@ def calculate_extendedspin( intervals, _, _ = build_energy_bins() # Get the energy ranges energy_ranges = get_binned_energy_ranges(intervals) - print("SDC energy_ranges:", energy_ranges) energy_bin_flags = get_energy_range_flags(energy_ranges) # Calculate the high energy quality flags energy_thresholds = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index 6d9bd34307..aff05be128 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -1015,11 +1015,6 @@ def flag_upstream_ion( f" independent flags and will be applied across all {mask.shape[0]} energy" f" channels." ) - # SDC - add inside flag_upstream_ion - # SDC - inside flag_upstream_ion, first call only - print("SDC upstream_1 total_scaled[:5]:", total_scaled[:5]) - print("SDC upstream_1 thresh:", thresh) - print("SDC upstream_1 flagged bins:", np.where(flagged)[0]) return flagged @@ -1167,14 +1162,12 @@ def get_valid_events_per_energy_range( ) valid_outliers = de_dataset["quality_outliers"].values == 0 valid_scattering = de_dataset["quality_scattering"].values == 0 - print("SDC valid_outliers:", valid_outliers) - print("SDC valid_scattering:", valid_scattering) # TODO what about species non-proton? For those psets dont cull based on # High energy? ebin = de_dataset["ebin"].values valid_ebin = np.isin(ebin, UltraConstants.TOFXPH_SPECIES_GROUPS["proton"]) for i in range(len(energy_ranges) - 1): - energy_mask = (event_energies > energy_ranges[i]) & ( + energy_mask = (event_energies >= energy_ranges[i]) & ( event_energies < energy_ranges[i + 1] ) if not np.any(energy_mask): @@ -1201,9 +1194,6 @@ def get_valid_events_per_energy_range( valid_earth_angle, ] ) - # SDC - inside get_valid_events_per_energy_range, after computing valid_events - for i in range(len(energy_ranges) - 1): - print(f"SDC energy bin {i} n_valid_events={np.sum(valid_events[i])}") return valid_events From c3bb966c76eb8d75542a920a8302b5d05693b09d Mon Sep 17 00:00:00 2001 From: Luisa Date: Mon, 4 May 2026 15:54:05 -0600 Subject: [PATCH 4/8] stat cull update --- imap_processing/ultra/l1b/ultra_l1b_culling.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index aff05be128..c12a0fa36c 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -844,7 +844,8 @@ def flag_statistical_outliers( quality_stats: np.ndarray = np.zeros((n_energy_bins, spin_bin_size), dtype=bool) # Initialize a mask to keep track of spin bins that have been flagged across all # energy bins - all_channel_mask: np.ndarray = np.zeros(spin_bin_size, dtype=bool) + # all_channel_mask: np.ndarray = np.zeros(spin_bin_size, dtype=bool) + all_channel_mask: np.ndarray = curr_mask[0, :].copy() # Initialize convergence array to keep track of poisson stats convergence = np.full(n_energy_bins, False) # Keep track of how many iterations we have done of flagging outliers and @@ -852,6 +853,7 @@ def flag_statistical_outliers( iterations = np.zeros(n_energy_bins) # keep track of the standard deviation difference from poisson stats per energy bin std_diff: np.ndarray = np.zeros(n_energy_bins, dtype=float) + count_summary = get_valid_de_count_summary( de_dataset, energy_ranges, spin_tbin_edges, sensor_id=sensor_id ) # shape (n_energy_bins, n_spin_bins) @@ -1194,6 +1196,7 @@ def get_valid_events_per_energy_range( valid_earth_angle, ] ) + return valid_events From 19a8c8c9f1eb3cb5c87b07f3e61e124f71909fba Mon Sep 17 00:00:00 2001 From: Luisa Date: Tue, 5 May 2026 16:11:35 -0600 Subject: [PATCH 5/8] clean up --- .../ultra/unit/test_ultra_l1b_culling.py | 17 +++++++---------- imap_processing/ultra/l1b/extendedspin.py | 7 ++++++- .../ultra/l1b/ultra_l1b_culling.py | 19 ++----------------- 3 files changed, 15 insertions(+), 28 deletions(-) diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py index 4945ab87bd..9eb37a9924 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py @@ -707,7 +707,6 @@ def test_flag_high_energy(): de_dataset, spin_tbin_edges, energy_range_edges, - None, cull_thresholds, 90, ) @@ -750,9 +749,7 @@ def test_validate_high_energy_cull(setup_repoint_47_data): de_ds, _, spin_tbin_edges = setup_repoint_47_data # Get the energy ranges energy_ranges = np.array([4.2, 9.4425, 21.2116, 47.2388, 105.202, 316.335]) - e_flags = flag_high_energy( - de_ds, spin_tbin_edges, energy_ranges, None, mock_thresholds - ) + e_flags = flag_high_energy(de_ds, spin_tbin_edges, energy_ranges, mock_thresholds) np.testing.assert_array_equal(e_flags, ~expected_qf.astype(bool)) @@ -812,7 +809,8 @@ def test_flag_statistical_outliers(): assert np.all(iterations[:-1] == 1) assert iterations[-1] == 2 # Check that all std_diff values are zero - assert np.all(std_diff[:-1] == 0) + assert np.all(std_diff[:-1] != 0) + assert std_diff[-1] == 0 def test_flag_statistical_outliers_invalid_events(): @@ -838,10 +836,9 @@ def test_flag_statistical_outliers_invalid_events(): energy_range_edges, mask, ) - # check that no flags are set because there were no valid events to calculate - # statistics on. + # check that all flags are set because the mask marks all events as invalid. np.testing.assert_array_equal( - quality_flags, np.zeros_like(quality_flags, dtype=bool) + quality_flags, np.ones_like(quality_flags, dtype=bool) ) # check that all energy bins are marked as converged (no valid events is not a # failure case for convergence since we just can't calculate statistics. @@ -872,11 +869,11 @@ def test_validate_stat_cull(setup_repoint_47_data): """Validate that statistical-outlier quality flags match expected results.""" # read test data from csv files results_df = pd.read_csv( - TEST_PATH / "validate_stat_culling_results_repoint00047_v2.csv" + TEST_PATH / "validate_stat_culling_results_repoint00047_v3.csv" ) de_ds, _, spin_tbin_edges = setup_repoint_47_data # Get the energy ranges - energy_ranges = np.array([4.2, 9.4425, 21.2116, 47.2388, 105.202, 316.335]) + energy_ranges = get_binned_energy_ranges(build_energy_bins()[0]) # Create a mask of flagged events to test that the stat cull algorithm # properly ignores these. The test data was created using this exact mask as well. diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index 67ea169325..d1749aac3c 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -98,7 +98,6 @@ def calculate_extendedspin( de_dataset, spin_tbin_edges, energy_ranges, - voltage_qf, energy_thresholds, instrument_id, ) @@ -176,6 +175,12 @@ def calculate_extendedspin( stop_per_spin[valid] = pulses.stop_per_spin[idx[valid]] coin_per_spin[valid] = pulses.coin_per_spin[idx[valid]] + # To be consistent with the ULTRA IT implementation, apply the low voltage mask + # to the high energy and upstream ion flags + upstream_ion_qf_2 |= voltage_qf + upstream_ion_qf_1 |= voltage_qf + high_energy_qf |= voltage_qf + spectral_qf |= voltage_qf # high energy and statistical outlier flags are energy dependent boolean arrays # with shape (n_energy_bins, n_spin_bins). We want to collapse the energy dimension # using a bitwise OR to get a single boolean flag per spin. diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index c12a0fa36c..da3d665a6c 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -687,7 +687,6 @@ def flag_high_energy( de_dataset: xr.Dataset, spin_tbin_edges: NDArray, energy_ranges: NDArray, - mask: NDArray = None, energy_thresholds: np.ndarray = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS, sensor_id: int = 90, ) -> NDArray: @@ -702,10 +701,6 @@ def flag_high_energy( Edges of the spin time bins. energy_ranges : numpy.ndarray Array of energy range edges. - mask : numpy.ndarray, optional - Mask indicating which events to consider for high energy flagging - (e.g., after low voltage culling). True indicates the spin bins that should - NOT be considered for high energy flagging. energy_thresholds : numpy.ndarray Array of count thresholds for flagging high energy events corresponding to each energy range. @@ -732,10 +727,6 @@ def flag_high_energy( f"HIGH_ENERGY_CULL_CHANNEL ({cull_channel}) is out of bounds" f" for {n_energy_bins} energy ranges." ) - - # Initialize all spin bins to have no high energy flag - spin_bin_size = len(spin_tbin_edges) - 1 - quality_flags: np.ndarray = np.zeros((n_energy_bins, spin_bin_size), dtype=bool) # Get valid events and counts at each spin bin for the # designated culling channel. de_counts = get_valid_de_count_summary( @@ -752,18 +743,12 @@ def flag_high_energy( cull_channel_counts[np.newaxis, :] >= energy_thresholds ) # (n_energy_bins, n_spin_bins) - if mask is not None: - quality_flags[:, ~mask] = flagged[:, ~mask] - else: - quality_flags = flagged - - num_culled: int = np.sum(quality_flags) logger.info( - f"High energy culling removed {num_culled} spin bins across {n_energy_bins} " + f"Found {np.sum(flagged)} high energy spin bins across {n_energy_bins} " f"energy channels. Energy thresholds: {energy_thresholds.flatten()}, " ) - return quality_flags + return flagged def flag_statistical_outliers( From e2a13005876558db99b1caaea2d4f0ea5a64f996 Mon Sep 17 00:00:00 2001 From: Luisa Date: Tue, 5 May 2026 16:31:51 -0600 Subject: [PATCH 6/8] add back masking logic to high energy --- .../ultra/unit/test_ultra_l1b_culling.py | 7 +++++-- .../ultra/l1b/ultra_l1b_culling.py | 21 +++++++++++++++---- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py index 9eb37a9924..94852f9b9b 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py @@ -707,6 +707,7 @@ def test_flag_high_energy(): de_dataset, spin_tbin_edges, energy_range_edges, + None, cull_thresholds, 90, ) @@ -749,7 +750,9 @@ def test_validate_high_energy_cull(setup_repoint_47_data): de_ds, _, spin_tbin_edges = setup_repoint_47_data # Get the energy ranges energy_ranges = np.array([4.2, 9.4425, 21.2116, 47.2388, 105.202, 316.335]) - e_flags = flag_high_energy(de_ds, spin_tbin_edges, energy_ranges, mock_thresholds) + e_flags = flag_high_energy( + de_ds, spin_tbin_edges, energy_ranges, None, mock_thresholds + ) np.testing.assert_array_equal(e_flags, ~expected_qf.astype(bool)) @@ -808,7 +811,7 @@ def test_flag_statistical_outliers(): # iterated twice. assert np.all(iterations[:-1] == 1) assert iterations[-1] == 2 - # Check that all std_diff values are zero + # Check that all std_diff values were set (not zero) except the last one assert np.all(std_diff[:-1] != 0) assert std_diff[-1] == 0 diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index da3d665a6c..3ec5c2a9b0 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -687,6 +687,7 @@ def flag_high_energy( de_dataset: xr.Dataset, spin_tbin_edges: NDArray, energy_ranges: NDArray, + mask: NDArray = None, energy_thresholds: np.ndarray = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS, sensor_id: int = 90, ) -> NDArray: @@ -701,6 +702,10 @@ def flag_high_energy( Edges of the spin time bins. energy_ranges : numpy.ndarray Array of energy range edges. + mask : numpy.ndarray, optional + Mask indicating which events to consider for high energy flagging + (e.g., after low voltage culling). True indicates the spin bins that should + NOT be considered for high energy flagging. energy_thresholds : numpy.ndarray Array of count thresholds for flagging high energy events corresponding to each energy range. @@ -727,6 +732,10 @@ def flag_high_energy( f"HIGH_ENERGY_CULL_CHANNEL ({cull_channel}) is out of bounds" f" for {n_energy_bins} energy ranges." ) + + # Initialize all spin bins to have no high energy flag + spin_bin_size = len(spin_tbin_edges) - 1 + quality_flags: np.ndarray = np.zeros((n_energy_bins, spin_bin_size), dtype=bool) # Get valid events and counts at each spin bin for the # designated culling channel. de_counts = get_valid_de_count_summary( @@ -743,12 +752,18 @@ def flag_high_energy( cull_channel_counts[np.newaxis, :] >= energy_thresholds ) # (n_energy_bins, n_spin_bins) + if mask is not None: + quality_flags[:, ~mask] = flagged[:, ~mask] + else: + quality_flags = flagged + + num_culled: int = np.sum(quality_flags) logger.info( - f"Found {np.sum(flagged)} high energy spin bins across {n_energy_bins} " + f"High energy culling removed {num_culled} spin bins across {n_energy_bins} " f"energy channels. Energy thresholds: {energy_thresholds.flatten()}, " ) - return flagged + return quality_flags def flag_statistical_outliers( @@ -829,7 +844,6 @@ def flag_statistical_outliers( quality_stats: np.ndarray = np.zeros((n_energy_bins, spin_bin_size), dtype=bool) # Initialize a mask to keep track of spin bins that have been flagged across all # energy bins - # all_channel_mask: np.ndarray = np.zeros(spin_bin_size, dtype=bool) all_channel_mask: np.ndarray = curr_mask[0, :].copy() # Initialize convergence array to keep track of poisson stats convergence = np.full(n_energy_bins, False) @@ -838,7 +852,6 @@ def flag_statistical_outliers( iterations = np.zeros(n_energy_bins) # keep track of the standard deviation difference from poisson stats per energy bin std_diff: np.ndarray = np.zeros(n_energy_bins, dtype=float) - count_summary = get_valid_de_count_summary( de_dataset, energy_ranges, spin_tbin_edges, sensor_id=sensor_id ) # shape (n_energy_bins, n_spin_bins) From b1021a6d70564afa39a71aa96f4e7b1131c07639 Mon Sep 17 00:00:00 2001 From: Luisa Date: Tue, 5 May 2026 16:34:06 -0600 Subject: [PATCH 7/8] add testfile --- imap_processing/ultra/l1b/extendedspin.py | 1 + 1 file changed, 1 insertion(+) diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index d1749aac3c..8fdd394999 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -98,6 +98,7 @@ def calculate_extendedspin( de_dataset, spin_tbin_edges, energy_ranges, + voltage_qf, energy_thresholds, instrument_id, ) From 0bfefbb7fc3feecefaf5e8a5876b62b63008d39d Mon Sep 17 00:00:00 2001 From: Luisa Date: Tue, 5 May 2026 16:34:56 -0600 Subject: [PATCH 8/8] add test --- imap_processing/tests/external_test_data_config.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imap_processing/tests/external_test_data_config.py b/imap_processing/tests/external_test_data_config.py index 3f6b3de60a..b962bb96d9 100644 --- a/imap_processing/tests/external_test_data_config.py +++ b/imap_processing/tests/external_test_data_config.py @@ -160,7 +160,7 @@ ("status_test_data_repoint00047.csv", "ultra/data/l1/"), ("voltage_culling_results_repoint00047.csv", "ultra/data/l1/"), ("validate_high_energy_culling_results_repoint00047_v2.csv", "ultra/data/l1/"), - ("validate_stat_culling_results_repoint00047_v2.csv", "ultra/data/l1/"), + ("validate_stat_culling_results_repoint00047_v3.csv", "ultra/data/l1/"), ("validate_upstream_ion_1_culling_results_repoint00047_v1.csv", "ultra/data/l1/"), ("validate_spectral_culling_results_repoint00047_v1.csv", "ultra/data/l1/"), ("de_test_data_repoint00047.csv", "ultra/data/l1/"),