fix(#304): implement ReferenceConstraint residuals and verify post-solve

CHANGES.md

@@ -7,6 +7,11 @@ issue tracker.  The initial project development roadmap is documented here:
 
 ## Unreleased
 
+### Fixed
+
+- `ReferenceConstraint.evaluate()` computes real residuals; verified post-solve. (#304)
+
+
 ## Release v0.11.2
 
 Released 2026-06-22

src/ad_hoc_diffractometer/forward.py

@@ -4234,6 +4234,7 @@ def _validate_solutions(
     from .display import precision_atol
     from .mode import BisectConstraint
     from .mode import ConstraintViolation
+    from .mode import ReferenceConstraint
     from .mode import SampleConstraint
 
     if not solutions:
@@ -4243,10 +4244,15 @@ def _validate_solutions(
 
     for i, sol in enumerate(solutions):
         for c in mode.constraints:
-            # Only validate constraints that have a numeric residual we can check.
-            # BisectConstraint and SampleConstraint are the checkable ones;
-            # DetectorConstraint is applied by the solver directly (ttheta_deg);
-            # ReferenceConstraint is not yet implemented.
+            # Validate constraints that expose a numeric per-solution residual.
+            # BisectConstraint and SampleConstraint are checked against the
+            # motor angles directly; DetectorConstraint is applied by the
+            # solver (ttheta_deg) so it never appears here.
+            # ReferenceConstraint pseudo-angle residuals
+            # (incidence/emergence/specular/omega/naz) are checked too — except
+            # "psi", which is enforced upstream as a validation filter in
+            # _solve_psi_mode (its motor-frame residual is subject to ±360°
+            # azimuthal wraparound and is not a reliable per-solution check).
             if isinstance(c, BisectConstraint | SampleConstraint):
                 try:
                     residual = c.evaluate(sol, geometry)
@@ -4272,6 +4278,29 @@ def _validate_solutions(
                         residual=residual,
                         tolerance=atol,
                     )
+            elif isinstance(c, ReferenceConstraint):
+                # "psi" is enforced by the validation filter, not here.
+                if c.name == "psi":
+                    continue
+                # Skip when the required reference vector is unset (evaluate()
+                # would raise ValueError); is_implemented() already gated the
+                # solve, but a mode could in principle reach here without one.
+                if not c.has_reference_vector(geometry):
+                    continue
+                try:
+                    residual = c.evaluate(sol, geometry)
+                except (KeyError, ValueError):
+                    continue
+                if abs(residual) > atol:
+                    raise ConstraintViolation(
+                        solution_index=i,
+                        constraint_repr=(
+                            f"ReferenceConstraint({c.name!r}, {c.value!r}): "
+                            f"residual {residual!s}° exceeds tolerance"
+                        ),
+                        residual=residual,
+                        tolerance=atol,
+                    )
 
 
 def _check_limits(

src/ad_hoc_diffractometer/mode.py

@@ -918,28 +918,77 @@ def evaluate(
         geometry: AdHocDiffractometer,
     ) -> float:
         """
-        Return the constraint residual (zero when satisfied).
+        Return the constraint residual in degrees (zero when satisfied).
+
+        The residual is the difference between the physical pseudo-angle
+        computed from the supplied motor angles (via the standalone
+        functions in :mod:`~ad_hoc_diffractometer.reference`) and the
+        target :attr:`value`:
+
+        - ``"incidence"`` — ``incidence_angle - value``
+        - ``"emergence"`` — ``emergence_angle - value``
+        - ``"specular"`` — ``incidence_angle - emergence_angle``
+          (the relational condition α_i = α_f; :attr:`value` is the
+          boolean flag ``True`` and does not enter the residual)
+        - ``"psi"`` — ``psi_angle - value``
+        - ``"naz"`` — ``naz_angle - value``
+        - ``"omega"`` — ``omega_pseudo - value``
+
+        Parameters
+        ----------
+        angles : dict[str, float]
+            Current motor angles in degrees, keyed by stage name.
+        geometry : AdHocDiffractometer
+            The diffractometer.  The required reference vector
+            (``surface_normal`` for incidence/emergence/specular,
+            ``azimuth`` for psi/naz) must be set, except for ``"omega"``
+            which needs no reference vector.
+
+        Returns
+        -------
+        float
+            Residual in degrees.  Zero means the constraint is satisfied.
 
-        Not yet implemented — all reference constraints require the
-        reference vector infrastructure (Issue J).  Raises
-        ``NotImplementedError`` when called.
+        Raises
+        ------
+        ValueError
+            If the required reference vector is not set on the geometry
+            (raised by the underlying
+            :mod:`~ad_hoc_diffractometer.reference` function).
         """
-        raise NotImplementedError(
-            f"ReferenceConstraint('{self._name}') solver is not yet implemented. "
-            "Reference constraint solvers require the reference vector "
-            "infrastructure (Issue J)."
-        )
+        # Local import to avoid module-load ordering issues (reference.py
+        # imports geometry types lazily; mirror the pattern used elsewhere
+        # in this package).
+        from . import reference
+
+        if self._name == "incidence":
+            return reference.incidence_angle(geometry, angles=angles) - float(
+                self._value
+            )
+        if self._name == "emergence":
+            return reference.emergence_angle(geometry, angles=angles) - float(
+                self._value
+            )
+        if self._name == "specular":
+            # Relational: incidence == emergence.  value is the flag True.
+            return reference.incidence_angle(
+                geometry, angles=angles
+            ) - reference.emergence_angle(geometry, angles=angles)
+        if self._name == "psi":
+            return reference.psi_angle(geometry, angles=angles) - float(self._value)
+        if self._name == "naz":
+            return reference.naz_angle(geometry, angles=angles) - float(self._value)
+        # omega
+        return reference.omega_pseudo(geometry, angles=angles) - float(self._value)
 
     def is_satisfied(
         self,
         angles: dict[str, float],
         geometry: AdHocDiffractometer,
         tol: float = 1e-6,
     ) -> bool:
-        """Return True when the constraint is satisfied (not yet implemented)."""
-        raise NotImplementedError(
-            f"ReferenceConstraint('{self._name}') is not yet implemented."
-        )
+        """Return True when ``|evaluate(angles, geometry)| < tol``."""
+        return abs(self.evaluate(angles, geometry)) < tol
 
     def has_reference_vector(self, geometry: AdHocDiffractometer) -> bool:
         """

tests/test_mode.py

@@ -797,39 +797,145 @@ def test_reference_constraint_to_dict_from_dict_a_eq_b():
 
 
 @pytest.mark.parametrize(
-    "name, value, surface_normal, expected",
+    "name, value, surface_normal, azimuth, expected",
     [
         # incidence/emergence/specular: implemented when surface_normal is set
-        pytest.param("incidence", 0.0, (0, 0, 1), True, id="incidence-with-sn"),
-        pytest.param("incidence", 0.0, None, False, id="incidence-no-sn"),
-        pytest.param("emergence", 0.0, (0, 0, 1), True, id="emergence-with-sn"),
-        pytest.param("specular", True, (0, 0, 1), True, id="specular-with-sn"),
-        pytest.param("specular", True, None, False, id="specular-no-sn"),
-        # psi/naz: never implemented (no forward solver yet)
-        pytest.param("psi", 0.0, (0, 0, 1), False, id="psi-not-implemented"),
-        pytest.param("naz", 0.0, (0, 0, 1), False, id="naz-not-implemented"),
+        pytest.param("incidence", 0.0, (0, 0, 1), None, True, id="incidence-with-sn"),
+        pytest.param("incidence", 0.0, None, None, False, id="incidence-no-sn"),
+        pytest.param("emergence", 0.0, (0, 0, 1), None, True, id="emergence-with-sn"),
+        pytest.param("specular", True, (0, 0, 1), None, True, id="specular-with-sn"),
+        pytest.param("specular", True, None, None, False, id="specular-no-sn"),
+        # psi: implemented (validation filter) when azimuth is set
+        pytest.param("psi", 0.0, None, (0, 0, 1), True, id="psi-with-azimuth"),
+        pytest.param("psi", 0.0, None, None, False, id="psi-no-azimuth"),
+        # naz: no forward solver yet
+        pytest.param("naz", 0.0, (0, 0, 1), None, False, id="naz-not-implemented"),
     ],
 )
-def test_reference_constraint_is_implemented(name, value, surface_normal, expected):
-    """ReferenceConstraint.is_implemented() reflects surface_normal presence and solver."""
+def test_reference_constraint_is_implemented(
+    name, value, surface_normal, azimuth, expected
+):
+    """ReferenceConstraint.is_implemented() reflects reference vector and solver."""
     g = _psic_like()
     g._surface_normal = surface_normal  # noqa: SLF001
+    g._azimuth = azimuth  # noqa: SLF001
     rc = ReferenceConstraint(name, value)
     assert rc.is_implemented(g) is expected
 
 
-def test_reference_constraint_evaluate_raises():
-    g = _psic_like()
-    rc = ReferenceConstraint("psi", 90.0)
-    with pytest.raises(NotImplementedError, match="not yet implemented"):
-        rc.evaluate({}, g)
+def _psic_reference_geometry(**kwargs):
+    """Real psic geometry with wavelength, cubic UB, and a reference vector."""
+    import ad_hoc_diffractometer as ahd
+    from ad_hoc_diffractometer import ub_identity
 
+    g = ahd.make_geometry("psic")
+    g.wavelength = 1.54
+    g.sample.lattice = ahd.Lattice(a=4.0)
+    ub_identity(g.sample)
+    for key, val in kwargs.items():
+        setattr(g, key, val)
+    return g
 
-def test_reference_constraint_is_satisfied_raises():
-    g = _psic_like()
-    rc = ReferenceConstraint("psi", 90.0)
-    with pytest.raises(NotImplementedError):
-        rc.is_satisfied({}, g)
+
+@pytest.mark.parametrize(
+    "name, value, kwargs, context",
+    [
+        pytest.param(
+            "incidence",
+            0.0,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="incidence-residual",
+        ),
+        pytest.param(
+            "emergence",
+            0.0,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="emergence-residual",
+        ),
+        pytest.param(
+            "specular",
+            True,
+            {"surface_normal": (0, 0, 1)},
+            does_not_raise(),
+            id="specular-residual",
+        ),
+        pytest.param(
+            "omega",
+            0.0,
+            {},
+            does_not_raise(),
+            id="omega-residual",
+        ),
+        pytest.param(
+            "incidence",
+            0.0,
+            {},
+            pytest.raises(ValueError, match=re.escape("surface_normal must be set")),
+            id="incidence-no-surface-normal-raises",
+        ),
+    ],
+)
+def test_reference_constraint_evaluate(name, value, kwargs, context):
+    """ReferenceConstraint.evaluate() returns a real pseudo-angle residual.
+
+    A solution returned by forward() satisfies its reference constraint, so
+    the residual evaluated at that solution must be ~0.  evaluate() raises
+    ValueError when the required reference vector is not set.
+    """
+    g = _psic_reference_geometry(**kwargs)
+    rc = ReferenceConstraint(name, value)
+    with context:
+        # Use the geometry's current (zeroed) angles for the raising case;
+        # for the satisfied cases use a real forward solution.
+        angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+        if "surface_normal" in kwargs or name == "omega":
+            g.mode_name = {
+                "incidence": "fixed_incidence_vertical",
+                "emergence": "fixed_emergence_vertical",
+                "specular": "specular_vertical",
+                "omega": "fixed_omega_vertical",
+            }[name]
+            sols = g.forward(1, 0, 1)
+            angles = sols[0]
+        residual = rc.evaluate(angles, g)
+        assert residual == pytest.approx(0.0, abs=1e-6)
+
+
+def test_reference_constraint_is_satisfied():
+    """is_satisfied() is True at a satisfying solution, False for a wrong target."""
+    g = _psic_reference_geometry(surface_normal=(0, 0, 1))
+    g.mode_name = "fixed_incidence_vertical"
+    sol = g.forward(1, 0, 1)[0]
+    assert ReferenceConstraint("incidence", 0.0).is_satisfied(sol, g) is True
+    assert ReferenceConstraint("incidence", 99.0).is_satisfied(sol, g) is False
+
+
+def test_reference_constraint_evaluate_psi_branch():
+    """evaluate() for the 'psi' branch returns psi_angle - target."""
+    import ad_hoc_diffractometer.reference as reference
+
+    g = _psic_reference_geometry(azimuth=(0, 0, 1))
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+    angles["delta"] = 30.0  # ensure a non-degenerate Q for psi
+    angles["eta"] = 15.0
+    expected = reference.psi_angle(g, angles=angles) - 12.0
+    assert ReferenceConstraint("psi", 12.0).evaluate(angles, g) == pytest.approx(
+        expected
+    )
+
+
+def test_reference_constraint_evaluate_naz_branch():
+    """evaluate() for the 'naz' branch returns naz_angle - target."""
+    import ad_hoc_diffractometer.reference as reference
+
+    g = _psic_reference_geometry(surface_normal=(1, 0, 0))
+    angles = {s.name: s.angle for s in g._stages.values()}  # noqa: SLF001
+    expected = reference.naz_angle(g, angles=angles) - 7.0
+    assert ReferenceConstraint("naz", 7.0).evaluate(angles, g) == pytest.approx(
+        expected
+    )
 
 
 # ---------------------------------------------------------------------------