Detect firmware event-buffer truncation and plot count rate, not raw sums

dosview/__init__.py

@@ -44,6 +44,7 @@
     OldLogParser,
     get_parser_for_file,
     parse_file,
+    TICK_SECONDS,
 )
 from .eeprom_widget import EepromManagerWidget
 from .rtc_widget import RTCManagerWidget
@@ -144,6 +145,23 @@ def __init__(self, parent=None, file_path=None):
         self._updating = False
         self._evo_region = None
         self._spec_region = None
+        self._durations = None
+        self._cur_display = np.array([])
+
+    def _extract_durations(self, data):
+        """Per-record interval length (seconds), when the source provides one.
+
+        Present for AIRDOS v2 logs/live sessions (fw 1.3+, see parsers.py);
+        absent for older logs, in which case the evolution plot falls back
+        to raw counts per record.
+        """
+        dur = data[6] if len(data) > 6 else None
+        if dur is None:
+            return None
+        dur = np.asarray(dur, dtype=float)
+        if dur.shape[0] != self._sums.shape[0]:
+            return None
+        return dur
 
     def plot(self, data):
         start_time = time.time()
@@ -163,6 +181,7 @@ def plot(self, data):
         self._sums = np.asarray(data[1], dtype=float)
         self._hist = np.asarray(data[2], dtype=float)
         self._channels = np.arange(len(self._hist), dtype=float)
+        self._durations = self._extract_durations(data)
         sm = data[5] if len(data) > 5 else None
         if sm is not None and hasattr(sm, "ndim") and sm.ndim == 2 and sm.shape[0] >= 1 and sm.shape[1] > 0:
             self._spectral_matrix = np.asarray(sm, dtype=float)
@@ -176,7 +195,15 @@ def plot(self, data):
                                    'bottom': SymLogAxisItem(orientation='bottom')})
 
         self.plot_evolution.showGrid(x=True, y=True)
-        self.plot_evolution.setLabel("left",  "Total count per exposition", units="#")
+        # Per-record intervals are no longer all the same length once a high-rate
+        # burst triggers an early firmware flush (see main.cpp flushDataOut()), so
+        # raw counts per record stop being comparable across records — plot a
+        # rate instead whenever a duration is known. Falls back to raw counts
+        # for sources with no duration info (older logs/firmware).
+        if self._durations is not None:
+            self.plot_evolution.setLabel("left", "Count rate", units="cps")
+        else:
+            self.plot_evolution.setLabel("left", "Total count per exposition", units="#")
         self.plot_evolution.setLabel("bottom", "Time", units="min")
 
         self._curve_evolution = self.plot_evolution.plot([], [],
@@ -240,8 +267,13 @@ def _on_mouse_moved(self, evt):
                 idx = int(np.argmin(np.abs(self._time_axis_min - x)))
                 t = float(self._time_axis_min[idx])
                 count = float(self._cur_sums[idx])
-                px, py = t, float(symlog(count))
-                text = f"t = {t:.2f} min\ncount = {count:.0f}"
+                value = float(self._cur_display[idx])
+                px, py = t, float(symlog(value))
+                if self._durations is not None and idx < len(self._durations):
+                    duration = float(self._durations[idx])
+                    text = f"t = {t:.2f} min\nrate = {value:.2f} cps\ncount = {count:.0f} in {duration:.2f} s"
+                else:
+                    text = f"t = {t:.2f} min\ncount = {count:.0f}"
             else:
                 if not len(self._channels):
                     continue
@@ -257,13 +289,25 @@ def _on_mouse_moved(self, evt):
                 item.setVisible(True)
 
     def _set_evolution_curve(self, sums):
-        """Draw the evolution curve (and rolling average) in symlog-y."""
+        """Draw the evolution curve (and rolling average) in symlog-y.
+
+        Plotted as counts/sec when per-record durations are known, since
+        records can now be shorter than the nominal interval (early flush on
+        a full event buffer — see flushDataOut() in firmware) and raw sums
+        of different-length records aren't directly comparable. Falls back
+        to raw counts per record otherwise.
+        """
         sums = np.asarray(sums, dtype=float)
-        self._cur_sums = sums  # currently displayed values (for the hover tooltip)
-        self._curve_evolution.setData(self._time_axis_min, symlog(sums))
+        self._cur_sums = sums  # raw per-record counts (for the hover tooltip)
+        if self._durations is not None and self._durations.shape[0] == sums.shape[0]:
+            display = sums / self._durations
+        else:
+            display = sums
+        self._cur_display = display  # currently displayed values (for the hover tooltip)
+        self._curve_evolution.setData(self._time_axis_min, symlog(display))
         window_size = self.WINDOW_SIZE
-        if len(sums) >= window_size:
-            rolling_avg = np.convolve(sums, np.ones(window_size) / window_size, mode='valid')
+        if len(display) >= window_size:
+            rolling_avg = np.convolve(display, np.ones(window_size) / window_size, mode='valid')
             self._curve_rolling_avg.setData(self._time_axis_min[window_size - 1:], symlog(rolling_avg))
         else:
             self._curve_rolling_avg.setData([], [])
@@ -393,6 +437,7 @@ def update_data(self, data):
         self._sums = np.asarray(data[1], dtype=float)
         self._hist = np.asarray(data[2], dtype=float)
         self._channels = np.arange(len(self._hist), dtype=float)
+        self._durations = self._extract_durations(data)
         sm = data[5] if len(data) > 5 else None
         if sm is not None and hasattr(sm, "ndim") and sm.ndim == 2 and sm.shape[0] >= 1 and sm.shape[1] > 0:
             self._spectral_matrix = np.asarray(sm, dtype=float)
@@ -778,6 +823,7 @@ def run(self):
         hist = np.zeros(65536, dtype=int)
         time_axis = _GrowableArray(dtype=float)
         sums = _GrowableArray(dtype=int)
+        durations = _GrowableArray(dtype=float)
         spectral_records = _GrowableMatrix(width=len(hist), dtype=int)
         metadata = {"log_runs_count": 0, "log_device_info": {}}
         fmt = None  # 'old' or 'v2'
@@ -788,6 +834,8 @@ def run(self):
 
         # v2-specific inter-record state
         current_hist = None
+        current_cap_start = 0
+        prev_abs_time = None
         current_counts = 0
 
         def _resolve_timestamp(parts):
@@ -950,6 +998,10 @@ def _build_telemetry(env_recs):
                         elif msg == "$START":
                             current_hist = np.zeros_like(hist)
                             current_counts = 0
+                            try:
+                                current_cap_start = int(parts[2]) if len(parts) > 2 else 0
+                            except ValueError:
+                                current_cap_start = 0
                         elif msg == "$E" and current_hist is not None and len(parts) >= 3:
                             try:
                                 ch = int(parts[2])
@@ -961,6 +1013,20 @@ def _build_telemetry(env_recs):
                                 logger.warning("Dropped malformed $E record: %r", line)
                         elif msg == "$STOP" and current_hist is not None:
                             try:
+                                if len(parts) > 4:
+                                    try:
+                                        ev_count = int(parts[4])
+                                    except ValueError:
+                                        ev_count = None
+                                    if ev_count is not None and ev_count > current_counts:
+                                        logger.warning(
+                                            "Device counted %d above-threshold events "
+                                            "but only %d were transmitted — the "
+                                            "firmware's per-interval event buffer was "
+                                            "exceeded, so the spectrum is undercounted "
+                                            "above THRESHOLD for this interval.",
+                                            ev_count, current_counts,
+                                        )
                                 for idx, val in enumerate(parts[5:]):
                                     try:
                                         current_hist[idx] += int(val)
@@ -973,13 +1039,38 @@ def _build_telemetry(env_recs):
                                 # fall back to a non-time field such as the run
                                 # index, which would silently compress the x-axis.
                                 t = _resolve_timestamp(parts)
+
+                                # Interval duration: _resolve_timestamp already gives a
+                                # monotonically increasing clock (device time when
+                                # GNSS-synced, host wall-clock otherwise), so a plain
+                                # delta against the previous record works for the
+                                # normal case. Only fall back to the free-running
+                                # TCNT1 ticks in captStart/stopSystime — and finally to
+                                # the old fixed cadence — for the first record, or if
+                                # a GNSS sync/unsync transition made time run backwards.
+                                if prev_abs_time is not None and t > prev_abs_time:
+                                    duration = t - prev_abs_time
+                                else:
+                                    try:
+                                        stop_ticks = int(parts[3]) if len(parts) > 3 else None
+                                    except ValueError:
+                                        stop_ticks = None
+                                    if stop_ticks is not None:
+                                        ticks_elapsed = (stop_ticks - current_cap_start) & 0xFFFF
+                                        duration = max(ticks_elapsed, 1) * TICK_SECONDS
+                                    else:
+                                        duration = 10.0
+                                prev_abs_time = t
+                                durations.append(duration)
+
                                 spectral_records.append(current_hist)
                                 hist += current_hist
                                 time_axis.append(t)
                                 sums.append(int(current_hist.sum()))
                                 self.dataUpdated.emit(
                                     [time_axis.view(), sums.view(), hist.copy(), metadata,
-                                     _build_telemetry(env_records), spectral_records.view()]
+                                     _build_telemetry(env_records), spectral_records.view(),
+                                     durations.view()]
                                 )
                             except (ValueError, IndexError):
                                 dropped_records += 1
@@ -1840,6 +1931,8 @@ def save_as(self):
                 arrays[f"telemetry_value_{key}"] = v
         if len(data) > 5 and data[5] is not None and hasattr(data[5], "shape") and data[5].ndim == 2:
             arrays["spectral_matrix"] = data[5]
+        if len(data) > 6 and data[6] is not None:
+            arrays["durations"] = data[6]
         np.savez_compressed(path, **arrays)
 
 

dosview/parsers.py

@@ -8,11 +8,40 @@
 
 from typing import List, Sequence, Tuple
 
+import logging
 import time
 from pathlib import Path
 
 import numpy as np
 
+logger = logging.getLogger(__name__)
+
+# Seconds per TCNT1 tick on the AIRDOS03 AVR firmware: Timer1 runs off the
+# 8 MHz clock with a /1024 prescaler (see fw/*/src/main.cpp, TCCR1B). Used to
+# turn the raw captStart/stopSystime tick counts in $START/$STOP into a
+# duration when the device has no GNSS fix (so tm.tm_s100 is always 0).
+TICK_SECONDS = 128e-6
+
+# TCNT1 is 16-bit, so tick-based durations are only unambiguous up to one
+# full wrap of the counter.
+MAX_TICK_DURATION_S = 65536 * TICK_SECONDS
+
+
+def _fill_invalid_durations(durations: np.ndarray) -> np.ndarray:
+    """Replace non-positive/non-finite duration entries with a sane fallback.
+
+    Individual $STOP records can fail to yield a duration (malformed line).
+    Rather than let a stray zero/NaN divide-by-zero corrupt a counts/sec
+    plot, backfill with the median of the valid durations in this log (or
+    1 s, i.e. "treat as raw counts", if none are valid at all).
+    """
+    valid = np.isfinite(durations) & (durations > 0)
+    if valid.all():
+        return durations
+    durations = durations.copy()
+    durations[~valid] = np.median(durations[valid]) if valid.any() else 1.0
+    return durations
+
 
 class BaseLogParser:
     """Base parser class."""
@@ -60,10 +89,13 @@ def parse(self):
         total_counts = 0
         sums: List[int] = []
         time_axis: List[float] = []
+        durations: List[float] = []
         spectral_records: List[np.ndarray] = []
         inside_run = False
         current_hist = None
         current_counts = 0
+        current_cap_start = 0
+        prev_abs_time: float | None = None
         device_type = "unknown"
         env_records: List[Tuple[float, ...]] = []
         batt_records: List[Tuple[float, ...]] = []
@@ -102,6 +134,10 @@ def parse(self):
                         inside_run = True
                         current_hist = np.zeros_like(hist)
                         current_counts = 0
+                        try:
+                            current_cap_start = int(parts[2]) if len(parts) > 2 else 0
+                        except ValueError:
+                            current_cap_start = 0
                     case "$E":
                         if inside_run and len(parts) >= 3:
                             channel = int(parts[2])
@@ -110,6 +146,20 @@ def parse(self):
                                 current_counts += 1
                     case "$STOP":
                         if inside_run:
+                            if len(parts) > 4:
+                                try:
+                                    ev_count = int(parts[4])
+                                except ValueError:
+                                    ev_count = None
+                                if ev_count is not None and ev_count > current_counts:
+                                    logger.warning(
+                                        "Run %d: device counted %d above-threshold "
+                                        "events but only %d were transmitted — the "
+                                        "firmware's per-interval event buffer was "
+                                        "exceeded, so the spectrum is undercounted "
+                                        "above THRESHOLD for this interval.",
+                                        metadata["log_runs_count"], ev_count, current_counts,
+                                    )
                             if len(parts) > 5:
                                 for idx, val in enumerate(parts[5:]):
                                     try:
@@ -120,7 +170,27 @@ def parse(self):
                             hist += current_hist
                             total_counts += current_counts
                             sums.append(current_counts)
-                            time_axis.append(float(parts[2]))
+                            cur_abs_time = float(parts[2])
+                            time_axis.append(cur_abs_time)
+
+                            # Interval duration: prefer the GNSS-synced wall-clock
+                            # delta between consecutive $STOP timestamps (works for
+                            # any length); fall back to the free-running TCNT1 tick
+                            # count in captStart/stopSystime (works without GNSS,
+                            # but only unambiguous up to one 16-bit timer wrap).
+                            if prev_abs_time and cur_abs_time > prev_abs_time:
+                                durations.append(cur_abs_time - prev_abs_time)
+                            else:
+                                try:
+                                    stop_ticks = int(parts[3]) if len(parts) > 3 else None
+                                except ValueError:
+                                    stop_ticks = None
+                                if stop_ticks is not None:
+                                    ticks_elapsed = (stop_ticks - current_cap_start) & 0xFFFF
+                                    durations.append(max(ticks_elapsed, 1) * TICK_SECONDS)
+                                else:
+                                    durations.append(float("nan"))
+                            prev_abs_time = cur_abs_time if cur_abs_time else None
                         inside_run = False
                         current_hist = None
                     case "$ENV":
@@ -179,8 +249,9 @@ def parse(self):
             telemetry["temperature"]        = (ba[:, 0], ba[:, 5])
 
         spectral_matrix = np.array(spectral_records) if spectral_records else np.zeros((0, hist.shape[0]), dtype=int)
+        durations_arr = _fill_invalid_durations(np.array(durations, dtype=float)) if durations else np.zeros(0)
         print("Parsed AIRDOS v2 format in", time.time() - start_time, "s")
-        return [np.array(time_axis), np.array(sums), hist, metadata, telemetry, spectral_matrix]
+        return [np.array(time_axis), np.array(sums), hist, metadata, telemetry, spectral_matrix, durations_arr]
 
 
 # Backwards-compatible alias
@@ -319,6 +390,8 @@ def parse(self):
             spectral_matrix = npz["spectral_matrix"]
         else:
             spectral_matrix = np.zeros((0, hist.shape[0]), dtype=int)
+        if "durations" in npz.files:
+            return [time_axis, sums, hist, metadata, telemetry, spectral_matrix, npz["durations"]]
         return [time_axis, sums, hist, metadata, telemetry, spectral_matrix]
 
 
@@ -345,4 +418,5 @@ def parse_file(file_path: str | Path):
     "NpzLogParser",
     "get_parser_for_file",
     "parse_file",
+    "TICK_SECONDS",
 ]