diff --git a/apps/autocal/src/orisys_autocal/multipoint.py b/apps/autocal/src/orisys_autocal/multipoint.py index d497919..9437878 100644 --- a/apps/autocal/src/orisys_autocal/multipoint.py +++ b/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 diff --git a/apps/autocal/tests/test_multipoint.py b/apps/autocal/tests/test_multipoint.py index ab85f6d..7ced388 100644 --- a/apps/autocal/tests/test_multipoint.py +++ b/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)