diff --git a/imap_processing/cdf/config/imap_hit_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_hit_l1b_variable_attrs.yaml index 2525ef8769..193eceb8e2 100644 --- a/imap_processing/cdf/config/imap_hit_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hit_l1b_variable_attrs.yaml @@ -520,6 +520,14 @@ dynamic_threshold_state: VALIDMIN: 0 VAR_TYPE: data +hdr_minute_cnt: + <<: *default + CATDESC: Raw header minute count per integration + DISPLAY_TYPE: no_plot + FIELDNAM: Header Minute Count + LABLAXIS: Minute count + VAR_TYPE: support_data + # <=== Standard Rates Data Variable Attributes ===> diff --git a/imap_processing/hit/l1a/hit_l1a.py b/imap_processing/hit/l1a/hit_l1a.py index 3d6c1f9411..28070af19d 100644 --- a/imap_processing/hit/l1a/hit_l1a.py +++ b/imap_processing/hit/l1a/hit_l1a.py @@ -152,7 +152,7 @@ def subcom_sectorates(sci_dataset: xr.Dataset) -> xr.Dataset: for i, mod_10 in enumerate(hdr_min_count_mod_10): # NOTE: this ignore is needed to avoid a mypy error in Git # tests. - data_by_species_and_energy_range[mod_10]["counts"][i] = updated_dataset[ # type: ignore[index] + data_by_species_and_energy_range[mod_10]["counts"][i] = updated_dataset[ # type: ignore[index, unused-ignore] "sectorates" ].values[i] @@ -336,46 +336,53 @@ def add_cdf_attributes( return dataset -def find_complete_mod10_sets(mod_vals: np.ndarray) -> np.ndarray: +def find_mod10_set_ids(hdr_minute_cnt: np.ndarray) -> np.ndarray: """ - Find start indices where mod values match [0,1,2,3,4,5,6,7,8,9] pattern. + Assign a set id to each frame of sectored counts. + + A set is nominally 10 frames (10 minutes) with mod 10 values 0-9. A new set + starts wherever the mod 10 value decreases. Sets can be incomplete, e.g. 9 + frames, if a frame was dropped. Parameters ---------- - mod_vals : np.ndarray - 1D array of mod 10 values from the hdr_minute_cnt field in the L1A counts data. + hdr_minute_cnt : numpy.ndarray + The header minute count data. Returns ------- np.ndarray - Indices in mod_vals where the complete pattern [0, 1, ..., 9] starts. + Set id for each frame in the dataset, starting at 0. """ - # The pattern to match is an array from 0-9 - window_size = 10 + # Identify 10-minute intervals of sectored counts + # by using the mod 10 values of the header minute counts. + # Mod 10 determines the species and energy bin the data belongs + # to. A mapping of mod 10 values to species and energy bins is + # provided in l0/constants.py for reference. + mod_10: np.ndarray = hdr_minute_cnt % 10 - if mod_vals.size < window_size: - logger.warning( - "Mod 10 array is smaller than the required window size for " - "pattern matching." - ) - return np.array([], dtype=int) - # Use sliding windows to find pattern matches - sw_view = np.lib.stride_tricks.sliding_window_view(mod_vals, window_size) - matches = np.all(sw_view == np.arange(window_size), axis=1) - return np.where(matches)[0] + # Assign frame a set ID + # E.g. if we have 0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9 + # The set_ids would be [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1] + # A set can be incomplete if a frame was dropped so we can not just look for + # complete sets [0-9]. + ids: np.ndarray = np.cumsum(np.diff(mod_10, prepend=0) < 0) + # TODO group by id and check if any epochs are greater than 10 mins apart. + # This means that 2 different sets were incorrectly identified as 1. + return ids def subset_sectored_counts( sectored_counts_dataset: xr.Dataset, packet_date: str ) -> xr.Dataset: """ - Subset data for complete sets of sectored counts and corresponding livetime values. + Subset data for sets of sectored counts and corresponding livetime values. A set of sectored data starts with hydrogen and ends with iron and correspond to - the mod 10 values 0-9. The livetime values from the previous 10 minutes are used + the mod 10 values 0-9. The livetime values from the previous set are used to calculate the rates for each set since those counts are transmitted 10 minutes - after they were collected. Therefore, only complete sets of sectored counts where - livetime from the previous 10 minutes are available are included in the output. + after they were collected. Sets with a dropped frame are kept. The first set is + only used for the livetime values of the next set. Parameters ---------- @@ -388,62 +395,55 @@ def subset_sectored_counts( Returns ------- xarray.Dataset - A dataset of complete sectored counts and corresponding livetime values + A dataset of sectored counts and corresponding livetime values for the processing day. """ # Modify livetime_counter to use a new epoch coordinate # that is aligned with the original epoch dimension. This # ensures that livetime doesn't get filtered when the original - # epoch dimension is filtered for complete sets. + # epoch dimension is filtered for sets. sectored_counts_dataset = update_livetime_coord(sectored_counts_dataset) - # Identify 10-minute intervals of complete sectored counts - # by using the mod 10 values of the header minute counts. - # Mod 10 determines the species and energy bin the data belongs - # to. A mapping of mod 10 values to species and energy bins is - # provided in l0/constants.py for reference. - bin_size = 10 - mod_10: np.ndarray = sectored_counts_dataset.hdr_minute_cnt.values % 10 - start_indices = find_complete_mod10_sets(mod_10) - - # Filter out start indices that are less than or equal to the bin size - # since the previous 10 minutes are needed for calculating rates + # Get set id for each frame + # The first set is set -1 and is only used for livetime values. + frame_set_ids: np.ndarray = ( + find_mod10_set_ids(sectored_counts_dataset.hdr_minute_cnt.values) - 1 + ) + start_indices = np.flatnonzero(np.diff(frame_set_ids) > 0) + 1 if start_indices.size == 0: raise ValueError( - "No data to process - valid start indices not found for " - "complete sectored counts." + "No data to process - valid start indices not found for sectored counts." ) - else: - start_indices = start_indices[start_indices >= bin_size] - # Subset data for complete sets of sectored counts. - # Each set of sectored counts is 10 minutes long, so we take the indices - # starting from the start indices and extending to the bin size of 10. - # This creates a 1D array of indices that correspond to the complete + # Subset data for sets of sectored counts. + # Each set of sectored counts is generally expected to be 10 minutes long, + # but we can have incomplete sets where a frame was dropped so we take the indices + # starting from the first start index to the end. + # This creates a 1D array of indices that correspond to the # sets of sectored counts which is used to filter the L1A dataset and # create the L1B sectored rates dataset. - data_indices = np.concatenate( - [np.arange(idx, idx + bin_size) for idx in start_indices] + complete_sectored_counts_dataset = sectored_counts_dataset.isel( + epoch=slice(start_indices[0], None) ) - complete_sectored_counts_dataset = sectored_counts_dataset.isel(epoch=data_indices) - epoch_per_complete_set: np.ndarray = np.repeat( - [ - complete_sectored_counts_dataset.epoch[idx : idx + bin_size].mean().item() - for idx in range(0, len(complete_sectored_counts_dataset.epoch), 10) - ], - bin_size, - ) + filtered_frame_epochs = complete_sectored_counts_dataset.epoch.values + filtered_frame_set_ids = frame_set_ids[start_indices[0] :] - # Filter dataset for data in the processing day + # Get the mean epoch of each set + set_counts = np.unique(filtered_frame_set_ids, return_counts=True)[1] + set_sums = np.bincount(filtered_frame_set_ids, weights=filtered_frame_epochs) + mean_epoch_per_set: np.ndarray = np.repeat(set_sums / set_counts, set_counts) - # Trim the sectored data to epoch_per_complete_set values in the processing day + # Filter dataset for data in the processing day + # Trim the sectored data to epoch_per_set values in the processing day filtered_dataset = filter_dataset_to_processing_day( - complete_sectored_counts_dataset, packet_date, epoch_vals=epoch_per_complete_set + complete_sectored_counts_dataset, packet_date, epoch_vals=mean_epoch_per_set ) - # Trim livetime to the size of the sectored data but shifted 10 minutes earlier. - filtered_dataset = subset_livetime(filtered_dataset) + # Keep the livetime values from the set before each set of sectored counts + filtered_dataset = subset_livetime( + filtered_dataset, sectored_counts_dataset.epoch.values, frame_set_ids + ) return filtered_dataset @@ -455,7 +455,7 @@ def update_livetime_coord(sectored_dataset: xr.Dataset) -> xr.Dataset: Assign a new epoch coordinate to the `livetime_counter` variable. This new coordinate is aligned with the original `epoch` dimension, ensuring that `livetime_counter` remains unaffected when the original - `epoch` dimension is filtered for complete sets in `subset_sectored_counts` + `epoch` dimension is filtered for sets in `subset_sectored_counts` function. Parameters @@ -486,47 +486,46 @@ def update_livetime_coord(sectored_dataset: xr.Dataset) -> xr.Dataset: return sectored_dataset -def subset_livetime(dataset: xr.Dataset) -> xr.Dataset: +def subset_livetime( + dataset: xr.Dataset, all_frame_epochs: np.ndarray, frame_set_ids: np.ndarray +) -> xr.Dataset: """ - Trim livetime to values shifted 10 minutes earlier than sectored data. + Trim livetime to the values from the set before each set of sectored counts. - This ensures that the livetime values correspond to the sectored counts correctly - since sectored data is collected 10 minutes before it's transmitted. + Sectored counts are collected over 10 minutes and transmitted over the next + 10 minutes, so the counts in set k use the livetime values from the frames + in set k-1. Livetime is matched by set number rather than by position so + that sets with missing frames are paired with the correct livetime values. Parameters ---------- dataset : xarray.Dataset The dataset containing sectored counts and livetime data. + all_frame_epochs : np.ndarray + Epoch values for all frames, including the set before the first set of + sectored counts. + frame_set_ids : np.ndarray + Set number for each value in `frame_epochs`. Returns ------- xarray.Dataset The updated dataset with trimmed livetime data. """ - # epoch values are per science frame which is 1 minute - epoch_vals = dataset["epoch"].values - epoch_livetime_vals = dataset["epoch_livetime"].values - - if not epoch_vals.size: + if not dataset["epoch"].size: raise ValueError( "Epoch values are empty. Cannot proceed with livetime subsetting." ) - - # Get index positions of epoch[0] and epoch[-1] in epoch_livetime - start_idx = np.where(epoch_livetime_vals == epoch_vals[0])[0][0] - end_idx = np.where(epoch_livetime_vals == epoch_vals[-1])[0][0] - - if start_idx < 10: - raise ValueError( - "Start index for livetime is less than 10. This indicates that the " - "dataset is too small to shift livetime correctly." - ) - - # Compute shifted indices by 10 minutes - start_trimmed = max(start_idx - 10, 0) - end_trimmed = max(end_idx - 10, 0) - - return dataset.isel(epoch_livetime=slice(start_trimmed, end_trimmed + 1)) + # Determine for each frame whether the epoch was kept (in the processing day) + kept_epochs = np.isin(all_frame_epochs, dataset["epoch"].values) + # get all the kept set ids + kept_sets = np.unique(frame_set_ids[kept_epochs]) + # Livetime values for the set comes from the frames in the previous set + livetime_sets = np.isin(frame_set_ids, kept_sets - 1) + livetime_ids = np.flatnonzero( + np.isin(dataset["epoch_livetime"].values, all_frame_epochs[livetime_sets]) + ) + return dataset.isel(epoch_livetime=livetime_ids) def filter_dataset_to_processing_day( diff --git a/imap_processing/hit/l1b/hit_l1b.py b/imap_processing/hit/l1b/hit_l1b.py index 05d5d2f214..ea2580bfbd 100644 --- a/imap_processing/hit/l1b/hit_l1b.py +++ b/imap_processing/hit/l1b/hit_l1b.py @@ -13,6 +13,7 @@ get_datasets_by_apid, process_housekeeping_data, ) +from imap_processing.hit.l1a.hit_l1a import find_mod10_set_ids from imap_processing.hit.l1b.constants import ( FILLVAL_FLOAT32, FILLVAL_INT64, @@ -298,27 +299,57 @@ def calculate_rates( return dataset -def sum_livetime_10min(livetime: xr.DataArray) -> xr.DataArray: +def sum_livetime_10min( + hdr_minute_cnt: xr.DataArray, livetime: xr.DataArray +) -> xr.DataArray: """ - Sum livetime values in 10-minute intervals. + Sum livetime values over the set before each set of sectored counts. + + In L1A, livetime is subset to the frames of the set before each set of + sectored counts, so the livetime set sizes are the sectored counts set sizes + shifted one set later. If the previous set is missing frames, the livetime + sum is scaled by 10 / (number of livetime frames in the set). Parameters ---------- + hdr_minute_cnt : xr.DataArray + 1D array of header minute count values on the epoch dimension. livetime : xr.DataArray - 1D array of livetime values. Shape equals the number of epochs in the dataset. + 1D array of livetime values on the epoch_livetime dimension. Returns ------- xr.DataArray - Livetime summed over 10-minute intervals. Values repeated for each epoch in the - 10-minute intervals to match the original livetime array shape. + Livetime summed over the previous set for each sectored counts frame, + on the epoch dimension. E.g. for sets of 10 frames [5,5,5,5,5,5,5,5,5,5, 6,6,6,6,6,6,6,6,6,6, 7,7,7,7,7,7,7,7,7,7]. """ - livetime_10min_sum = [ - livetime[i : i + 10].sum().item() for i in range(0, len(livetime) - 9, 10) - ] - livetime_expanded: np.ndarray = np.repeat(livetime_10min_sum, 10) - return xr.DataArray(livetime_expanded, dims=livetime.dims, coords=livetime.coords) + # Set ids for each sectored counts frame + counts_set_ids = find_mod10_set_ids(hdr_minute_cnt.values) + set_counts = np.bincount(counts_set_ids) + + # Livetime for set k comes from set k-1, so the livetime set sizes are the + # counts set sizes shifted one set later. The first set's livetime is the + # set before the first set of counts. The first set's size is the remainder counts. + first_set_count = len(livetime) - set_counts[:-1].sum() + if not 0 < first_set_count <= 10: + raise ValueError("Livetime does not line up with the sectored counts sets.") + livetime_counts = np.insert(set_counts[:-1], 0, first_set_count) + # Once you have the shifted counts, recreate the set ids for each frame. + livetime_set_ids: np.ndarray = np.repeat( + np.arange(set_counts.size), livetime_counts + ) + + # Sum livetime for each set and scale for missing frames + livetime_sum = np.bincount(livetime_set_ids, weights=livetime.values) + livetime_sum = livetime_sum * (10 / livetime_counts) + + # Array with one scaled livetime sum per counts frame + return xr.DataArray( + livetime_sum[counts_set_ids], + dims=hdr_minute_cnt.dims, + coords=hdr_minute_cnt.coords, + ) def process_summed_rates_data( @@ -433,9 +464,8 @@ def process_sectored_rates_data( for var in l1a_counts_dataset.data_vars if any(str(var).startswith(f"{p}_") for p in particles) ] - # Sum livetime over 10 minute intervals - livetime_10min = sum_livetime_10min(livetime) + livetime_10min = sum_livetime_10min(l1a_counts_dataset.hdr_minute_cnt, livetime) # Initialize the L1B dataset with coordinates from the subset L1A dataset l1b_sectored_rates_dataset = initialize_l1b_dataset( @@ -451,6 +481,8 @@ def process_sectored_rates_data( "fe_energy_mean", ], ) + # Keep the header minute count so sets of sectored data can be identified in L2 + l1b_sectored_rates_dataset["hdr_minute_cnt"] = l1a_counts_dataset["hdr_minute_cnt"] # Dictionary to store variable rename mappings for L1B dataset rename_map = {} diff --git a/imap_processing/hit/l2/hit_l2.py b/imap_processing/hit/l2/hit_l2.py index 557a82be3b..cbadd91180 100644 --- a/imap_processing/hit/l2/hit_l2.py +++ b/imap_processing/hit/l2/hit_l2.py @@ -11,6 +11,7 @@ add_summed_particle_data_to_dataset, get_attribute_manager, ) +from imap_processing.hit.l1a.hit_l1a import find_mod10_set_ids from imap_processing.hit.l2.constants import ( FILLVAL_FLOAT32, N_AZIMUTH, @@ -810,16 +811,17 @@ def transform_to_10_minute_chunks(macropixel_dataset: xr.Dataset) -> xr.Dataset: ("nemgsi", 2), ("fe", 1), ] - - # Use the first record in each 10-record group as the output template. + hr_minute_count = macropixel_dataset.hdr_minute_cnt.values + set_ids = find_mod10_set_ids(hr_minute_count) + set_starts = np.flatnonzero(np.diff(set_ids, prepend=-1)) + # Use the first record in each set as the output template. transformed_dataset = macropixel_dataset.isel( - epoch=slice(None, None, 10), + epoch=set_starts, ).copy(deep=True) - # Each minute in a 10-record group contains one species/energy combination, - # ordered as described by species_energy. Track that minute's packet offset. - species_i = 0 - for species, num_energy_levels in species_energy: + # Each minute in a set contains one species/energy combination and every + # other species/energy is a fill value. + for species, _ in species_energy: energy_dim = f"{species}_energy_mean" # Gather the intensity and uncertainty variables that share this # species' energy dimension. @@ -829,21 +831,34 @@ def transform_to_10_minute_chunks(macropixel_dataset: xr.Dataset) -> xr.Dataset: if macropixel_dataset[var].dims[:2] == ("epoch", energy_dim) ] - for energy_i in range(num_energy_levels): - for var in species_variables: - # Select this species/energy packet from every 10-record group - # and place it in the corresponding output energy plane. - data_i = macropixel_dataset[var].values[species_i::10, energy_i] - transformed_dataset[var].values[:, energy_i] = data_i - species_i += 1 + for var in species_variables: + # Replace fill vals with nans so the skipna works. + da = macropixel_dataset[var] + da = xr.where(da == FILLVAL_FLOAT32, np.nan, da) + # Group by set_ids and take the first non NaN value in each group. + # (There will only be one row that isn't fill so taking the first is + # safe). If a frame was dropped, that value stays as fill. + transformed_dataset[var].values = ( + da.groupby(xr.DataArray(set_ids, dims="epoch")) + .first(skipna=True) + .fillna(FILLVAL_FLOAT32) + .values + ) minute_cadence_epochs = macropixel_dataset["epoch"].values - ten_minute_cadence_epochs = minute_cadence_epochs.reshape(-1, 10) nanoseconds_per_10_min = SECONDS_PER_10_MIN * 1_000_000_000 + # The header minute count mod 10 gives the minutes past the set start. + # Convert this to nanoseconds. + nanoseconds_since_set_start = (hr_minute_count % 10) * 60 * 1_000_000_000 nanoseconds_per_5_min = nanoseconds_per_10_min // 2 - start_times = ten_minute_cadence_epochs[:, 0] - end_times = ten_minute_cadence_epochs[:, -1] - new_epochs = start_times + (end_times - start_times) // 2 - nanoseconds_per_10_min + # To find the actual set starts in case of a dropped frame, subtract the time + # since the set started. E.g. if the first frame in a set is dropped, the first + # remaining frame is not the actual set start time. If the first frame wasn't + # dropped, its time since set start is 0. + calculated_set_start_times = minute_cadence_epochs - nanoseconds_since_set_start + new_epochs = ( + calculated_set_start_times + nanoseconds_per_5_min - nanoseconds_per_10_min + )[set_starts] transformed_dataset = transformed_dataset.assign_coords(epoch=np.array(new_epochs)) transformed_dataset["epoch_delta"] = xr.DataArray( @@ -851,4 +866,5 @@ def transform_to_10_minute_chunks(macropixel_dataset: xr.Dataset) -> xr.Dataset: dims=["epoch"], ) - return transformed_dataset + # hdr_minute_cnt was only needed to find the sets, so drop it from the L2 product + return transformed_dataset.drop_vars("hdr_minute_cnt") diff --git a/imap_processing/tests/hit/test_hit_l1a.py b/imap_processing/tests/hit/test_hit_l1a.py index 1bf7ab8156..7b9633cbeb 100644 --- a/imap_processing/tests/hit/test_hit_l1a.py +++ b/imap_processing/tests/hit/test_hit_l1a.py @@ -16,7 +16,7 @@ calculate_uncertainties, decom_hit, filter_dataset_to_processing_day, - find_complete_mod10_sets, + find_mod10_set_ids, hit_l1a, subcom_sectorates, subset_livetime, @@ -140,31 +140,31 @@ def test_update_livetime_coord(): def test_subset_livetime(): """Test the subset_livetime function.""" - # Test case 1: Normal case with valid epoch and epoch_livetime sizes + # Test case 1: Complete sets - # epoch_livetime goes from 0 to 59 (simulates 60 minutes) - epoch_livetime = np.arange(60) + # 60 frames (60 minutes) in 6 sets of 10 frames + all_frame_epochs = np.arange(60) + frame_set_ids = np.repeat(np.arange(6), 10) - # epoch is 30 minute subset of epoch_livetime - epoch = epoch_livetime[10:40] + # Sectored counts kept for sets 2 and 3 + epoch = all_frame_epochs[20:40] dataset = xr.Dataset( { "livetime_counter": xr.DataArray( np.arange(60), dims=["epoch_livetime"], - coords={"epoch_livetime": epoch_livetime}, + coords={"epoch_livetime": all_frame_epochs}, ) }, - coords={"epoch": epoch, "epoch_livetime": epoch_livetime}, + coords={"epoch": epoch, "epoch_livetime": all_frame_epochs}, ) # Call the function - trimmed = subset_livetime(dataset) + trimmed = subset_livetime(dataset, all_frame_epochs, frame_set_ids) - # The trimmed livetime should match the slice from index 0 to 29 (inclusive) - expected_indices = slice(0, 30) - expected_livetime = dataset.isel(epoch_livetime=expected_indices) + # The livetime should come from the previous sets (sets 1 and 2) + expected_livetime = dataset.isel(epoch_livetime=slice(10, 30)) # Validate the result np.testing.assert_array_equal( @@ -174,7 +174,22 @@ def test_subset_livetime(): trimmed["livetime_counter"].values, expected_livetime.livetime_counter.values ) - # Test case 2: epoch is empty, the function should raise a ValueError + # Test case 2: Dropped frame in set 1. Set 2 is complete. + all_frame_epochs_dropped = np.delete(all_frame_epochs, 14) + frame_set_ids_dropped = np.delete(frame_set_ids, 14) + dataset_dropped = dataset.sel(epoch_livetime=all_frame_epochs_dropped) + + trimmed = subset_livetime( + dataset_dropped, all_frame_epochs_dropped, frame_set_ids_dropped + ) + + # Set 2 uses the 9 frames from set 1 and set 3 uses the 10 frames from set 2 + expected_epoch_livetime = np.delete(np.arange(10, 30), 4) + np.testing.assert_array_equal( + trimmed["epoch_livetime"].values, expected_epoch_livetime + ) + + # Test case 3: epoch is empty, the function should raise a ValueError dataset = xr.Dataset( { "livetime_counter": ("epoch_livetime", np.array([10, 20, 30])), @@ -188,60 +203,23 @@ def test_subset_livetime(): ValueError, match="Epoch values are empty. Cannot proceed with livetime subsetting.", ): - subset_livetime(dataset) - - # Test case 3: Not enough livetime values, the function should raise a ValueError - dataset = xr.Dataset( - { - "livetime_counter": ("epoch_livetime", np.array([10, 20, 30, 40, 50])), - }, - coords={ - "epoch": np.array([1, 2, 3]), - "epoch_livetime": np.array([0, 1, 2, 3, 4]), - }, - ) - with pytest.raises( - ValueError, - match="Start index for livetime is less than 10. This indicates that the " - "dataset is too small to shift livetime correctly.", - ): - subset_livetime(dataset) + subset_livetime(dataset, np.array([0, 1, 2]), np.array([0, 0, 0])) -def test_find_complete_mod10_sets(): - """Test the find_complete_mod10_sets function.""" +def test_find_mod10_set_ids(): + """Test the find_mod10_set_ids function.""" - # Test case 1: Valid pattern exists throughout - mod_vals = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) - expected_indices = np.array([0, 10]) - result = find_complete_mod10_sets(mod_vals) - assert np.array_equal(result, expected_indices), ( - f"Expected {expected_indices}, got {result}" - ) + # Test with complete sets + result = find_mod10_set_ids(np.arange(100, 120)) + np.testing.assert_array_equal(result, np.repeat([0, 1], 10)) - # Test case 2: No valid pattern - mod_vals = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6]) - expected_indices = np.array([]) - result = find_complete_mod10_sets(mod_vals) - assert np.array_equal(result, expected_indices), ( - f"Expected {expected_indices}, got {result}" - ) - - # Test case 3: Pattern in the middle - mod_vals = np.array([5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4]) - expected_indices = np.array([5]) - result = find_complete_mod10_sets(mod_vals) - assert np.array_equal(result, expected_indices), ( - f"Expected {expected_indices}, got {result}" - ) + # Test with dropped frame in the second set + result = find_mod10_set_ids(np.delete(np.arange(100, 130), 14)) + np.testing.assert_array_equal(result, np.repeat([0, 1, 2], [10, 9, 10])) - # Test case 4: Empty input - mod_vals = np.array([]) - expected_indices = np.array([]) - result = find_complete_mod10_sets(mod_vals) - assert np.array_equal(result, expected_indices), ( - f"Expected {expected_indices}, got {result}" - ) + # Test with data starting in the middle of a set + result = find_mod10_set_ids(np.arange(105, 120)) + np.testing.assert_array_equal(result, np.repeat([0, 1], [5, 10])) def test_calculate_uncertainties(): @@ -382,6 +360,8 @@ def test_filter_dataset_to_processing_day(): def test_subset_sectored_counts(): """Test the subset_sectored_counts function.""" + epoch_start = 316008084684000000 + def create_l1a_counts_dataset(hdr_minute_cnt_values): """Helper to create L1A counts dataset.""" return xr.Dataset( @@ -393,27 +373,34 @@ def create_l1a_counts_dataset(hdr_minute_cnt_values): }, coords={ "epoch": np.arange( - 316008084684000000, 316008084684000000 + len(hdr_minute_cnt_values) + epoch_start, epoch_start + len(hdr_minute_cnt_values) ), }, ) - def validate_subset(l1a_counts_dataset): + def validate_subset(l1a_counts_dataset, expected_epoch, expected_livetime): """Helper to validate the subset results.""" # Define the packet date packet_date = "20100106" subset_dataset = subset_sectored_counts(l1a_counts_dataset, packet_date) - assert subset_dataset.sizes["epoch"] == 10 - assert subset_dataset.sizes["epoch_livetime"] == 10 - assert np.all(subset_dataset["hdr_minute_cnt"].values % 10 == np.arange(10)) - assert np.all( - subset_dataset["epoch_livetime"].values - == subset_dataset["epoch"].values - 10 - ), "epoch_livetime values are not shifted by 10 from epoch values" + np.testing.assert_array_equal( + subset_dataset["epoch"].values, + expected_epoch, + err_msg="epoch values are incorrect.", + ) + np.testing.assert_array_equal( + subset_dataset["epoch_livetime"].values, + expected_livetime, + err_msg="epoch_livetime values are incorrect.", + ) # Test with partial data at the start and end of the dataset + # Sets: [0-4] [5-14] [15-24] [25-29] + # Counts are every set but the first; livetime is every set but the last. l1a_counts_dataset = create_l1a_counts_dataset(np.arange(105, 135)) - validate_subset(l1a_counts_dataset) + expected_epoch = epoch_start + np.arange(5, 30) + expected_livetime = epoch_start + np.arange(0, 25) + validate_subset(l1a_counts_dataset, expected_epoch, expected_livetime) # Test with partial data in the middle of the dataset l1a_counts_dataset = create_l1a_counts_dataset( @@ -450,7 +437,10 @@ def validate_subset(l1a_counts_dataset): 139, ] ) - validate_subset(l1a_counts_dataset) + # Sets: [0-9] [10-14] [15-19] [20-29] + expected_epoch = epoch_start + np.arange(10, 30) + expected_livetime = epoch_start + np.arange(0, 20) + validate_subset(l1a_counts_dataset, expected_epoch, expected_livetime) # Test with partial data at the start, middle, and end of the dataset l1a_counts_dataset = create_l1a_counts_dataset( @@ -487,16 +477,86 @@ def validate_subset(l1a_counts_dataset): 141, ] ) - validate_subset(l1a_counts_dataset) + # Sets: [0-4] [5-14] [15-17] [18-27] [28-29] + expected_epoch = epoch_start + np.arange(5, 30) + expected_livetime = epoch_start + np.arange(0, 28) + validate_subset(l1a_counts_dataset, expected_epoch, expected_livetime) - # Test with only partial data in the dataset + # Test with only partial data in the dataset (every other frame dropped) + # Sets: [0-4] [5-9] [10-14] [15-19] [20-24] [25-29] l1a_counts_dataset = create_l1a_counts_dataset(np.arange(100, 160, 2)) - with pytest.raises( - ValueError, - match="No data to process - valid start indices not found for " - "complete sectored counts.", - ): - subset_sectored_counts(l1a_counts_dataset, packet_date="20100106") + expected_epoch = epoch_start + np.arange(5, 30) + expected_livetime = epoch_start + np.arange(0, 25) + validate_subset(l1a_counts_dataset, expected_epoch, expected_livetime) + + +def test_subset_sectored_counts_dropped_frame(): + """Test that a dropped frame in the sectored counts data is handled correctly.""" + # Create a mock dataset with a dropped frame + # Define 4 complete sets + hdr_minute_cnt_values = np.arange(100, 140) + # Simulate a dropped frame by removing a value in the third set + hdr_minute_cnt_values = np.delete(hdr_minute_cnt_values, 25) + epoch = np.concatenate( + [ + np.arange(20), # First 20 are outside the processing date + np.arange( + 316008084684000000, 316008084684000000 + 19 + ), # Next 19 are within + ] + ) + l1a_counts_dataset = xr.Dataset( + { + "hdr_minute_cnt": ("epoch", hdr_minute_cnt_values), + "h_sectored_counts": ("epoch", np.arange(len(hdr_minute_cnt_values))), + "he4_sectored_counts": ("epoch", np.arange(len(hdr_minute_cnt_values))), + "livetime_counter": ("epoch", np.arange(len(hdr_minute_cnt_values))), + }, + coords={"epoch": ("epoch", epoch)}, + ) + date = "20100106" + + # Call the subset_sectored_counts function + subset_dataset = subset_sectored_counts(l1a_counts_dataset, date) + + # Assert that the dataset has the correct number of epochs + # Should be one full set (10) plus a set with a missing frame (9) + assert subset_dataset.sizes["epoch"] == 19 + assert subset_dataset.sizes["epoch_livetime"] == 19 + + # Assert that the hdr_minute_cnt values are correct + # - the first two sets are skipped because they not in the processing date. + # - The third set is included but has a dropped frame, so it should have 9 values. + # - The final complete set is included + expected_hdr_minute_cnt = hdr_minute_cnt_values[20:] + assert np.array_equal( + subset_dataset["hdr_minute_cnt"].values, expected_hdr_minute_cnt + ) + # The livetime values should come from the middle two sets. + expected_livetime_epoch = epoch[10:29] + assert np.array_equal( + subset_dataset["epoch_livetime"].values, expected_livetime_epoch + ) + + +def test_subset_sectored_counts_not_enough_data(): + """Test that an error is raised when there are not enough sets.""" + # Create a mock dataset with a dropped frame + # Define 1 complete set (not enough to process) + hdr_minute_cnt_values = np.arange(100, 110) + epoch = np.arange(316008084684000000, 316008084684000000 + 10) + l1a_counts_dataset = xr.Dataset( + { + "hdr_minute_cnt": ("epoch", hdr_minute_cnt_values), + "h_sectored_counts": ("epoch", np.arange(len(hdr_minute_cnt_values))), + "he4_sectored_counts": ("epoch", np.arange(len(hdr_minute_cnt_values))), + "livetime_counter": ("epoch", np.arange(len(hdr_minute_cnt_values))), + }, + coords={"epoch": ("epoch", epoch)}, + ) + date = "20100106" + with pytest.raises(ValueError, match="No data to process"): + subset_sectored_counts(l1a_counts_dataset, date) def test_validate_l1a_housekeeping_data(hk_packet_filepath): @@ -609,13 +669,13 @@ def test_validate_l1a_counts_data(sci_packet_filepath, validation_data): # Copy the validation data for sectored data validation. # The first complete set of sectored data with sufficient livetime - # data available is from index 17 to -4. Slice the validation data. - sectored_validation_data = validation_data.iloc[17:-4].copy().reset_index(drop=True) + # data available is from index 7 to the end. Slice the validation data. + sectored_validation_data = validation_data.iloc[7:].copy().reset_index(drop=True) - # The corresponding livetime values for the sectored data is from index 7 to -14 + # The corresponding livetime values for the sectored data is from index 0 to -4 # (i.e. 10 minutes before the first complete set of sectored data). livetime_validation_data = ( - validation_data[["livetime_counter"]].iloc[7:-14].copy().reset_index(drop=True) + validation_data[["livetime_counter"]].iloc[0:-4].copy().reset_index(drop=True) ) # NOTE: slicing indices are specific to the sci_sample_raw.csv validation file diff --git a/imap_processing/tests/hit/test_hit_l1b.py b/imap_processing/tests/hit/test_hit_l1b.py index 84b5ed1dfb..3b75a2378b 100644 --- a/imap_processing/tests/hit/test_hit_l1b.py +++ b/imap_processing/tests/hit/test_hit_l1b.py @@ -134,14 +134,17 @@ def test_calculate_rates(): ) -def test_sum_livetime_10min(): - """Test the sum_livetime_10min function.""" - # Create a sample livetime DataArray +def test_sum_livetime_10min_3_complete_sets(): + """Test the sum_livetime_10min function with 3 complete sets.""" + # Create a sample livetime DataArray of 3 complete sets livetime_values = np.arange(1, 31) # 30 epochs with values 1 to 30 livetime = xr.DataArray( livetime_values, dims=["epoch"], coords={"epoch": np.arange(30)} ) - + # Create a sample hdr_minute_cnt array that has 3 complete sets + hdr_minute_cnt = xr.DataArray( + np.arange(10, 40), dims="epoch", coords={"epoch": np.arange(30)} + ) # Expected result: sum of every 10 values repeated 10 times expected_values = np.repeat( [sum(livetime_values[i : i + 10]) for i in range(0, 30, 10)], 10 @@ -149,14 +152,55 @@ def test_sum_livetime_10min(): expected_livetime = xr.DataArray( expected_values, dims=["epoch"], coords={"epoch": np.arange(30)} ) - # Call the function - result = sum_livetime_10min(livetime) + result = sum_livetime_10min(hdr_minute_cnt, livetime) + # Assert the result is as expected + xr.testing.assert_equal(result, expected_livetime) + +def test_sum_livetime_10min_incomplete_set(): + """Test the sum_livetime_10min function with an incomplete set.""" + # Create a sample livetime DataArray of 2 complete and 1 incomplete set + livetime_values = np.ones(29) # 29 epochs all ones so we can check that + # the livetime sum scaling was handled correctly. + livetime = xr.DataArray( + livetime_values, dims=["epoch"], coords={"epoch": np.arange(29)} + ) + # Create a sample hdr_minute_cnt array + hdr_minute_cnt = xr.DataArray(np.arange(10, 40), dims="epoch") + # Drop the third frame to simulate an incomplete set + hdr_minute_cnt = hdr_minute_cnt.drop_isel(epoch=3) + hdr_minute_cnt = hdr_minute_cnt.assign_coords({"epoch": np.arange(29)}) + # Expected result: sum of every 10 values (one) repeated 10 times except for + # the second set. Since the first set had a dropped frame, and the second + # set's livetime comes from the first set that livetime sum should be scaled + # correctly but still have a missing frame. + expected_values = np.repeat(np.array([10, 10, 10]), np.array([10, 9, 10])) + expected_livetime = xr.DataArray( + expected_values, dims=["epoch"], coords={"epoch": np.arange(29)} + ) + # Call the function + result = sum_livetime_10min(hdr_minute_cnt, livetime) # Assert the result is as expected xr.testing.assert_equal(result, expected_livetime) +def test_sum_livetime_10min_livetime_mismatch(): + """Test the sum_livetime_10min function when there is a shape mismatch.""" + # Create a sample livetime DataArray of 2 complete sets + livetime_values = np.ones(20) + livetime = xr.DataArray( + livetime_values, dims=["epoch"], coords={"epoch": np.arange(20)} + ) + # Create a sample hdr_minute_cnt array that conflicts with the livetime array. + hdr_minute_cnt = xr.DataArray(np.arange(10, 42)) + with pytest.raises( + ValueError, match="Livetime does not line up with the sectored counts sets." + ): + # Call the function + sum_livetime_10min(hdr_minute_cnt, livetime) + + def test_process_summed_rates_data(dependencies): """Test the variables in the summed rates dataset""" diff --git a/imap_processing/tests/hit/test_hit_l2.py b/imap_processing/tests/hit/test_hit_l2.py index facb82f5cb..ef811e1a1f 100644 --- a/imap_processing/tests/hit/test_hit_l2.py +++ b/imap_processing/tests/hit/test_hit_l2.py @@ -812,20 +812,21 @@ def test_transform_to_10_minute_chunks(): ("fe", 1), ] + # Like the L1B data, each minute only has data in its own species/energy + # slot and every other slot is filled with FILLVAL_FLOAT32. data_vars = {} slot_for = {} species_i = 0 for species, num_energy_levels in species_energy: energy_dim = f"{species}_energy_mean" - values = np.array( - [[m * 100 + e for e in range(num_energy_levels)] for m in range(n_minutes)], - dtype=np.float32, - ) - data_vars[f"{species}_macropixel_intensity"] = (("epoch", energy_dim), values) + values = np.full((n_minutes, num_energy_levels), FILLVAL_FLOAT32, np.float32) for energy_i in range(num_energy_levels): slot_for[(species, energy_i)] = species_i + minutes = np.arange(species_i, n_minutes, 10) + values[minutes, energy_i] = minutes * 100 + energy_i species_i += 1 - + data_vars[f"{species}_macropixel_intensity"] = (("epoch", energy_dim), values) + data_vars["hdr_minute_cnt"] = ("epoch", np.arange(n_minutes)) macropixel_dataset = xr.Dataset(data_vars, coords={"epoch": epochs}) result = transform_to_10_minute_chunks(macropixel_dataset) @@ -839,12 +840,15 @@ def test_transform_to_10_minute_chunks(): ) np.testing.assert_array_equal(result["epoch_delta"].values, expected_epoch_delta) - # Each new epoch is centered on its 10-minute group, then shifted back by + # Each new epoch is centered on its 10-minute set, then shifted back by # a full 10-minute chunk. for chunk in range(n_chunks): start = epochs[chunk * 10] - end = epochs[chunk * 10 + 9] - expected_epoch = start + (end - start) // 2 - SECONDS_PER_10_MIN * 1_000_000_000 + expected_epoch = ( + start + + SECONDS_PER_10_MIN * 1_000_000_000 // 2 + - SECONDS_PER_10_MIN * 1_000_000_000 + ) assert result["epoch"].values[chunk] == expected_epoch # Each species/energy variable should pull its value from the minute slot