增加多点标定输出三个合力平均值

apps/autocal/src/orisys_autocal/multipoint.py

@@ -20,6 +20,7 @@ import csv
 import json
 import math
 import re
+import statistics
 import time
 from dataclasses import dataclass
 from datetime import date, datetime, timezone
@@ -124,6 +125,126 @@ class SiteSpec:
         return {"x": self.x, "y": self.y, "z": self.z}
 
 
+#: Columns tried, in order, for the three separate force means. The first set
+#: whose three values are all present is the one averaged. The combined
+#: magnitude is not one of them.
+_COMPONENT_SETS = (
+    ("N", ("FN_output_N", "FX_output_N", "FY_output_N")),
+    ("N", ("FN_N", "FX_N", "FY_N")),
+    ("raw", ("FN_raw", "FX_raw", "FY_raw")),
+)
+STABLE_MEANS_HEADER = (
+    "site_id",
+    "repeat_index",
+    "step_index",
+    "branch",
+    "target_load_n",
+    "samples",
+    "FN_mean",
+    "FX_mean",
+    "FY_mean",
+    "unit",
+)
+
+
+def _cell(row: dict, key: str) -> float | None:
+    text = str(row.get(key) or "").strip()
+    if not text:
+        return None
+    try:
+        value = float(text)
+    except ValueError:
+        return None
+    return value if math.isfinite(value) else None
+
+
+def _components(row: dict) -> tuple[str, tuple[float, float, float]] | None:
+    """Fn, Fx, Fy from one sample. A combined magnitude is never returned."""
+    if row.get("valid") not in (None, "", "1"):
+        return None
+    for unit, keys in _COMPONENT_SETS:
+        values = tuple(_cell(row, key) for key in keys)
+        if all(value is not None for value in values):
+            return unit, values  # type: ignore[return-value]
+    return None
+
+
+#: A point this many robust standard deviations from the median is an outlier.
+_OUTLIER_SIGMAS = 3.0
+_MAD_TO_SIGMA = 1.4826
+
+
+def _without_outliers(values: list[float]) -> list[float]:
+    """Drop values far from the median. Fewer than three samples are left as they are."""
+    if len(values) < 3:
+        return values
+    centre = statistics.median(values)
+    mad = statistics.median([abs(value - centre) for value in values])
+    if mad == 0.0:
+        kept = [value for value in values if value == centre]
+        return kept or values
+    limit = _OUTLIER_SIGMAS * _MAD_TO_SIGMA * mad
+    kept = [value for value in values if abs(value - centre) <= limit]
+    return kept or [centre]
+
+
+def stable_component_rows(repeat_dir: Path, *, site_id: str, repeat_index: int) -> list[dict]:
+    """Mean Fn, Fx and Fy over the latter half of each completed hold.
+
+    The first half of a hold is the settling transient. Each component then
+    drops its own outliers before the mean. The file this feeds does not
+    contain a combined force.
+    """
+    reference_path = repeat_dir / "reference.csv"
+    if not reference_path.is_file():
+        return []
+    with reference_path.open(encoding="utf-8-sig", newline="") as handle:
+        holds = list(csv.DictReader(handle))
+    stamps: dict[int, int] = {}
+    frames_path = repeat_dir / "capture" / "frames.csv"
+    if frames_path.is_file():
+        with frames_path.open(encoding="utf-8-sig", newline="") as handle:
+            for row in csv.DictReader(handle):
+                if row.get("dropped") not in (None, "", "0") or not row.get("t_host_mono_ns"):
+                    continue
+                stamps[int(row["frame_id"])] = int(row["t_host_mono_ns"])
+    samples: list[tuple[int, str, tuple[float, float, float]]] = []
+    forces_path = repeat_dir / "capture" / "forces.csv"
+    if forces_path.is_file():
+        with forces_path.open(encoding="utf-8-sig", newline="") as handle:
+            for row in csv.DictReader(handle):
+                components = _components(row)
+                stamp = stamps.get(int(row["frame_id"])) if row.get("frame_id") else None
+                if components is None or stamp is None:
+                    continue
+                samples.append((stamp, components[0], components[1]))
+    rows = []
+    for hold in holds:
+        if hold.get("hold_complete") not in (None, "", "1"):
+            continue
+        start = int(hold["t_start_host_mono_ns"])
+        end = int(hold["t_end_host_mono_ns"])
+        stable = start + (end - start) // 2
+        chosen = [item for item in samples if stable <= item[0] < end]
+        unit = chosen[0][1] if chosen else ""
+        used = [item[2] for item in chosen if item[1] == unit]
+        cleaned = [_without_outliers(list(axis)) for axis in zip(*used)] if used else []
+        means = tuple(sum(axis) / len(axis) for axis in cleaned) if cleaned else ("", "", "")
+        rows.append({
+            "site_id": site_id,
+            "repeat_index": repeat_index,
+            "step_index": hold.get("step_index", ""),
+            "branch": hold.get("branch", ""),
+            "target_load_n": hold.get("target_load_n", ""),
+            "samples": len(used),
+            "FN_mean": "" if means[0] == "" else f"{means[0]:.6g}",
+            "FX_mean": "" if means[1] == "" else f"{means[1]:.6g}",
+            "FY_mean": "" if means[2] == "" else f"{means[2]:.6g}",
+            "unit": unit,
+        })
+    return rows
+
+
 def load_steps(
     levels: list[float],
     max_force_n: float,
@@ -813,7 +934,25 @@ class Campaign:
         self._document["status"] = status
         self._document["reason"] = self.reason
         self._document["ended_at"] = datetime.now(timezone.utc).isoformat()
+        self._write_stable_means()
         self._save()
+
+    def _write_stable_means(self) -> None:
+        """Write Fn, Fx and Fy means. No combined force."""
+        if self.directory is None or self._document is None:
+            return
+        rows = []
+        for trial in self._document.get("trials") or []:
+            rows.extend(stable_component_rows(
+                self.directory / trial["directory"],
+                site_id=str(trial.get("site_id", "")),
+                repeat_index=int(trial.get("repeat_index") or 0),
+            ))
+        path = self.directory / "stable_means.csv"
+        with path.open("w", newline="", encoding="utf-8") as handle:
+            writer = csv.DictWriter(handle, fieldnames=list(STABLE_MEANS_HEADER))
+            writer.writeheader()
+            writer.writerows(rows)
         if self._session_events_file is not None:
             self._session_events_file.close()
             self._session_events_file = None

apps/autocal/tests/test_multipoint.py

@@ -244,3 +244,97 @@ def test_solved_pose_is_stored_and_a_point_can_be_recorded_again(tmp_path):
     session = json.loads((directory / "session.json").read_text(encoding="utf-8"))
     assert session["completion"][0]["repeats_complete"] == 2
     assert session["completion"][0]["complete"] is True
+
+
+def test_stable_means_are_the_three_components_of_the_latter_half(tmp_path):
+    """The summary is Fn, Fx, Fy after the hold settles. No combined magnitude."""
+    from orisys_autocal.multipoint import STABLE_MEANS_HEADER, stable_component_rows
+
+    repeat = tmp_path / "repeat-001"
+    (repeat / "capture").mkdir(parents=True)
+    start, end = 1_000, 1_000 + 100
+    with (repeat / "reference.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=[
+            "step_index", "branch", "target_load_n", "t_start_host_mono_ns",
+            "t_end_host_mono_ns", "hold_complete",
+        ])
+        writer.writeheader()
+        writer.writerow({
+            "step_index": 1, "branch": "load", "target_load_n": 5,
+            "t_start_host_mono_ns": start, "t_end_host_mono_ns": end, "hold_complete": 1,
+        })
+    with (repeat / "capture" / "frames.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=["frame_id", "t_host_mono_ns", "dropped"])
+        writer.writeheader()
+        for frame, stamp in ((1, start + 10), (2, start + 60), (3, start + 80)):
+            writer.writerow({"frame_id": frame, "t_host_mono_ns": stamp, "dropped": 0})
+    with (repeat / "capture" / "forces.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=[
+            "frame_id", "valid", "FN_raw", "FX_raw", "FY_raw", "Fmag_raw",
+            "FN_output_N", "FX_output_N", "FY_output_N",
+        ])
+        writer.writeheader()
+        # The first sample is in the settling half and must not enter the mean.
+        writer.writerow({
+            "frame_id": 1, "valid": 1, "FN_raw": 9, "FX_raw": 9, "FY_raw": 9, "Fmag_raw": 15.588,
+            "FN_output_N": 100, "FX_output_N": 100, "FY_output_N": 100,
+        })
+        writer.writerow({
+            "frame_id": 2, "valid": 1, "FN_raw": 1, "FX_raw": 1, "FY_raw": 1, "Fmag_raw": 1.732,
+            "FN_output_N": 2, "FX_output_N": 0.4, "FY_output_N": -0.2,
+        })
+        writer.writerow({
+            "frame_id": 3, "valid": 1, "FN_raw": 3, "FX_raw": 3, "FY_raw": 3, "Fmag_raw": 5.196,
+            "FN_output_N": 4, "FX_output_N": 0.8, "FY_output_N": -0.6,
+        })
+
+    rows = stable_component_rows(repeat, site_id="P1", repeat_index=1)
+
+    assert list(rows[0]) == list(STABLE_MEANS_HEADER)
+    assert rows[0]["samples"] == 2
+    assert rows[0]["unit"] == "N"
+    assert float(rows[0]["FN_mean"]) == pytest.approx(3.0)
+    assert float(rows[0]["FX_mean"]) == pytest.approx(0.6)
+    assert float(rows[0]["FY_mean"]) == pytest.approx(-0.4)
+    assert "Fmag" not in rows[0]
+
+
+def test_an_outlier_in_the_stable_half_is_left_out_of_its_own_mean(tmp_path):
+    from orisys_autocal.multipoint import stable_component_rows
+
+    repeat = tmp_path / "repeat-001"
+    (repeat / "capture").mkdir(parents=True)
+    start, end = 1_000, 1_100
+    with (repeat / "reference.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=[
+            "step_index", "branch", "target_load_n", "t_start_host_mono_ns",
+            "t_end_host_mono_ns", "hold_complete",
+        ])
+        writer.writeheader()
+        writer.writerow({
+            "step_index": 1, "branch": "load", "target_load_n": 5,
+            "t_start_host_mono_ns": start, "t_end_host_mono_ns": end, "hold_complete": 1,
+        })
+    stamps = (1060, 1070, 1080, 1090)
+    with (repeat / "capture" / "frames.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=["frame_id", "t_host_mono_ns", "dropped"])
+        writer.writeheader()
+        for frame, stamp in enumerate(stamps, start=1):
+            writer.writerow({"frame_id": frame, "t_host_mono_ns": stamp, "dropped": 0})
+    normals = ((2.0, 0.4, -0.2), (4.0, 0.8, -0.6), (3.0, 0.6, -0.4), (80.0, 0.5, -0.3))
+    with (repeat / "capture" / "forces.csv").open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=[
+            "frame_id", "valid", "FN_output_N", "FX_output_N", "FY_output_N",
+        ])
+        writer.writeheader()
+        for frame, (fn, fx, fy) in enumerate(normals, start=1):
+            writer.writerow({
+                "frame_id": frame, "valid": 1,
+                "FN_output_N": fn, "FX_output_N": fx, "FY_output_N": fy,
+            })
+
+    row = stable_component_rows(repeat, site_id="P1", repeat_index=1)[0]
+
+    assert float(row["FN_mean"]) == pytest.approx(3.0)
+    assert float(row["FX_mean"]) == pytest.approx(0.575)
+    assert float(row["FY_mean"]) == pytest.approx(-0.375)