experiments: signals lab round 1, which extra phone signals sharpen clearances

experiments/drift-anatomy/METHODS.md

@@ -0,0 +1,32 @@
+# Drift anatomy: methods and full results
+
+The README is the front door. This file keeps the method as it was pre-registered in 397507f, and the full result as the run reported it.
+
+## Question
+
+The server allows 0.16 ft of error per foot along the wall. That allowance dominates a check's error bar beyond a few feet, so which part of ARKit's error is worth attacking? The data is MARViN (iPhone 14 Pro Max, bar and church, walks within 10% of the reference's scale) and ADVIO 20–22 (iPhone 6s), using the evals' windows.
+
+- **A1, return to the meter.** When a walk comes back within 0.3 m of an earlier point more than 30 s later, ARKit's closure residual r is the error in the displacement between the two visits. Does r bound the largest error between the visits?
+- **A2, scale.** Is ARKit's scale one number per walk? This measures the spread of scale over sliding 5 m windows within a walk.
+- **A3, the wall plane.** What share of the error runs across the direction of travel, which a wall plane would remove, compared with along it?
+- **A4, a range to the meter.** A second phone left at the meter ranges the walking phone by UWB (σ 10 cm, simulated) at every window end. What is the p90 left over?
+
+## Pass criteria, fixed before the run
+
+- **A1.** Across loops, the largest error between visits is at most max(2r, 2 in) in at least 90% of loops, so r can set a per-session error bar. Needs at least 8 loops.
+- **A2.** One reference object near the meter is enough if the within-walk spread is under 1% (1 SD). Drop that idea if it is over 3%.
+- **A3.** A wall-plane constraint is worth building if the along-travel p90 is at most 60% of the total p90 at 20 ft.
+- **A4.** The p90 at 20 and 30 ft is at most 6 in.
+
+## Result
+
+All four fail. Tables: [results/drift_anatomy.md](results/drift_anatomy.md).
+
+- **A1.** 30 of 56 MARViN loops (54%). ADVIO has no revisit under 0.3 m while ARKit tracks. Scale and heading error cancel on return, so r misses them.
+- **A2.** 5.7% SD on MARViN, 19% on ADVIO: drop. A few outlier windows inflate MARViN's SD; its robust SD is 1.8%.
+- **A3.** At 20 ft, along is 0.92 of the p90 on MARViN and 0.76 on ADVIO.
+- **A4.** MARViN goes from 10.4 and 13.6 in to 7.7 and 8.5 in. At σ 10 cm the range's own p90 is 6.5 in, so the bar was nearly out of reach. At 5 cm: 4.9 and 6.1 in.
+
+**Run:** `uv run python run.py`. Needs the evals checkout at 190ed33 (`EVALS_HARNESS`) and MARViN and ADVIO in `~/house-scanning-data`.
+
+**What it changes:** the error worth attacking runs along travel. A range to the meter reduces it, a wall plane wouldn't, and r can't bound a session.

experiments/drift-anatomy/README.md

@@ -0,0 +1,25 @@
+# Drift anatomy
+
+## Question
+
+The server allows 0.16 ft of error per foot along the wall, which dominates a check's error bar beyond a few feet. Which part of ARKit's error is worth attacking? A1 asks whether returning to the meter bounds it, A2 whether scale is one number per walk, A3 how much runs across travel, where a wall plane would help, and A4 what a UWB range to the meter leaves. The data is MARViN (iPhone 14 Pro Max) and ADVIO (iPhone 6s); [METHODS.md](METHODS.md) has the details.
+
+## Pass criteria, fixed before the run (397507f)
+
+- **A1.** Across loops, the largest error between visits is at most max(2r, 2 in) in at least 90% of loops, so r can set a per-session error bar. Needs at least 8 loops.
+- **A2.** One reference object near the meter is enough if the within-walk spread is under 1% (1 SD). Drop that idea if it is over 3%.
+- **A3.** A wall-plane constraint is worth building if the along-travel p90 is at most 60% of the total p90 at 20 ft.
+- **A4.** The p90 at 20 and 30 ft is at most 6 in.
+
+## Result
+
+All four fail ([tables](results/drift_anatomy.md)).
+
+- **A1.** The closure residual bounds the peak error in 30 of 56 loops. Scale and heading error cancel on return.
+- **A2.** Scale spread within a walk is 1.8% robust SD (5.7% SD) on MARViN.
+- **A3.** At 20 ft, the along-travel p90 is 9.6 of 10.4 in.
+- **A4.** UWB cuts 20 ft from 10.4 to 7.7 in at σ 10 cm, and to 4.9 in at σ 5 cm.
+
+**Run:** `uv run python run.py`: 20 s. Needs the evals checkout at 190ed33 and MARViN and ADVIO in `~/house-scanning-data`.
+
+**What it changes:** the error worth attacking runs along travel. Measuring a span inside one photo avoids it. A wall plane, a return to the meter or one ruler at the meter don't.

experiments/drift-anatomy/metrics.py

@@ -0,0 +1,185 @@
+"""Geometry for the drift-anatomy questions. Pure numpy, no dataset access; every function here has
+a strict test on a synthetic trajectory in tests/test_metrics.py.
+
+Conventions: positions are (N, 3) in a y-up frame, meters. "Horizontal" is the (x, z) part. A
+heading offset delta turns a y-up vector about +y the way `evals.modern_arkit.turn` and
+`evals.geometry.rot_y` do, so ARKit displacements turned by the offset at a window's start are in the
+reference's frame.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+MAD_TO_SD = 1.4826  # median absolute deviation to standard deviation, for normal data
+
+
+def turn(v: np.ndarray, delta: np.ndarray | float) -> np.ndarray:
+    """Turn y-up vectors (N, 3) about +y by delta (radians); the same map as evals.modern_arkit.turn."""
+    c, s = np.cos(delta), np.sin(delta)
+    return np.c_[v[:, 0] * c + v[:, 2] * s, v[:, 1], v[:, 2] * c - v[:, 0] * s]
+
+
+def yaw_offset(R_from: np.ndarray, R_to: np.ndarray) -> np.ndarray:
+    """Heading angle (N,) of the rotation about +y closest to R_to @ R_from^T, radians. Turning a
+    vector by it with `turn` equals applying evals.drift.yaw_between(R_from, R_to)."""
+    M = R_to @ np.swapaxes(R_from, -1, -2)
+    return np.arctan2(M[..., 0, 2] - M[..., 2, 0], M[..., 0, 0] + M[..., 2, 2])
+
+
+def at(p: np.ndarray, j: np.ndarray) -> np.ndarray:
+    """Positions at fractional sample indices, linearly interpolated (as evals.drift.at)."""
+    lo = np.floor(j).astype(int)
+    hi = np.minimum(lo + 1, len(p) - 1)
+    f = (j - lo)[:, None]
+    return p[lo] * (1 - f) + p[hi] * f
+
+
+def path_length(p: np.ndarray) -> np.ndarray:
+    """Cumulative horizontal distance travelled, starting at 0."""
+    steps = np.linalg.norm(np.diff(p[:, [0, 2]], axis=0), axis=1)
+    return np.concatenate([[0.0], np.cumsum(steps)])
+
+
+def horizontal_norm(v: np.ndarray) -> np.ndarray:
+    return np.linalg.norm(v[..., [0, 2]], axis=-1)
+
+
+# A1: return to the meter
+
+
+def find_loops(
+    times: np.ndarray,
+    ref: np.ndarray,
+    min_gap_s: float = 30.0,
+    max_dist_m: float = 0.3,
+    separation_s: float = 5.0,
+) -> list[tuple[int, int]]:
+    """Revisits (i, j): times[j] - times[i] > min_gap_s and |ref[j] - ref[i]| < max_dist_m.
+
+    Each start i keeps its closest qualifying j. Starts less than separation_s apart form one
+    cluster (so their separation_s-wide neighbourhoods overlap), and each cluster keeps only its
+    closest pair. A walker who pauses at the meter, or retraces a stretch of path, therefore
+    yields one loop, not one per sample. Returned in time order.
+    """
+    if len(times) != len(ref):
+        raise ValueError(f"{len(times)} times for {len(ref)} positions")
+    if np.any(np.diff(times) <= 0):
+        raise ValueError("times must increase")
+    candidates = []
+    for i in range(len(times)):
+        later = np.flatnonzero(times - times[i] > min_gap_s)
+        if len(later) == 0:
+            break
+        d = np.linalg.norm(ref[later] - ref[i], axis=1)
+        k = int(np.argmin(d))
+        if d[k] < max_dist_m:
+            candidates.append((float(d[k]), i, int(later[k])))
+    loops: list[tuple[int, int]] = []
+    cluster: list[tuple[float, int, int]] = []
+    for c in candidates:  # already in time order of i
+        if cluster and times[c[1]] - times[cluster[-1][1]] >= separation_s:
+            loops.append(min(cluster)[1:])
+            cluster = []
+        cluster.append(c)
+    if cluster:
+        loops.append(min(cluster)[1:])
+    return loops
+
+
+def loop_errors(
+    ark: np.ndarray, ref: np.ndarray, delta: np.ndarray, i: int, j: int
+) -> dict[str, float]:
+    """For a loop from sample i to sample j: the horizontal error of ARKit's displacement from i
+    (turned by the heading offset at i) against the reference's, at j (the residual r) and at its
+    largest over the samples strictly between (the peak). Also the reference's horizontal path
+    length from i to j and its largest horizontal distance from the point at i (the excursion)."""
+    if not 0 <= i < j - 1 < len(ref) - 1:
+        raise ValueError(f"need a sample strictly between i={i} and j={j}")
+    seg = slice(i, j + 1)
+    da = turn(ark[seg] - ark[i], delta[i])
+    dr = ref[seg] - ref[i]
+    e = horizontal_norm(da - dr)
+    return {
+        "r": float(e[-1]),
+        "peak": float(e[1:-1].max()),
+        "path_m": float(path_length(ref[seg])[-1]),
+        "excursion_m": float(horizontal_norm(dr).max()),
+    }
+
+
+# A2: local scale
+
+
+def path_windows(ref: np.ndarray, length_m: float, every_m: float) -> tuple[np.ndarray, np.ndarray]:
+    """Windows of exactly `length_m` of horizontal reference path. Starts are the first sample at
+    or past each multiple of `every_m` along the path (duplicates removed); ends are fractional
+    sample indices, interpolated to the exact length. Windows that run off the end are dropped."""
+    walked = path_length(ref)
+    marks = np.arange(0.0, walked[-1] - length_m, every_m)
+    i = np.unique(np.searchsorted(walked, marks))
+    target = walked[i] + length_m
+    j0 = np.searchsorted(walked, target)
+    keep = j0 < len(ref)
+    i, j0, target = i[keep], j0[keep], target[keep]
+    frac = (target - walked[j0 - 1]) / (walked[j0] - walked[j0 - 1])
+    return i, (j0 - 1) + frac
+
+
+def local_scale(ark: np.ndarray, ref: np.ndarray, i: np.ndarray, j: np.ndarray) -> np.ndarray:
+    """ARKit's straight-line displacement over the reference's, per window."""
+    return np.linalg.norm(at(ark, j) - ark[i], axis=1) / np.linalg.norm(at(ref, j) - ref[i], axis=1)
+
+
+def robust_sd(x: np.ndarray) -> float:
+    """1.4826 x the median absolute deviation from the median."""
+    return float(MAD_TO_SD * np.median(np.abs(x - np.median(x))))
+
+
+def split_by_length(
+    sd_short: float, sd_long: float, short_m: float, long_m: float
+) -> tuple[float, float]:
+    """Split a scale spread measured on windows of two lengths into a part that falls as 1/length
+    and a part that does not: sd(L)^2 = persistent^2 + (c / L)^2. Endpoint noise, from either
+    track, falls as 1/L. Returns (the 1/L part at short_m, the persistent part), each clipped at 0
+    when the two spreads are inconsistent with the model."""
+    k = (long_m / short_m) ** 2
+    falling_short = (sd_short**2 - sd_long**2) * k / (k - 1)
+    persistent = (k * sd_long**2 - sd_short**2) / (k - 1)
+    return float(np.sqrt(max(falling_short, 0.0))), float(np.sqrt(max(persistent, 0.0)))
+
+
+# A3: along and across travel
+
+
+def along_across(err: np.ndarray, ref_disp: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    """Horizontal error (N, 3) split into its signed component along the reference's horizontal
+    displacement (N, 3) and the signed component across it. With u the unit travel direction and e
+    the error, both as (x, z): along = u . e and across = u_z e_x - u_x e_z."""
+    u = ref_disp[:, [0, 2]]
+    n = np.linalg.norm(u, axis=1)
+    if np.any(n == 0):
+        raise ValueError("a window with no horizontal reference displacement has no direction")
+    u = u / n[:, None]
+    e = err[:, [0, 2]]
+    along = (e * u).sum(axis=1)
+    across = u[:, 1] * e[:, 0] - u[:, 0] * e[:, 1]
+    return along, across
+
+
+# A4: a range to the meter
+
+
+def range_correct(est: np.ndarray, rng_m: np.ndarray) -> np.ndarray:
+    """Move each estimated displacement from the anchor (N, 3) radially in the horizontal plane,
+    keeping its height, until its 3-D distance from the anchor equals the measured range. A range
+    shorter than the height difference puts the point directly above or below the anchor."""
+    h = horizontal_norm(est)
+    if np.any(h == 0):
+        raise ValueError(
+            "an estimate on the anchor's vertical has no horizontal direction to move along"
+        )
+    target = np.sqrt(np.maximum(np.asarray(rng_m) ** 2 - est[:, 1] ** 2, 0.0))
+    out = est.copy()
+    out[:, [0, 2]] *= (target / h)[:, None]
+    return out

experiments/drift-anatomy/pyproject.toml

@@ -0,0 +1,24 @@
+[project]
+name = "drift-anatomy"
+version = "0.1.0"
+description = "Which part of ARKit's position error is worth attacking: loop closure, local scale, along versus across travel, and a range to the meter"
+requires-python = ">=3.12"
+dependencies = ["numpy>=2.1", "scipy>=1.14"]
+
+[dependency-groups]
+dev = ["pytest>=8.3", "ruff>=0.8"]
+
+[tool.uv]
+package = false
+
+[tool.ruff]
+line-length = 100
+
+[tool.ruff.lint]
+select = ["E", "W", "F", "I", "B", "UP", "SIM", "RUF"]
+ignore = ["E501"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+pythonpath = ["."]
+addopts = ["--strict-markers", "--strict-config"]