diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 8e1bafc49..673286b20 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -18,6 +18,7 @@ import re import sys from abc import ABC, abstractmethod +from datetime import datetime, timedelta from pathlib import Path from typing import final @@ -1345,16 +1346,48 @@ def do_processing( # noqa: PLR0912 ] if self.data_level == "l1c": - science_files = dependencies.get_file_paths(source="mag", data_type="l1b") + start_datetime = datetime.strptime(self.start_date, "%Y%m%d") + science_files = dependencies.get_valid_inputs_for_start_date( + start_datetime + ).get_file_paths(source="mag", data_type="l1b") input_data = [load_cdf(dep) for dep in science_files] - # Input datasets can be in any order, and are validated within mag_l1c + + # The previous day's L1C may arrive as an extra dependency (delivered by + # sds-data-manager orchestration) so today's L1C timeline can continue + # the previous day's cadence and phase across the day boundary. + previous_day_files = [ + path + for path in dependencies.get_valid_inputs_for_start_date( + start_datetime - timedelta(days=1) + ).get_file_paths(source="mag", data_type="l1c") + if self.descriptor in path.name + ] + previous_day_dataset = ( + load_cdf(previous_day_files[0]) if previous_day_files else None + ) + + # input_data is the burst mode L1B file, normal mode L1B file, or both, + # and appears in any order if len(input_data) == 1: - datasets = [mag_l1c(input_data[0], current_day)] + datasets = [ + mag_l1c( + input_data[0], + current_day, + previous_day_dataset=previous_day_dataset, + ) + ] elif len(input_data) == 2: - datasets = [mag_l1c(input_data[0], current_day, input_data[1])] + datasets = [ + mag_l1c( + input_data[0], + current_day, + input_data[1], + previous_day_dataset=previous_day_dataset, + ) + ] else: raise ValueError( - f"Invalid dependencies found for MAG L1C:" + f"Invalid current-day dependencies found for MAG L1C:" f"{dependencies}. Expected one or two dependencies." ) if self.data_level == "l1d": diff --git a/imap_processing/mag/l1c/mag_l1c.py b/imap_processing/mag/l1c/mag_l1c.py index 861a4e047..c6263b7e0 100644 --- a/imap_processing/mag/l1c/mag_l1c.py +++ b/imap_processing/mag/l1c/mag_l1c.py @@ -22,6 +22,7 @@ def mag_l1c( first_input_dataset: xr.Dataset, day_to_process: np.datetime64, second_input_dataset: xr.Dataset = None, + previous_day_dataset: xr.Dataset = None, ) -> xr.Dataset: """ Will process MAG L1C data from L1A data. @@ -40,6 +41,12 @@ def mag_l1c( The second input dataset to process. This should be burst if first_input_dataset was norm, or norm if first_input_dataset was burst. It should match the instrument - both inputs should be mago or magi. + previous_day_dataset : xr.Dataset, optional + The previous day's L1C dataset for the same sensor. When the current day + opens with a gap, timestamps generated for that gap continue the + previous day's cadence and phase so the L1C timeline stays continuous across + the day boundary. If not provided (or not usable), gaps at the start of the + day are filled with timestamps counted from the window boundary, as before. Returns ------- @@ -49,8 +56,6 @@ def mag_l1c( # TODO: # find missing sequences and output them # Missing burst file - just pass through norm file - # Missing norm file - go back to previous L1C file to find timestamps, then - # interpolate the entire day from burst input_logical_source_1 = first_input_dataset.attrs["Logical_source"] if isinstance(first_input_dataset.attrs["Logical_source"], list): @@ -63,10 +68,17 @@ def mag_l1c( first_input_dataset, second_input_dataset ) + if previous_day_dataset is not None: + previous_day_dataset = _validated_previous_day(previous_day_dataset, sensor) + interp_function = InterpolationFunction[configuration.L1C_INTERPOLATION_METHOD] if burst_mode_dataset is not None: full_interpolated_timeline: np.ndarray = process_mag_l1c( - normal_mode_dataset, burst_mode_dataset, interp_function, day_to_process + normal_mode_dataset, + burst_mode_dataset, + interp_function, + day_to_process, + previous_day_dataset=previous_day_dataset, ) elif normal_mode_dataset is not None: full_interpolated_timeline = fill_normal_data(normal_mode_dataset) @@ -272,11 +284,128 @@ def select_datasets( return normal_mode_dataset, burst_mode_dataset +def _validated_previous_day( + previous_day_dataset: xr.Dataset, sensor: str +) -> xr.Dataset | None: + """ + Validate the previous day's dataset, returning None if it is not usable. + + The previous day's dataset must be a MAG L1C dataset for the same sensor as the + current day's inputs, with at least one epoch. An unusable dataset is + ignored with a warning rather than raised, so processing still succeeds with only + one day of data. + + Parameters + ---------- + previous_day_dataset : xr.Dataset + The candidate previous day dataset. + sensor : str + The sensor of the current day's inputs, "o" (mago) or "i" (magi). + + Returns + ------- + xr.Dataset or None + The validated dataset, or None if it should be ignored. + """ + logical_source = previous_day_dataset.attrs["Logical_source"] + if isinstance(logical_source, list): + logical_source = logical_source[0] + + if "l1c" not in logical_source or logical_source[-1] != sensor: + logger.warning( + f"Ignoring previous day dataset with logical source {logical_source}; " + f"expected L1C data for sensor mag{sensor}." + ) + return None + if ( + "epoch" not in previous_day_dataset + or previous_day_dataset["epoch"].data.size == 0 + ): + logger.warning("Ignoring previous day dataset with no epochs.") + return None + return previous_day_dataset + + +def _expected_day_ns(day_to_process: np.datetime64) -> tuple[int, int]: + """ + Return the expected L1C processing window in TTJ2000 nanoseconds. + + The window is the 24-hour day extended by 30 minutes on each side. + + Parameters + ---------- + day_to_process : np.datetime64 + The day to process, in np.datetime64[D] format. + + Returns + ------- + tuple[int, int] + The (start, end) of the processing window in TTJ2000 nanoseconds. + """ + day_start = day_to_process.astype("datetime64[s]") - np.timedelta64(30, "m") + day_end = ( + day_to_process.astype("datetime64[s]") + + np.timedelta64(1, "D") + + np.timedelta64(30, "m") + ) + return ( + int(et_to_ttj2000ns(str_to_et(str(day_start)))), + int(et_to_ttj2000ns(str_to_et(str(day_end)))), + ) + + +def _get_last_timestamp_and_rate_from_previous_day_in_ns( + previous_day_dataset: xr.Dataset, midnight_ns: int +) -> tuple[int, int] | None: + """ + Get the previous day's last vector timestamp and its vector rate. + + Only samples within the previous 24-hour day count: the last timestamp is the + last one before ``midnight_ns``, and the rate is the sample spacing there, + matched against the known MAG rates. + + Parameters + ---------- + previous_day_dataset : xr.Dataset + The previous day's L1C dataset. + midnight_ns : int + Start of the current 24-hour day in TTJ2000 nanoseconds. + + Returns + ------- + tuple[int, int] or None + ``(last_timestamp_ns, rate)``, or None when the previous day has fewer than + two samples before midnight or its final spacing matches no known MAG rate. + """ + previous_epochs = previous_day_dataset["epoch"].data + last_index = int(np.searchsorted(previous_epochs, midnight_ns, side="left")) - 1 + if last_index < 1: + logger.warning( + "Previous day dataset has fewer than two samples before the current day; " + "not continuing its timeline." + ) + return None + + last_timestamp_ns = int(previous_epochs[last_index]) + spacing = float(previous_epochs[last_index] - previous_epochs[last_index - 1]) + + for vecsec in VecSec: + if _is_expected_rate(spacing, vecsec.value): + return last_timestamp_ns, vecsec.value + + logger.warning( + f"Previous day dataset ends with sample spacing {spacing} ns, which matches " + f"no known MAG rate; not continuing its timeline." + ) + return None + + def process_mag_l1c( normal_mode_dataset: xr.Dataset | None, burst_mode_dataset: xr.Dataset, interpolation_function: InterpolationFunction, day_to_process: np.datetime64 | None = None, + previous_day_dataset: xr.Dataset | None = None, ) -> np.ndarray: """ Create MAG L1C data from L1B datasets. @@ -307,6 +436,12 @@ def process_mag_l1c( The day to process, in np.datetime64[D] format. This is used to fill gaps at the beginning or end of the day if needed. If not included, these gaps will not be filled. + previous_day_dataset : xr.Dataset, optional + The previous day's L1C dataset. When the current day opens with a gap, the + timestamps generated for that gap continue the previous day's cadence + and phase instead of counting from the window boundary, keeping the timeline + continuous across the day boundary. Requires day_to_process; ignored without + it. Returns ------- @@ -315,20 +450,34 @@ def process_mag_l1c( """ day_start_ns = None day_end_ns = None + continued_gap_start_ns = None + previous_day_rate = None if day_to_process is not None: - day_start = day_to_process.astype("datetime64[s]") - np.timedelta64(30, "m") - - # get the end of the day plus 30 minutes - day_end = ( - day_to_process.astype("datetime64[s]") - + np.timedelta64(1, "D") - + np.timedelta64(30, "m") - ) - - day_start_ns = et_to_ttj2000ns(str_to_et(str(day_start))) - day_end_ns = et_to_ttj2000ns(str_to_et(str(day_end))) - + day_start_ns, day_end_ns = _expected_day_ns(day_to_process) + + previous_day_timeline = None + if previous_day_dataset is not None: + # The previous day's 24-hour day ends at the current day's midnight, + # which is the window start without its 30 minute buffer. + midnight_ns = int( + et_to_ttj2000ns(str_to_et(str(day_to_process.astype("datetime64[s]")))) + ) + previous_day_timeline = ( + _get_last_timestamp_and_rate_from_previous_day_in_ns( + previous_day_dataset, midnight_ns + ) + ) + if previous_day_timeline is not None: + last_timestamp_ns, previous_day_rate = previous_day_timeline + period_ns = int(1e9 // previous_day_rate) + # One period before the first continued timestamp at or after the window + # start: interpolate_gaps only fills points strictly inside a gap, and + # this extra leading point is removed after generate_timeline. + steps = max(1, -((last_timestamp_ns - day_start_ns) // period_ns)) + continued_gap_start_ns = last_timestamp_ns + (steps - 1) * period_ns + + inherited_gap_start_ns = None if normal_mode_dataset: norm_epoch = normal_mode_dataset["epoch"].data if "vectors_per_second" in normal_mode_dataset.attrs: @@ -339,6 +488,34 @@ def process_mag_l1c( normal_vecsec_dict = None gaps = find_all_gaps(norm_epoch, normal_vecsec_dict, day_start_ns, day_end_ns) + if ( + continued_gap_start_ns is not None + and gaps.shape[0] > 0 + and gaps[0][0] == day_start_ns + ): + logger.info( + f"MAG L1C filling the gap at the start of the day by continuing the " + f"previous day's timeline (rate {previous_day_rate} vectors/second)." + ) + inherited_gap_start_ns = continued_gap_start_ns + gaps[0] = [inherited_gap_start_ns, gaps[0][1], previous_day_rate] + elif continued_gap_start_ns is not None: + logger.info( + f"MAG L1C has no normal mode data; generating the full timeline by " + f"continuing the previous day's timeline (rate {previous_day_rate} " + f"vectors/second)." + ) + inherited_gap_start_ns = continued_gap_start_ns + norm_epoch = [inherited_gap_start_ns, day_end_ns] + gaps = np.array( + [ + [ + inherited_gap_start_ns, + day_end_ns, + previous_day_rate, + ] + ] + ) else: norm_epoch = [day_start_ns, day_end_ns] gaps = np.array( @@ -353,6 +530,10 @@ def process_mag_l1c( new_timeline = generate_timeline(norm_epoch, gaps) + if inherited_gap_start_ns is not None: + # Drop the extra leading point; see the continued gap start computation above. + new_timeline = new_timeline[new_timeline > inherited_gap_start_ns] + if normal_mode_dataset: norm_filled: np.ndarray = fill_normal_data(normal_mode_dataset, new_timeline) else: @@ -757,15 +938,23 @@ def generate_missing_timestamps(gap: np.ndarray) -> np.ndarray: ------- full_timeline: numpy.ndarray Completed timeline. + + Raises + ------ + ValueError + If the gap bounds are not integers. """ + if not np.issubdtype(np.asarray(gap).dtype, np.integer): + # float64 cannot represent TTJ2000 nanoseconds exactly. + raise ValueError(f"Gap bounds must be integer nanoseconds, got {gap}.") difference_ns = int(0.5 * 1e9) # Support both legacy (start, end) gaps, which use the historical 0.5 s cadence, # and newer (start, end, rate) gaps, which use the declared cadence. if len(gap) > 2: difference_ns = int(1e9 / int(gap[2])) output: np.ndarray = np.arange( - int(np.rint(gap[0])), - int(np.rint(gap[1])), + int(gap[0]), + int(gap[1]), difference_ns, dtype=np.int64, ) diff --git a/imap_processing/tests/mag/test_mag_l1c.py b/imap_processing/tests/mag/test_mag_l1c.py index 0f0020e35..e671f6d3d 100644 --- a/imap_processing/tests/mag/test_mag_l1c.py +++ b/imap_processing/tests/mag/test_mag_l1c.py @@ -10,6 +10,7 @@ estimate_rate, ) from imap_processing.mag.l1c.mag_l1c import ( + _expected_day_ns, fill_normal_data, find_all_gaps, find_gaps, @@ -20,7 +21,6 @@ process_mag_l1c, vectors_per_second_from_string, ) -from imap_processing.spice.time import et_to_ttj2000ns, str_to_et from imap_processing.tests.mag.conftest import ( generate_test_epoch, mag_l1a_dataset_generator, @@ -188,10 +188,12 @@ def test_process_mag_l1c(norm_dataset, burst_dataset): def test_interpolate_gaps(norm_dataset, mag_l1b_dataset): # np.array([0, 0.5, 1, 1.5, 2, 4, 4.25, 5.5, 5.75, 6]) * 1e9 - gaps = np.array([[2 * 1e9, 4 * 1e9, 2], [4.25 * 1e9, 5.5 * 1e9, 2]]) + gaps = np.array( + [[2_000_000_000, 4_000_000_000, 2], [4_250_000_000, 5_500_000_000, 2]] + ) generated_timeline = generate_timeline(norm_dataset["epoch"].data, gaps) norm_timeline: np.ndarray = fill_normal_data(norm_dataset, generated_timeline) - gaps = np.array([[2 * 1e9, 4 * 1e9, 2]]) + gaps = np.array([[2_000_000_000, 4_000_000_000, 2]]) output = interpolate_gaps( mag_l1b_dataset, gaps, norm_timeline, InterpolationFunction.linear ) @@ -271,17 +273,6 @@ def test_mag_attributes(): assert attr in output.attrs -def _ttj2000_day_bounds(day: np.datetime64) -> tuple[int, int]: - day_start = day.astype("datetime64[s]") - np.timedelta64(30, "m") - day_end = ( - day.astype("datetime64[s]") + np.timedelta64(1, "D") + np.timedelta64(30, "m") - ) - return ( - int(et_to_ttj2000ns(str_to_et(str(day_start)))), - int(et_to_ttj2000ns(str_to_et(str(day_end)))), - ) - - def _build_mag_l1b( epochs: np.ndarray, logical_source: str, vps_attr: str ) -> xr.Dataset: @@ -299,7 +290,7 @@ def _build_mag_l1b( def _build_day_aligned_mixed_l1b(day: np.datetime64) -> tuple[xr.Dataset, xr.Dataset]: - day_start_ns, _ = _ttj2000_day_bounds(day) + day_start_ns, _ = _expected_day_ns(day) nm_start = day_start_ns + 300 * 1_000_000_000 nm_end = nm_start + 120 * 1_000_000_000 @@ -326,7 +317,7 @@ def test_process_mag_l1c_leading_burst_only_coverage(): silently dropped. """ day = np.datetime64("2025-01-01") - day_start_ns, _ = _ttj2000_day_bounds(day) + day_start_ns, _ = _expected_day_ns(day) # NM starts 5 minutes into the day window and lasts 2 minutes at 2 vec/s. nm_start = day_start_ns + 300 * 1_000_000_000 @@ -371,7 +362,7 @@ def test_process_mag_l1c_leading_burst_only_coverage(): def test_process_mag_l1c_trailing_burst_only_coverage(): """Burst coverage after the NM window must be interpolated as BURST.""" day = np.datetime64("2025-01-01") - day_start_ns, day_end_ns = _ttj2000_day_bounds(day) + day_start_ns, day_end_ns = _expected_day_ns(day) # NM sits mid-day, spanning 2 minutes at 2 vec/s. nm_center = (day_start_ns + day_end_ns) // 2 @@ -410,7 +401,7 @@ def test_mag_l1c_mixed_input_uses_day_to_process(): mag_l1c() against re-introduction of the old short-circuit to None. """ day = np.datetime64("2025-01-01") - day_start_ns, _ = _ttj2000_day_bounds(day) + day_start_ns, _ = _expected_day_ns(day) nm_start = day_start_ns + 600 * 1_000_000_000 nm_end = nm_start + 60 * 1_000_000_000 @@ -437,7 +428,7 @@ def test_mag_l1c_mixed_input_uses_day_to_process(): def test_mag_l1c_burst_only_generates_l1c_output(): """Burst-only L1C generation still fills the day window from BM data.""" day = np.datetime64("2025-01-01") - day_start_ns, _ = _ttj2000_day_bounds(day) + day_start_ns, _ = _expected_day_ns(day) # Burst covers the first 2 minutes of the day window at 8 vec/s. burst_epochs = np.arange( @@ -687,7 +678,7 @@ def test_generate_timeline(): 3, [VecSec.FOUR_VECS_PER_S], gaps=[[0.5, 1], [2, 3]] ) - gaps = np.array([[0.5, 1], [2, 3]]) * 1e9 + gaps = np.array([[500_000_000, 1_000_000_000], [2_000_000_000, 3_000_000_000]]) expected_output = np.array([0, 0.25, 0.5, 1, 1.25, 1.5, 1.75, 2, 2.5, 3]) * 1e9 output = generate_timeline(epoch_test, gaps) assert np.array_equal(output, expected_output) @@ -701,7 +692,7 @@ def test_generate_timeline(): epoch_test = generate_test_epoch( 5, [VecSec.TWO_VECS_PER_S], starting_point=1, gaps=[[3, 5]] ) - gaps = np.array([[3, 5]]) * 1e9 + gaps = np.array([[3_000_000_000, 5_000_000_000]]) expected_output = np.array([1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5]) * 1e9 output = generate_timeline(epoch_test, gaps) @@ -711,7 +702,7 @@ def test_generate_timeline(): # Timeline starts at 2s but day should start at 0s epoch_beginning_gap = np.array([2, 2.5, 3, 3.5, 4]) * 1e9 # Gap from 0s to 2s (beginning of day gap) - gaps_beginning = np.array([[0, 2 * 1e9, 2]]) + gaps_beginning = np.array([[0, 2_000_000_000, 2]]) output_beginning = generate_timeline(epoch_beginning_gap, gaps_beginning) expected_beginning = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]) * 1e9 @@ -721,7 +712,7 @@ def test_generate_timeline(): # Timeline ends at 3s but day should end at 5s epoch_end_gap = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3]) * 1e9 # Gap from 3s to 5s (end of day gap) - gaps_end = np.array([[3 * 1e9, 5 * 1e9, 2]]) + gaps_end = np.array([[3_000_000_000, 5_000_000_000, 2]]) output_end = generate_timeline(epoch_end_gap, gaps_end) # Expected: 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5 @@ -731,7 +722,9 @@ def test_generate_timeline(): # Test Case: Both beginning and end of day gaps epoch_middle_only = np.array([2, 2.5, 3]) * 1e9 # Gaps at beginning (0-2s) and end (3-5s) - gaps_both_ends = np.array([[0, 2 * 1e9, 2], [3 * 1e9, 5 * 1e9, 2]]) + gaps_both_ends = np.array( + [[0, 2_000_000_000, 2], [3_000_000_000, 5_000_000_000, 2]] + ) output_both = generate_timeline(epoch_middle_only, gaps_both_ends) # Expected: 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5 @@ -741,7 +734,9 @@ def test_generate_timeline(): # Test Case: Adjacent gaps that cause sorting issues (reproduces validation bug) epoch_edge = np.array([0, 0.5, 1, 2, 3, 3.5, 4]) * 1e9 gaps = find_all_gaps(epoch_edge, vectors_per_second_from_string("0:2")) - gaps_edge = np.array([[1 * 1e9, 2 * 1e9, 2], [2 * 1e9, 3 * 1e9, 2]]) # Adjacent + gaps_edge = np.array( + [[1_000_000_000, 2_000_000_000, 2], [2_000_000_000, 3_000_000_000, 2]] + ) # Adjacent # This test case reproduces the sorting bug from the validation test # The function should work but currently fails due to sorting issue @@ -788,7 +783,9 @@ def test_generate_timeline(): # np.arange() is end-exclusive for each generated segment, and the # epoch-copy/final-append logic preserves the real boundary sample once. epoch_adj = np.array([0, 0.5, 1, 2, 3, 3.5, 4]) * 1e9 - gaps_adj = np.array([[1e9, 2e9, 2], [2e9, 3e9, 4]]) + gaps_adj = np.array( + [[1_000_000_000, 2_000_000_000, 2], [2_000_000_000, 3_000_000_000, 4]] + ) output_adj = generate_timeline(epoch_adj, gaps_adj) expected_adj = np.array([0, 0.5, 1, 1.5, 2, 2.25, 2.5, 2.75, 3, 3.5, 4]) * 1e9 assert np.array_equal(output_adj, expected_adj) @@ -1147,3 +1144,258 @@ def test_cic_filter_delay_compensation(): f"{len(input_filtered_case2)} elements, vectors_filtered has " f"{len(vectors_filtered_case2)} elements" ) + + +def _build_cross_day_datasets( + day1: np.datetime64, day2: np.datetime64 +) -> tuple[xr.Dataset, xr.Dataset, xr.Dataset, dict]: + """Build a previous-day L1C (day1) and current-day L1B (day2) scenario. + + Day 1's L1C timeline ends at 4 vec/s carrying a sub-cadence phase offset, so it + does not line up with day 2's own 2 vec/s timeline (requirement 1). + Day 2's NM data does not begin until 10 minutes into the day window, leaving a + gap at the start of the day (requirement 2). Day 2 burst covers that leading + gap so the simulated NM timestamps can be filled. + + Returns l1c_day1, norm_day2, burst_day2, and a dict of the key constants the + continuity assertions are written against. + """ + cadence_day1 = 250_000_000 # 4 vec/s + cadence_day2 = 500_000_000 # 2 vec/s + phase_ns = 73_000_000 # sub-cadence offset, < cadence_day1 + + day2_start_ns, _ = _expected_day_ns(day2) + + # Day 1's L1C ends just before day 2's window, at 4 vec/s offset by phase. + day1_last = day2_start_ns - cadence_day1 + phase_ns + day1_epochs = day1_last - np.arange(9, -1, -1) * cadence_day1 + l1c_day1 = _build_mag_l1b( + day1_epochs.astype(np.int64), "imap_mag_l1c_norm-mago", "" + ) + + # Day 2 NM starts 10 minutes into the window at 2 vec/s (phase 0 vs boundary). + day2_nm_start = day2_start_ns + 600 * 1_000_000_000 + day2_nm_epochs = np.arange( + day2_nm_start, + day2_nm_start + 120 * 1_000_000_000 + 1, + cadence_day2, + dtype=np.int64, + ) + norm_day2 = _build_mag_l1b(day2_nm_epochs, "imap_mag_l1b_norm-mago", "0:2") + + # Day 2 burst covers the leading gap (and a little past NM start) at 8 vec/s. + burst_day2_epochs = np.arange( + day2_start_ns, + day2_nm_start + 30 * 1_000_000_000, + 125_000_000, + dtype=np.int64, + ) + burst_day2 = _build_mag_l1b(burst_day2_epochs, "imap_mag_l1b_burst-mago", "0:8") + + meta = { + "day2_start_ns": day2_start_ns, + "day1_last": int(day1_last), + "cadence_day1": cadence_day1, + "cadence_day2": cadence_day2, + "phase_ns": phase_ns, + "day2_nm_start": int(day2_nm_start), + "day2_nm_epochs": day2_nm_epochs, + } + return l1c_day1, norm_day2, burst_day2, meta + + +def test_process_mag_l1c_continues_previous_day_timeline(): + """Leading-gap timestamps must continue the previous day's L1C timeline. + + When day 2 opens with a gap, the simulated + NM timeline that fills it should adopt the *previous day's* ending rate and + phase (algorithm doc 7.3.4 step 3, "regular and consistent, across boundaries + between days"), not a fresh timeline anchored at day 2's window boundary. + """ + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + l1c_day1, norm_day2, burst_day2, meta = _build_cross_day_datasets(day1, day2) + + result = process_mag_l1c( + norm_day2, + burst_day2, + InterpolationFunction.linear, + day2, + previous_day_dataset=l1c_day1, + ) + epochs_out = result[:, 0] + + # epochs in the first 10 minute gap are filled + leading = epochs_out[epochs_out < meta["day2_nm_start"]] + assert leading.size > 0 + + # The fill starts one day-1 cadence after day 1's last sample and continues at + # day 1's ending rate (4 vec/s), not day 2's 2 vec/s. Timestamps are compared + # at the ~1 us precision used across the MAG validation tests; the float64 + # epoch columns cannot hold 2025-era TTJ2000 nanoseconds exactly. + expected_leading = meta["day1_last"] + ( + np.arange(1, leading.size + 1, dtype=np.int64) * meta["cadence_day1"] + ) + assert np.allclose(leading, expected_leading, rtol=0, atol=1e3) + # Day 2's real NM samples are still used where they exist + assert np.all(np.isin(meta["day2_nm_epochs"], epochs_out)) + + +def test_mag_l1c_continues_previous_day_timeline(): + """mag_l1c() must thread the previous day's L1C timeline into the output. + + End-to-end companion to test_process_mag_l1c_continues_previous_day_timeline: the + public entry point should accept previous_day_dataset and produce an L1C epoch + axis whose start-of-day gap fill continues day 1's rate and phase, while still + using day 2's own NM samples where they exist. + """ + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + l1c_day1, norm_day2, burst_day2, meta = _build_cross_day_datasets(day1, day2) + + output = mag_l1c(norm_day2, day2, burst_day2, previous_day_dataset=l1c_day1) + epochs_out = output["epoch"].data + + leading = epochs_out[epochs_out < meta["day2_nm_start"]] + assert leading.size > 0 + + # The gap fill continues day 1's ending timeline - rate and phase - starting one + # day-1 cadence after its last sample. Compared at the ~1 us precision used + # across the MAG validation tests. + expected_leading = meta["day1_last"] + ( + np.arange(1, leading.size + 1, dtype=np.int64) * meta["cadence_day1"] + ) + assert np.allclose(leading, expected_leading, rtol=0, atol=1e3) + # Day 2's real NM samples are still used where they exist + assert np.all(np.isin(meta["day2_nm_epochs"], epochs_out)) + + +@pytest.mark.parametrize( + "logical_source", + [ + "imap_mag_l1b_norm-mago", # L1B, not L1C + "imap_mag_l1c_norm-magi", # wrong sensor + "imap_mag_l1b_burst-mago", # burst L1B + ], +) +def test_mag_l1c_ignores_unusable_previous_day(logical_source): + """An unusable previous day dataset is ignored, not raised.""" + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + l1c_day1, norm_day2, burst_day2, _ = _build_cross_day_datasets(day1, day2) + l1c_day1.attrs["Logical_source"] = logical_source + + baseline = mag_l1c(norm_day2, day2, burst_day2) + output = mag_l1c(norm_day2, day2, burst_day2, previous_day_dataset=l1c_day1) + + assert np.array_equal(output["epoch"].data, baseline["epoch"].data) + assert np.array_equal(output["vectors"].data, baseline["vectors"].data) + + +def test_mag_l1c_ignores_previous_day_with_unknown_cadence(): + """A previous day whose final spacing matches no MAG rate is ignored.""" + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + l1c_day1, norm_day2, burst_day2, meta = _build_cross_day_datasets(day1, day2) + # Rebuild day 1 with an irregular 0.7 s spacing, matching no known rate. + day1_epochs = meta["day1_last"] - np.arange(9, -1, -1) * 700_000_000 + l1c_day1 = _build_mag_l1b( + day1_epochs.astype(np.int64), "imap_mag_l1c_norm-mago", "" + ) + + baseline = mag_l1c(norm_day2, day2, burst_day2) + output = mag_l1c(norm_day2, day2, burst_day2, previous_day_dataset=l1c_day1) + + assert np.array_equal(output["epoch"].data, baseline["epoch"].data) + + +def test_mag_l1c_no_norm_continues_previous_day_timeline(): + """With no NM data, the whole output continues the previous day's timeline.""" + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + l1c_day1, _, burst_day2, meta = _build_cross_day_datasets(day1, day2) + + output = mag_l1c(burst_day2, day2, previous_day_dataset=l1c_day1) + epochs_out = output["epoch"].data + + assert epochs_out.size > 0 + # Every output timestamp continues day 1's ending timeline, one cadence after + # its last sample - not the day-aligned burst-only fallback timeline. Compared at + # the ~1 us precision used across the MAG validation tests. + assert abs(epochs_out[0] - (meta["day1_last"] + meta["cadence_day1"])) <= 1e3 + residual = (epochs_out - meta["day1_last"]) % meta["cadence_day1"] + timeline_distance = np.minimum(residual, meta["cadence_day1"] - residual) + assert np.all(timeline_distance <= 1e3) + + +def test_process_mag_l1c_previous_day_anchor_ignores_buffer_samples(): + """The continued timeline anchors on the last sample in the previous 24-hour day. + + Samples in the previous day's trailing buffer (at or past the current day's + midnight) are on a shifted timeline here; if the anchor used them, the leading fill + would carry their phase instead of the pre-midnight phase. + """ + day2 = np.datetime64("2025-01-02") + window_start_ns, _ = _expected_day_ns(day2) + midnight_ns = window_start_ns + 1800 * 1_000_000_000 + cadence = 500_000_000 + + # Pre-midnight samples on a timeline offset by 128 ns; buffer samples past midnight + # shifted by an extra quarter second. + pre_midnight = np.arange( + midnight_ns - 20 * 1_000_000_000 + 128, midnight_ns, cadence, dtype=np.int64 + ) + buffer_samples = np.arange( + midnight_ns + 128 + 250_000_000, + midnight_ns + 10 * 1_000_000_000, + cadence, + dtype=np.int64, + ) + previous_day = _build_mag_l1b( + np.concatenate([pre_midnight, buffer_samples]), + "imap_mag_l1c_norm-mago", + "", + ) + anchor = int(pre_midnight[-1]) + + nm_epochs = np.arange( + midnight_ns + 60 * 1_000_000_000, + midnight_ns + 120 * 1_000_000_000, + cadence, + dtype=np.int64, + ) + norm_day2 = _build_mag_l1b(nm_epochs, "imap_mag_l1b_norm-mago", "0:2") + burst_day2 = _build_mag_l1b( + np.arange( + midnight_ns, midnight_ns + 30 * 1_000_000_000, 125_000_000, dtype=np.int64 + ), + "imap_mag_l1b_burst-mago", + "0:8", + ) + + result = process_mag_l1c( + norm_day2, + burst_day2, + InterpolationFunction.linear, + day2, + previous_day_dataset=previous_day, + ) + epochs_out = result[:, 0] + leading = epochs_out[epochs_out < nm_epochs[0]] + + assert leading.size > 0 + assert leading[0] == anchor + cadence + assert np.all((leading - anchor) % cadence == 0) + + +def test_mag_l1c_ignores_previous_day_without_epochs(): + """A previous day dataset with no epoch variable is ignored, not raised.""" + day1 = np.datetime64("2025-01-01") + day2 = np.datetime64("2025-01-02") + _, norm_day2, burst_day2, _ = _build_cross_day_datasets(day1, day2) + no_epochs = xr.Dataset(attrs={"Logical_source": "imap_mag_l1c_norm-mago"}) + + baseline = mag_l1c(norm_day2, day2, burst_day2) + output = mag_l1c(norm_day2, day2, burst_day2, previous_day_dataset=no_epochs) + + assert np.array_equal(output["epoch"].data, baseline["epoch"].data) diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index f2027522b..a0b752e2c 100644 --- a/imap_processing/tests/test_cli.py +++ b/imap_processing/tests/test_cli.py @@ -28,6 +28,7 @@ Hit, Idex, Lo, + Mag, Spacecraft, Swe, Ultra, @@ -902,3 +903,80 @@ def test_post_processing( "naif0012.tls", "imap_sclk_0001.tsc", ] + + +@mock.patch("imap_processing.cli.check_epochs_within_day_offsets") +@mock.patch("imap_processing.cli.mag_l1c") +def test_mag_l1c_previous_day_routing( + mock_mag_l1c, mock_check_epochs, mock_instrument_dependencies +): + """A previous-day L1C dependency routes to mag_l1c's previous_day_dataset.""" + mocks = mock_instrument_dependencies + + norm_file = "imap_mag_l1b_norm-mago_20251215_v001.cdf" + burst_file = "imap_mag_l1b_burst-mago_20251215_v001.cdf" + previous_file = "imap_mag_l1c_norm-mago_20251214_v001.cdf" + input_collection = ProcessingInputCollection( + ScienceInput(norm_file), ScienceInput(burst_file), ScienceInput(previous_file) + ) + mocks["mock_pre_processing"].return_value = input_collection + + datasets_by_name = { + norm_file: xr.Dataset({}, attrs={"cdf_filename": norm_file}), + burst_file: xr.Dataset({}, attrs={"cdf_filename": burst_file}), + previous_file: xr.Dataset({}, attrs={"cdf_filename": previous_file}), + } + mocks["mock_load_cdf"].side_effect = lambda path: datasets_by_name[Path(path).name] + mock_mag_l1c.return_value = xr.Dataset( + {"epoch": ("epoch", np.array([0, 1], dtype=np.int64))}, + attrs={"cdf_filename": "l1c_out", "Logical_source": "imap_mag_l1c_norm-mago"}, + ) + mocks["mock_write_cdf"].side_effect = ["/path/to/l1c_out"] + + dependency_str = json.dumps( + [{"type": "science", "files": [norm_file, burst_file, previous_file]}] + ) + instrument = Mag( + "l1c", "norm-mago", dependency_str, "20251215", None, "v001", False + ) + instrument.process() + + assert mock_mag_l1c.call_count == 1 + call_args, call_kwargs = mock_mag_l1c.call_args + # Current-day files are the positional inputs; the previous day's L1C file is + # passed separately and never treated as a current-day input. + assert call_args[0] is datasets_by_name[norm_file] + assert call_args[2] is datasets_by_name[burst_file] + assert call_kwargs["previous_day_dataset"] is datasets_by_name[previous_file] + + +@mock.patch("imap_processing.cli.check_epochs_within_day_offsets") +@mock.patch("imap_processing.cli.mag_l1c") +def test_mag_l1c_without_previous_day( + mock_mag_l1c, mock_check_epochs, mock_instrument_dependencies +): + """With only current-day dependencies, previous_day_dataset is None.""" + mocks = mock_instrument_dependencies + + norm_file = "imap_mag_l1b_norm-mago_20251215_v001.cdf" + input_collection = ProcessingInputCollection(ScienceInput(norm_file)) + mocks["mock_pre_processing"].return_value = input_collection + + norm_dataset = xr.Dataset({}, attrs={"cdf_filename": norm_file}) + mocks["mock_load_cdf"].side_effect = lambda path: norm_dataset + mock_mag_l1c.return_value = xr.Dataset( + {"epoch": ("epoch", np.array([0, 1], dtype=np.int64))}, + attrs={"cdf_filename": "l1c_out", "Logical_source": "imap_mag_l1c_norm-mago"}, + ) + mocks["mock_write_cdf"].side_effect = ["/path/to/l1c_out"] + + dependency_str = json.dumps([{"type": "science", "files": [norm_file]}]) + instrument = Mag( + "l1c", "norm-mago", dependency_str, "20251215", None, "v001", False + ) + instrument.process() + + assert mock_mag_l1c.call_count == 1 + call_args, call_kwargs = mock_mag_l1c.call_args + assert call_args[0] is norm_dataset + assert call_kwargs["previous_day_dataset"] is None