Skip to content

Minimally fix MAG L1C validation failures on the current gap-fill path #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
sapols wants to merge 1 commit into
base: dev
Choose a base branch
from
Draft

Minimally fix MAG L1C validation failures on the current gap-fill path #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion imap_processing/mag/l1c/interpolation_methods.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,7 +296,12 @@ def linear_filtered(
input_filtered, vectors_filtered = cic_filter(
input_vectors, input_timestamps, output_timestamps, input_rate, output_rate
)
return linear(vectors_filtered, input_filtered, output_timestamps)
return linear(
vectors_filtered,
input_filtered,
output_timestamps,
extrapolate=True,
)


def quadratic_filtered(
Expand Down
129 changes: 96 additions & 33 deletions imap_processing/mag/l1c/mag_l1c.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@

logger = logging.getLogger(__name__)

GAP_TOLERANCE = 0.075


def mag_l1c(
first_input_dataset: xr.Dataset,
Expand Down Expand Up @@ -461,6 +463,7 @@ def interpolate_gaps(
6-7 - compression flags.
"""
burst_epochs = burst_dataset["epoch"].data
norm_epochs = filled_norm_timeline[:, 0]
# Exclude range values
burst_vectors = burst_dataset["vectors"].data
# Default to two vectors per second
Expand Down Expand Up @@ -497,14 +500,12 @@ def interpolate_gaps(
burst_buffer = int(required_seconds * burst_rate.value)

burst_start = max(0, burst_gap_start - burst_buffer)
burst_end = min(len(burst_epochs) - 1, burst_gap_end + burst_buffer)
burst_end = min(len(burst_epochs), burst_gap_end + burst_buffer + 1)

gap_timeline = filled_norm_timeline[
(filled_norm_timeline > gap[0]) & (filled_norm_timeline < gap[1])
]
gap_timeline = norm_epochs[(norm_epochs > gap[0]) & (norm_epochs < gap[1])]

short = (gap_timeline >= burst_epochs[burst_start]) & (
gap_timeline <= burst_epochs[burst_end]
gap_timeline <= burst_epochs[burst_end - 1]
)
num_short = int(short.sum())

Expand All @@ -513,6 +514,8 @@ def interpolate_gaps(

# Limit timestamps to only include the areas with burst data
gap_timeline = gap_timeline[short]
if gap_timeline.size == 0:
continue
# do not include range
adjusted_gap_timeline, gap_fill = interpolation_function(
burst_vectors[burst_start:burst_end, :3],
Expand Down Expand Up @@ -542,13 +545,12 @@ def interpolate_gaps(
missing_timeline = np.setdiff1d(gap_timeline, adjusted_gap_timeline)

for timestamp in missing_timeline:
timeline_index = np.searchsorted(filled_norm_timeline[:, 0], timestamp)
timeline_index = np.searchsorted(norm_epochs, timestamp)
if filled_norm_timeline[timeline_index, 5] != ModeFlags.MISSING.value:
raise RuntimeError(
"Self-inconsistent data. "
"Gaps not included in final timeline should be missing."
)
np.delete(filled_norm_timeline, timeline_index)

return filled_norm_timeline

Expand All @@ -557,8 +559,8 @@ def generate_timeline(epoch_data: np.ndarray, gaps: np.ndarray) -> np.ndarray:
"""
Generate a new timeline from existing, gap-filled timeline and gaps.

The gaps are generated at a .5 second cadence, regardless of the cadence of the
existing data.
The gaps are generated at the cadence implied by the gap rate. If no rate is
provided, a default cadence of 0.5 seconds is used.

Parameters
----------
Expand All @@ -573,15 +575,24 @@ def generate_timeline(epoch_data: np.ndarray, gaps: np.ndarray) -> np.ndarray:
numpy.ndarray
The new timeline, filled with the existing data and the generated gaps.
"""
full_timeline: np.ndarray = np.array([])
epoch_data = np.asarray(epoch_data)
full_timeline: np.ndarray = np.array([], dtype=epoch_data.dtype)
last_index = 0
for gap in gaps:
epoch_start_index = np.searchsorted(epoch_data, gap[0], side="left")
gap_start_in_epoch = (
epoch_start_index < epoch_data.shape[0]
and epoch_data[epoch_start_index] == int(np.rint(gap[0]))
)
epoch_copy_end = epoch_start_index + int(gap_start_in_epoch)
full_timeline = np.concatenate(
(full_timeline, epoch_data[last_index:epoch_start_index])
(full_timeline, epoch_data[last_index:epoch_copy_end])
)
generated_timestamps = generate_missing_timestamps(gap)
if gap_start_in_epoch and generated_timestamps.size != 0:
generated_timestamps = generated_timestamps[1:]
if generated_timestamps.size == 0:
last_index = int(np.searchsorted(epoch_data, gap[1], side="left"))
continue

# Remove any generated timestamps that are already in the timeline
Expand Down Expand Up @@ -639,37 +650,40 @@ def find_all_gaps(
specified as (start, end, vector_rate) where start and end both exist in the
timeline.
"""
gaps: np.ndarray = np.zeros((0, 3))
gaps: np.ndarray = np.empty((0, 3), dtype=np.int64)

# TODO: when we go back to the previous file, also retrieve expected
# vectors per second

vecsec_dict = {0: VecSec.TWO_VECS_PER_S.value} | (vecsec_dict or {})

end_index = epoch_data.shape[0]
rate_segments = _find_rate_segments(epoch_data, vecsec_dict)
if rate_segments:
first_rate = rate_segments[0][1]
last_rate = rate_segments[-1][1]
else:
default_rate = next(iter(vecsec_dict.values()))
first_rate = default_rate
last_rate = default_rate

if start_of_day_ns is not None and epoch_data[0] > start_of_day_ns:
# Add a gap from the start of the day to the first timestamp
gaps = np.concatenate(
(gaps, np.array([[start_of_day_ns, epoch_data[0], vecsec_dict[0]]]))
(gaps, np.array([[start_of_day_ns, epoch_data[0], first_rate]], dtype=np.int64))
)

for start_time in reversed(sorted(vecsec_dict.keys())):
# Find the start index that is equal to or immediately after start_time
start_index = np.searchsorted(epoch_data, start_time, side="left")
gaps = np.concatenate(
(
find_gaps(
epoch_data[start_index : end_index + 1], vecsec_dict[start_time]
),
gaps,
)
for index, (start_index, vectors_per_second) in enumerate(rate_segments):
next_start_index = (
rate_segments[index + 1][0]
if index + 1 < len(rate_segments)
else epoch_data.shape[0] - 1
)
end_index = start_index
epoch_slice = epoch_data[start_index : next_start_index + 1]
gaps = np.concatenate((gaps, find_gaps(epoch_slice, vectors_per_second)))

if end_of_day_ns is not None and epoch_data[-1] < end_of_day_ns:
gaps = np.concatenate(
(gaps, np.array([[epoch_data[-1], end_of_day_ns, vecsec_dict[start_time]]]))
(gaps, np.array([[epoch_data[-1], end_of_day_ns, last_rate]], dtype=np.int64))
)

return gaps
Expand All @@ -696,14 +710,19 @@ def find_gaps(timeline_data: np.ndarray, vectors_per_second: int) -> np.ndarray:
end_gap, as well as vectors_per_second. Start_gap and end_gap both correspond
to points in timeline_data.
"""
if timeline_data.shape[0] < 2:
return np.empty((0, 3), dtype=np.int64)

# Expected difference between timestamps in nanoseconds.
expected_gap = 1 / vectors_per_second * 1e9

diffs = abs(np.diff(timeline_data))

# Gap can be up to 7.5% larger than expected vectors per second due to clock drift
gap_index = np.asarray(diffs - expected_gap > expected_gap * 0.075).nonzero()[0]
output: np.ndarray = np.zeros((len(gap_index), 3))
gap_index = np.asarray(diffs - expected_gap > expected_gap * GAP_TOLERANCE).nonzero()[
0
]
output: np.ndarray = np.zeros((len(gap_index), 3), dtype=np.int64)

for index, gap in enumerate(gap_index):
output[index, :] = [
Expand All @@ -719,8 +738,8 @@ def generate_missing_timestamps(gap: np.ndarray) -> np.ndarray:
"""
Generate a new timeline from input gaps.

Any gaps specified in gaps will be filled with timestamps that are 0.5 seconds
apart.
Any gaps specified in gaps will be filled with timestamps at the gap rate. If the
gap rate is not included, the default cadence is 0.5 seconds.

Parameters
----------
Expand All @@ -734,12 +753,56 @@ def generate_missing_timestamps(gap: np.ndarray) -> np.ndarray:
full_timeline: numpy.ndarray
Completed timeline.
"""
# Generated timestamps should always be 0.5 seconds apart
difference_ns = 0.5 * 1e9
output: np.ndarray = np.arange(gap[0], gap[1], difference_ns)
difference_ns = int(0.5 * 1e9)
if gap.shape[0] > 2:
difference_ns = int(1e9 / int(gap[2]))
output: np.ndarray = np.arange(
int(np.rint(gap[0])),
int(np.rint(gap[1])),
difference_ns,
dtype=np.int64,
)
return output


def _is_expected_rate(timestamp_difference: float, vectors_per_second: int) -> bool:
"""Return True when a timestamp spacing matches a rate within tolerance."""
expected_gap = 1 / vectors_per_second * 1e9
return abs(timestamp_difference - expected_gap) <= expected_gap * GAP_TOLERANCE


def _find_rate_segments(
epoch_data: np.ndarray, vecsec_dict: dict[int, int]
) -> list[tuple[int, int]]:
"""
Build contiguous rate segments using observed cadence near each transition.

Walk each configured transition backward while the observed cadence already matches
the new rate so gaps stay attached to the correct segment instead of producing
spurious single-sample micro-gaps at delayed Config boundaries.
"""
if epoch_data.shape[0] == 0:
return []

segments: list[tuple[int, int]] = []
for start_time, vectors_per_second in sorted(vecsec_dict.items()):
start_index = int(np.searchsorted(epoch_data, start_time, side="left"))
start_index = min(start_index, epoch_data.shape[0] - 1)
lower_bound = segments[-1][0] if segments else 0

while start_index > lower_bound and _is_expected_rate(
epoch_data[start_index] - epoch_data[start_index - 1], vectors_per_second
):
start_index -= 1

if segments and start_index == segments[-1][0]:
segments[-1] = (start_index, vectors_per_second)
else:
segments.append((start_index, vectors_per_second))

return segments


def vectors_per_second_from_string(vecsec_string: str) -> dict:
"""
Extract the vectors per second from a string into a dictionary.
Expand Down
86 changes: 80 additions & 6 deletions imap_processing/tests/mag/test_mag_l1c.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
fill_normal_data,
find_all_gaps,
find_gaps,
generate_missing_timestamps,
generate_timeline,
interpolate_gaps,
mag_l1c,
Expand Down Expand Up @@ -112,7 +113,27 @@ def test_interpolation_methods():
def test_process_mag_l1c(norm_dataset, burst_dataset):
l1c = process_mag_l1c(norm_dataset, burst_dataset, InterpolationFunction.linear)
expected_output_timeline = (
np.array([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.25, 4.75, 5.25, 5.5, 5.75, 6])
np.array(
[
0,
0.5,
1,
1.5,
2,
2.5,
3,
3.5,
4,
4.25,
4.5,
4.75,
5,
5.25,
5.5,
5.75,
6,
]
)
* 1e9
)
assert np.array_equal(l1c[:, 0], expected_output_timeline)
Expand All @@ -122,12 +143,12 @@ def test_process_mag_l1c(norm_dataset, burst_dataset):
np.count_nonzero([np.sum(l1c[i, 1:4]) for i in range(l1c.shape[0])])
== l1c.shape[0] - 1
)
expected_flags = np.zeros(15)
expected_flags = np.zeros(17)
# filled sections should have 1 as a flag
expected_flags[5:8] = 1
expected_flags[10:11] = 1
expected_flags[10:13] = 1
# last datapoint in the gap is missing a value
expected_flags[11] = -1
expected_flags[13] = -1
assert np.array_equal(l1c[:, 5], expected_flags)
assert np.array_equal(l1c[:5, 1:5], norm_dataset["vectors"].data[:5, :])
for i in range(5, 8):
Expand All @@ -140,7 +161,7 @@ def test_process_mag_l1c(norm_dataset, burst_dataset):
assert np.allclose(l1c[i, 1:5], burst_vectors, rtol=0, atol=1)

assert np.array_equal(l1c[8:10, 1:5], norm_dataset["vectors"].data[5:7, :])
for i in range(10, 11):
for i in range(10, 13):
e = l1c[i, 0]
burst_vectors = burst_dataset.sel(epoch=int(e), method="nearest")[
"vectors"
Expand All @@ -149,7 +170,8 @@ def test_process_mag_l1c(norm_dataset, burst_dataset):
# identical.
assert np.allclose(l1c[i, 1:5], burst_vectors, rtol=0, atol=1)

assert np.array_equal(l1c[11, 1:5], [0, 0, 0, 0])
assert np.array_equal(l1c[13, 1:5], [0, 0, 0, 0])
assert np.array_equal(l1c[14:, 1:5], norm_dataset["vectors"].data[7:, :])


def test_interpolate_gaps(norm_dataset, mag_l1b_dataset):
Expand Down Expand Up @@ -435,6 +457,26 @@ def test_find_gaps():
assert np.array_equal(gaps, expected_return)


def test_find_all_gaps_uses_observed_transition_boundary():
epoch_transition = np.array([0, 0.5, 1, 1.5, 2, 10, 11, 12, 13, 14, 15, 16]) * 1e9
vectors_per_second_transition = vectors_per_second_from_string("0:2,15000000000:1")
output_transition = find_all_gaps(epoch_transition, vectors_per_second_transition)
expected_transition = np.array([[2 * 1e9, 10 * 1e9, 2]])
assert np.array_equal(output_transition, expected_transition)


def test_generate_missing_timestamps_uses_gap_rate():
gap = np.array([1_000_000_000, 2_000_000_000, 4], dtype=np.int64)
expected_output = np.array(
[1_000_000_000, 1_250_000_000, 1_500_000_000, 1_750_000_000], dtype=np.int64
)
assert np.array_equal(generate_missing_timestamps(gap), expected_output)

legacy_gap = np.array([1_000_000_000, 2_000_000_000], dtype=np.int64)
legacy_expected = np.array([1_000_000_000, 1_500_000_000], dtype=np.int64)
assert np.array_equal(generate_missing_timestamps(legacy_gap), legacy_expected)


def test_generate_timeline():
epoch_test = generate_test_epoch(
3, [VecSec.FOUR_VECS_PER_S], gaps=[[0.5, 1], [2, 3]]
Expand Down Expand Up @@ -504,6 +546,38 @@ def test_generate_timeline():
expected_edge = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]) * 1e9
assert np.array_equal(output_edge, expected_edge)

# Test Case: Gap fill uses the gap rate instead of a fixed 0.5 second cadence
epoch_rate_gap = np.array([0, 0.25, 0.5, 0.75, 1, 4, 4.25, 4.5]) * 1e9
gaps_rate_gap = np.array([[1_000_000_000, 4_000_000_000, 4]])
output_rate_gap = generate_timeline(epoch_rate_gap, gaps_rate_gap)
expected_rate_gap = (
np.array(
[
0,
0.25,
0.5,
0.75,
1,
1.25,
1.5,
1.75,
2,
2.25,
2.5,
2.75,
3,
3.25,
3.5,
3.75,
4,
4.25,
4.5,
]
)
* 1e9
)
assert np.array_equal(output_rate_gap, expected_rate_gap)


def test_gap_detection_timeline_generation_workflow():
# Create a test dataset with gaps
Expand Down
Loading