Fix MAG L1C gap-fill bugs for validation tests T015 and T016#2899
Fix MAG L1C gap-fill bugs for validation tests T015 and T016#2899sapols wants to merge 10 commits intoIMAP-Science-Operations-Center:devfrom
Conversation
Six localized bug fixes to the existing MAG L1C gap-filling algorithm, plus targeted regression tests. Keeps the current production algorithm shape — no architectural rewrite. Fixes: - Rate-aware gap cadence in generate_missing_timestamps() - 2D array indexing in interpolate_gaps() (column 0 explicitly) - Off-by-one burst slice endpoint - CIC filter edge clipping via extrapolate=True in linear_filtered() - Spurious micro-gaps at Config transitions via _find_rate_segments() - Burst-only trailing-edge trim for T024 exact row count - Remove dead np.delete() no-op Refs IMAP-Science-Operations-Center#1993 Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com> Co-Authored-By: Codex <noreply@openai.com>
|
One scope clarification: My read is that this PR is enough to close the issue, but we haven't completed every historical MAG L1C TODO that came up in the discussion. Relative to @maxinelasp’s earlier checklist:
Relative to @alastairtree’s broader “L1C should” list:
|
|
Another clarification: we're probably sweeping a lot under the rug by not counting day-boundary / inherited-timeline work (T017-T019) as in-scope for this ticket. E.g. Alastair's comment about gaps in the Dec 5-6 data appears to be directly related to L1C not treating leading/trailing burst-only intervals on mixed norm+burst days as fillable gaps. In those windows the parent burst L1B data does exist, but the current mixed-input L1C path collapses the output window to the existing normal-mode coverage instead of extending it across those burst-only spans. |
|
@sapols Agreed, would you be able to open a follow-up ticket for the day boundary coverage if one does not exist? |
|
Initial thoughts from reading through the PR description only:
Nice work! I'll review the code more thoroughly in a bit. |
|
| input_filtered, vectors_filtered = cic_filter( | ||
| input_vectors, input_timestamps, output_timestamps, input_rate, output_rate | ||
| ) | ||
| return linear(vectors_filtered, input_filtered, output_timestamps) |
There was a problem hiding this comment.
Probably, we should add extrapolate=true to all the _filtered methods if we're doing it for one.
There was a problem hiding this comment.
Good catch, you're right, the same post-CIC tail-clipping bug exists in quadratic_filtered and cubic_filtered. It's dormant because L1C_INTERPOLATION_METHOD = "linear_filtered", but it'd bite us if the config changed. I'll add it to the other methods, and add an edge-case test that hits all three at the tail boundary.
|
|
||
| logger = logging.getLogger(__name__) | ||
|
|
||
| GAP_TOLERANCE = 0.075 |
There was a problem hiding this comment.
Should go in mag/constants.py. Might need a clearer name too (L1C_GAP_TOLERANCE) and a comment describing what this is.
There was a problem hiding this comment.
Done. Moved to constants.py under the new name L1C_CADENCE_TOLERANCE (to make clear that it applies to cadence matching generally, not just gap detection) and I added this comment:
# Relative tolerance (7.5%) for L1C cadence checks; allows small clock-drift
# variation around the expected cadence before a spacing is treated as a gap.
(If you like L1C_GAP_TOLERANCE better we can use it just lmk)
| 6-7 - compression flags. | ||
| """ | ||
| burst_epochs = burst_dataset["epoch"].data | ||
| norm_epochs = filled_norm_timeline[:, 0] |
There was a problem hiding this comment.
Could you rename this variable to make it clear it's the full or filled timeline?
There was a problem hiding this comment.
Renamed norm_epochs to filled_timeline_epochs
| # In the burst-only fallback, CIC delay compensation shortens the usable | ||
| # filtered range at the trailing edge by roughly one output cadence. | ||
| short_end -= int(1e9 / norm_rate.value) |
There was a problem hiding this comment.
See my comment in the main section about my concerns about this section.
There was a problem hiding this comment.
Assuming this refers to your question 3, I replied above. Happy to expand here too if useful.
| ] | ||
| gap_timeline = norm_epochs[(norm_epochs > gap[0]) & (norm_epochs < gap[1])] | ||
|
|
||
| short_end = burst_epochs[burst_end - 1] |
There was a problem hiding this comment.
The name of this variable is kind of confusing to me
There was a problem hiding this comment.
Totally fair. It was called short_end because it’s the shortened upper bound used for the short mask after the optional one-cadence trim in the burst-only fallback, but that name wasn’t clear about what was being shortened or that it was an epoch. So I renamed it to usable_burst_end_epoch.
| ) | ||
| generated_timestamps = generate_missing_timestamps(gap) | ||
| if generated_timestamps.size == 0: | ||
| last_index = int(np.searchsorted(epoch_data, gap[1], side="left")) |
| difference_ns = 0.5 * 1e9 | ||
| output: np.ndarray = np.arange(gap[0], gap[1], difference_ns) | ||
| difference_ns = int(0.5 * 1e9) | ||
| if len(gap) > 2: |
There was a problem hiding this comment.
What is the reasoning for this if statement? Add a comment plz
There was a problem hiding this comment.
The if len(gap) > 2 is there to support both legacy (start, end) gaps, which keep the historical 0.5 s cadence, and newer (start, end, rate) gaps, which use the declared cadence.
I added this comment:
# 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.
|
@maxinelasp I addressed all comments. Ready for review again! (I also updated the PR description to better explain exactly which changes fixed which tests) |
| Not required for this interpolation method. | ||
| extrapolate : bool, optional | ||
| Whether to allow extrapolation of output timestamps outside the range of input | ||
| timestamps. Default is False. |
There was a problem hiding this comment.
Can you convince/reassure me this does not lead to any extrapolated science. We are only allowed to filter/decimate science between any two MAG timestamps. those timestamps could be organic/real ones or ones where we have created a timestamp between 2 real MAG timestamps because we are downsampling/decimating the data.
Wherever we are gap filling we should usually have a period where we have data in both BM and NM streams. this allows you throw away data in the overlap period when filtering (because you need to throw away points like in the CIC filter) and still have a complete series across the gap because of the overlap.
Sometimes we change mode and do not overlap because if instrument buffer requirements. In these cases there is a short gap (up to a few secs) between streams. This should appear as a real gap and not be extrapolated.
So at least the science should never be extrapolated, the overlapping data streams allow that to be possible without gaps most of the time, and gaps should not be filled with extrapolation usually. Adding an extrapolation option makes me nervous as it sounds like making up science data which is obliviously a bad look.
There was a problem hiding this comment.
I agree with this concern - from our discussions, I believe Shawn that this will never extrapolate timestamps that end up in the final output, but maybe we could add another step in here that tracks and removes all extrapolated timestamps before returning anything, to ensure they are always cropped out. It's my understanding that in this method, the extrapolated values would never be passed out of the interpolation step - so it shouldn't affect anything to call remove_invalid_output_timestamps as a final step, right?
There was a problem hiding this comment.
@alastairtree @maxinelasp Okay, I think the requirement Alastair is describing here is stricter than I had been interpreting. So unfortunately, we need to stop here and discuss before I know how to move forward.
The current code is not extrapolating outside the original burst timestamp range, but it is technically extrapolating science beyond the CIC-filtered support a tiny amount (by less than one output cadence). For T016, the difference is exactly one L1C output row in each affected product. Without extrapolate=True, T016 is short by one row (and therefore fails); with it, that final boundary row is kept (so it passes).
That extrapolation is bounded by the CIC group delay. The filtered range is shortened at the trailing edge by decimation_factor - 1 burst-rate samples, which is (decimation_factor - 1) / decimation_factor output cadences. That is always less than one output cadence, so at most one output timestamp can land in that tail interval. In the T016 case, the extrapolated row is ~21.4 ms past the last CIC-filtered MAGO sample and ~5.8 ms past the last CIC-filtered MAGI sample.
That seems to put two expectations in conflict. The T016 validation case expects a full gapless Normal product across this Config-mode transition, including that final boundary row. But the requirement Alastair is describing says that an output row should only be produced where the CIC-filtered burst stream has real support; if the BM/NM streams do not overlap enough after filtering, the remaining boundary interval should stay a real gap.
So I think that leaves us with two clean paths:
-
If MAG considers this one sub-cadence CIC-edge row acceptable, then the current
extrapolate=Trueapproach is doing what T016 expects, and I can tighten the comments/tests around that bounded behavior. -
If MAG says absolutely no science may be produced outside the CIC-filtered support, then I should remove the
extrapolate=Truechanges from the filtered interpolation methods. That would shrink this PR and guarantee no post-CIC extrapolated samples leave interpolation, but T016 will go back to failing by one row in each affected product because the available burst data does not contain enough CIC-supported samples to produce that final boundary row.
So I need MAG direction before changing this further: should the T016 gapless-output expectation win for this one bounded edge sample, or should the no-extrapolated-science rule win even though that means T016 is one row short with the current validation inputs/expected output?
| FILLVAL = -1e31 | ||
| # Relative tolerance (7.5%) for L1C cadence checks; allows small clock-drift | ||
| # variation around the expected cadence before a spacing is treated as a gap. | ||
| L1C_CADENCE_TOLERANCE = 0.075 |
There was a problem hiding this comment.
What are the units of this? Percentage of what? Might be useful to see what this actually is for 1/2/4/8Hz products. We are pretty sensitive to timestamp jitter (appears as sin waves in our data due to s/c spin)
There was a problem hiding this comment.
This is a percentage of gap, which is defined as expected_gap = 1 / vectors_per_second * 1e9. So this attempts to compensate for very small amounts of timestamp jitter which may appear as a data gap to a more strict gap-checking system.
This is already in the existing code as a magic number. I believe we landed on that percentage when meeting in person about L1C to compensate for potential clock drift.
@sapols, can you expand the comment a little to address Alastair's question?
There was a problem hiding this comment.
And possibly rename to something like, L1C_TIMESTAMP_GAP_TOLERANCE?
There was a problem hiding this comment.
Sure thing! Renamed to L1C_TIMESTAMP_GAP_TOLERANCE and updated the comment to:
# Relative tolerance for L1C timestamp-gap checks; allows small clock-drift
# variation around the expected cadence before a spacing is treated as a gap.
# This is 7.5% of expected_gap = 1e9 / vectors_per_second ns
# (75, 37.5, 18.75, or 9.375 ms at 1, 2, 4, or 8 Hz, respectively).
L1C_TIMESTAMP_GAP_TOLERANCE = 0.075
|
Thanks for this. Have added a couple of comments inline. Feel free to send over some CDFs generated with this to the MAG team and we can take a look at the results if that helps. |
maxinelasp
left a comment
maxinelasp
left a comment
There was a problem hiding this comment.
This looks good to me, as long as we're able to address @alastairtree's concerns. I'll approve once we get Alastair's thumbs up, but no major changes requested from me :)
Thanks Shawn!!
| Not required for this interpolation method. | ||
| extrapolate : bool, optional | ||
| Whether to allow extrapolation of output timestamps outside the range of input | ||
| timestamps. Default is False. |
There was a problem hiding this comment.
I agree with this concern - from our discussions, I believe Shawn that this will never extrapolate timestamps that end up in the final output, but maybe we could add another step in here that tracks and removes all extrapolated timestamps before returning anything, to ensure they are always cropped out. It's my understanding that in this method, the extrapolated values would never be passed out of the interpolation step - so it shouldn't affect anything to call remove_invalid_output_timestamps as a final step, right?
| FILLVAL = -1e31 | ||
| # Relative tolerance (7.5%) for L1C cadence checks; allows small clock-drift | ||
| # variation around the expected cadence before a spacing is treated as a gap. | ||
| L1C_CADENCE_TOLERANCE = 0.075 |
There was a problem hiding this comment.
This is a percentage of gap, which is defined as expected_gap = 1 / vectors_per_second * 1e9. So this attempts to compensate for very small amounts of timestamp jitter which may appear as a data gap to a more strict gap-checking system.
This is already in the existing code as a magic number. I believe we landed on that percentage when meeting in person about L1C to compensate for potential clock drift.
@sapols, can you expand the comment a little to address Alastair's question?
| ) | ||
|
|
||
|
|
||
| def _find_rate_segments( |
| FILLVAL = -1e31 | ||
| # Relative tolerance (7.5%) for L1C cadence checks; allows small clock-drift | ||
| # variation around the expected cadence before a spacing is treated as a gap. | ||
| L1C_CADENCE_TOLERANCE = 0.075 |
There was a problem hiding this comment.
And possibly rename to something like, L1C_TIMESTAMP_GAP_TOLERANCE?
| ) | ||
| assert np.array_equal(output_rate_gap, expected_rate_gap) | ||
|
|
||
| # Test Case: Adjacent gaps share a real epoch boundary at gap[0]. |
There was a problem hiding this comment.
Nice, I like this test case
3 participants
Summary
Closes issue #1993
Six localized bug fixes to the existing MAG L1C gap-filling code. These fixes un-skip and pass validation tests T015 and T016 (both mago and magi) while preserving the existing passes on T013, T014, and T024. All 10 L1C validation cases now pass with exact row-count matches against the checked-in expected outputs.
This PR intentionally keeps the current production algorithm shape. It does not rewrite the gap-fill pipeline, add new dataclasses, or replace the simplified L1C flow that was developed collaboratively between SDC and the MAG team. It only patches the specific correctness bugs that were causing T015 and T016 to fail.
Background
What was failing and why
Current dev skips T015 and T016 because they produce wrong output. Running them un-skipped shows three concrete failures:
T015-mago and T024 (both sensors) already passed on dev.
The validation tests exercise increasingly complex mode transitions:
Why T015 magi failed (+5 rows): When MAG transitions through Config mode, the
vectors_per_secondmetadata timestamp doesn't perfectly align with where the cadence actually changes in the data. The oldfind_all_gaps()split the timeline at the declared metadata timestamp, andfind_gaps()misidentified the boundary misalignment as 1-sample phantom gaps, generating 5 spurious gap-fill rows. T015 mago was unaffected because its Config transition boundary happened to align cleanly. Fix 5 (_find_rate_segmentswalkback) eliminates these phantom gaps at the root.Why T016 failed (-663 mago / +330 magi rows): Three bugs interacted to produce the T016 failures. No single fix alone was sufficient:
Reverting all three together reproduces the exact dev failures (-663 mago, +330 magi), confirming they interact non-linearly.
Why T024 was already passing but needed Fix 6: T024 passed on dev because it's burst-only and doesn't exercise the broken transition paths. However, Fix 4 (
extrapolate=True) causes T024 to produce one extra trailing sample. Fix 6 (burst-only trailing-edge trim) compensates for this, preventing a regression.How this was investigated
This issue was investigated using two AI coding agents (Claude Code and Codex) in parallel, guided by manual review. The process was:
Relationship to the algorithm document
The MAG Science Algorithm Document (section 7.3.4) describes a six-step L1C process. The current code implements a simplified version of these steps that was developed collaboratively between SDC and the MAG team. This PR preserves that simplified implementation. Specifically:
Step 3 (Generate a new NM time sequence) is the area with the largest gap between the document and the code. The document describes a burst-driven approach: find the BM anchor time tC closest to the gap start tA, decimate BM timestamps, shift by (tA - tC), and trim. The current code instead generates gap timestamps arithmetically from the gap boundaries at the appropriate NM cadence, then interpolates burst data onto those timestamps.
Both approaches produce a regular NM time sequence that bridges the gap. The current approach is simpler and was explicitly approved by the MAG team. The key fix in this PR is that the cadence used for gap-fill timestamps now respects the actual NM rate declared in the gap metadata, rather than always using 0.5s (2 Hz). This aligns with the document's intent that "the rate depends on the rate of the NM data before the gap."
Step 5 (Identify remaining gaps with >1.1s threshold) is not changed by this PR. The current implementation uses a per-rate tolerance check rather than a fixed 1.1s threshold, which is arguably more correct for non-2 Hz rates.
The six bugs and their fixes
Fix 1: Rate-aware gap cadence (
generate_missing_timestamps)Bug: Gap-fill timestamps were always generated at 0.5s intervals (2 Hz) regardless of the actual normal mode rate. For non-2 Hz gaps, this was wrong by construction.
Fix: Read the rate from
gap[2](the third element of each gap tuple) and compute the cadence asint(1e9 / rate). Fall back to 0.5s if no rate is provided. Use integer arithmetic throughout to avoid float precision loss at TTJ2000-scale nanosecond timestamps.Which tests this fixed: T016 (jointly with Fixes 3 and 4). Reverting this fix alone produces T016 mago -662 rows and T016 magi +332 rows -- the dominant contributor to the T016 failure. Also covered by the new
test_generate_missing_timestamps_uses_gap_rate()and the updated rate-aware case intest_generate_timeline().Spec alignment: This directly implements the document's instruction that the interpolated rate should depend on "the rate of the NM data before the gap."
Fix 2: 2D array indexing in
interpolate_gapsBug:
filled_norm_timelineis an(n, 8)array, but the gap-timeline extraction compared the entire 2D array against scalar timestamps without specifying a column. The boolean mask broadcasted across all 8 columns instead of just the epoch column (index 0).Fix: Extract
norm_epochs = filled_norm_timeline[:, 0]once at the top of the function and filter against it explicitly.Which tests this fixed: No isolated validation delta reproduced when this line was reverted in the current patch set. This is still a real correctness fix: it makes the gap-timeline mask operate on the epoch column by construction instead of relying on non-epoch columns being numerically unrelated to TTJ2000 timestamps.
Fix 3: Off-by-one burst slice endpoint
Bug:
burst_endwas computed as an inclusive index (viamin(len-1, ...)), but then used in a Python sliceburst_epochs[burst_start:burst_end]which is end-exclusive, silently dropping the last intended burst sample from the interpolation window.Fix: Make
burst_endexclusive from the start:burst_end = min(len(burst_epochs), burst_gap_end + burst_buffer + 1).Which tests this fixed: T016. In isolated replay, reverting this fix preserved row counts but changed the first mismatched interpolated vector at the burst-resume boundary for both mago and magi.
Fix 4: CIC filter edge clipping (
linear_filtered)Bug: After the CIC filter compensates for group delay, the filtered timestamp range is shorter than the input by delay samples at the tail.
linear()then callsremove_invalid_output_timestamps(), which clips output timestamps to this shortened filtered range, dropping valid gap timestamps near the burst boundary.Fix: Pass
extrapolate=Truetolinear()insidelinear_filtered(). This matches the behavior of the spec's reference CIC implementation, which usesInterpolatedUnivariateSpline(extrapolates by default).Which tests this fixed: T016. In isolated replay, reverting this fix alone causes both T016-mago and T016-magi to lose the final interpolated row.
Interaction with T024: Adding
extrapolate=Truealone would cause a +1 row regression on T024 (burst-only). Fix 6 below compensates for this.Fix 5: Spurious micro-gaps at Config mode transitions (
_find_rate_segments)Bug: When MAG transitions between modes via Config mode, the
vectors_per_secondmetadata records a transition timestamp that does not perfectly align with where the cadence actually changes in the data. The existingfind_all_gaps()split the timeline at the declared metadata timestamps, andfind_gaps()misidentified the boundary misalignment as 1-sample phantom gaps. These caused extra rows in the T015 output.Fix: New
_find_rate_segments()function that walks each configured transition backward while the observed cadence already matches the new rate. This aligns rate segment boundaries to the actual cadence in the data, preventing spurious micro-gaps from being generated in the first place. Also extracts the existing 0.075 tolerance magic number into a sharedGAP_TOLERANCEconstant used byfind_gaps().Which tests this fixed: T015-magi (had 5 extra rows due to phantom gaps at the Config boundary). Reverting the walkback logic alone reintroduces that exact +5 row failure.
Why not a post-hoc filter: An earlier approach filtered out gaps shorter than 2 cadence periods after detection. This is simpler (around 5 lines vs around 35 lines) but cannot distinguish Config-transition artifacts from genuine single-sample data loss. The root-cause fix via
_find_rate_segments()preserves detection of legitimate short gaps while eliminating only the artifacts.Fix 6: Burst-only trailing-edge trim for T024
Bug: Fix 4 (
extrapolate=True) allows interpolation slightly past the CIC-filtered range, which is correct for T016 but causes T024 (burst-only, no normal data) to produce one extra trailing sample beyond the expected output.Fix: Detect the burst-only case by checking whether any samples in the filled timeline have
ModeFlags.NORMstatus (has_norm_context). When there is no normal-mode context, shorten the burst coverage upper bound by one output cadence to compensate for the CIC filter delay. This removes the extra trailing sample without affecting the normal+burst path.Which tests this fixed: T024-mago and T024-magi (regressed to 687 rows with Fix 4 alone, expected 686). Reverting this trim alone reproduces that exact +1 row failure.
Additional cleanup: Dead
np.deleteremovalThe existing code called
np.delete(filled_norm_timeline, timeline_index)without assigning the result. Sincenp.delete()returns a new array, this was a no-op. Removed.Defensive additions
find_gaps(): Returns empty whentimeline_datahas fewer than 2 elements, preventing crashes on single-timestamp rate segments.dtype=np.int64throughout to prevent silent float precision loss on nanosecond timestamps.generate_timeline()bug fix: Updatelast_indexwhen generated timestamps are empty, preventing duplicate epoch data copies with adjacent gaps.Test changes
pytest.skip()guard intest_mag_l1c_validation.assertcheckinglen(expected_output.index) == len(l1c["epoch"].data)before the per-row value checks. This catches off-by-one regressions immediately rather than silently passing when the output has extra rows.test_process_mag_l1cexpectations: The existing unit test now expects rate-aware cadence (17 timestamps instead of 15 in the 4 Hz gap region, with updated flag indices).test_find_all_gaps_uses_observed_transition_boundary(): Regression test for the_find_rate_segments()micro-gap fix.test_generate_missing_timestamps_uses_gap_rate(): Regression test for rate-aware cadence, including the 2 Hz fallback.test_generate_timeline(): Verifies 4 Hz gap timestamps are spaced at 0.25s.test_generate_timeline(): Verifies no duplicate timestamps when adjacent gaps share a real epoch boundary, with a mixed-rate transition.Files changed
imap_processing/mag/l1c/mag_l1c.py+131/-32numpydocimap_processing/mag/l1c/interpolation_methods.py+1/-1extrapolate=True)imap_processing/tests/mag/test_mag_l1c.py+73/-6imap_processing/tests/mag/test_mag_validation.py+1/-4poetry.lock+10/-0poetry-lockhook refreshWhat this PR does NOT change
GapFillPlandataclass orBM-anchor/shift/trimstep-3 implementationValidation results
Rerun locally on
mag-l1c-fix:Exact row-count validation (previously skipped T015/T016 now pass):