Fix mathematical correctness across modules; add errata + verification anchors

CHANGELOG.md

@@ -1,5 +1,29 @@
 # Changelog
 
+## [0.1.1] - 2026-05-31
+
+### Fixed (mathematical correctness audit)
+Independent adversarial re-audit of every quantitative claim; full record in
+[ERRATA.md](ERRATA.md). Highlights:
+- **Critical:** corrected Fefferman–Graham `g₄` (spurious `1/(n−4)` on the
+  algebraic term; `n=4` trace coefficient `⅛ → 1/16`) and renormalized-volume
+  `v₂` (`−1/(n−2)J → −½J`, and the `n=2` sign). (`C1`, `C2`, `M11`, `M12`)
+- **Major:** the conformal Laplacian `P₂` now includes its curvature term
+  `−(n−2)/(4(n−1))R` (was the bare Laplacian); the `n=3` holographic stress
+  tensor is now actually traceless; the discrete Bach proxy uses the
+  mass-inverse `M⁻¹L M⁻¹L` bi-Laplacian; removed a false discrete `Q₄` identity;
+  relabeled discrete `H²−K` as the Willmore integrand (not 4D GJMS `Q₄`);
+  corrected the `P₆`/README `−Δ³ → +Δ³` sign. (`M2`, `M3`, `M13`–`M17`)
+- **Documented honestly (incomplete math flagged in-code, open for
+  contribution):** extrinsic `Q₄`/`P₂` (`M4`, `M8`), `L₂` Fialkow terms (`M7`),
+  weight-4 invariant under-count (`M5`, `M6`), Carroll connection scope and
+  symmetry predicate (`M9`, `M10`).
+- **Minor/style:** tractor `μ`-slot weight label `+1`; `bach()` `1/(n−3)`
+  docstring; `is_valid` rank check; removed dead/incorrect symmetrization line;
+  `Q₂` normalization caveat. (`m1`–`m6`, `s1`)
+- Added [ERRATA.md](ERRATA.md), [docs/TOOLING_GAPS.md](docs/TOOLING_GAPS.md), and
+  README "Verify it yourself" / "Contributing corrections" sections.
+
 ## [0.1.0] - 2026-05-11
 
 ### Added
@@ -12,4 +36,4 @@
 - Discrete conformal features (PyTorch): mesh utilities, Gaussian/mean curvature, discrete Q-curvature, discrete Bach norm, Willmore density, cross-ratios, conformal factor via Yamabe flow
 - Feature extraction pipeline: mesh_conformal_features with rotation-invariant features
 - Benchmark scripts: ShapeNet classification, SHREC retrieval, FAUST correspondence
-- 125 tests across both tracks
+- 160 tests across both tracks (157 original + 3 errata regression anchors)

README.md

@@ -9,25 +9,31 @@
 
 Two Python packages for computing conformal invariants — geometric quantities unchanged by local stretching — both symbolically (exact formulas via SageMath) and numerically on triangle meshes (GPU-ready via PyTorch). Implements tractor calculus, GJMS operators, Q-curvature, Blitz's conformal fundamental forms, Carroll geometry, and Fefferman-Graham holographic data.
 
+> **This library was independently re-audited for mathematical correctness in May 2026.** Every confirmed error — and exactly how we caught it — is documented in **[ERRATA.md](ERRATA.md)**. We publish the full record on purpose: a teaching tool earns trust by showing its work, mistakes included. See **["Verify it yourself"](#verify-it-yourself)** and **["Contributing corrections"](#contributing-corrections)**.
+
 ---
 
 ## Quick Start
 
-Compute Q-curvature exactly on the round 4-sphere, then verify the discrete approximation converges:
+Compute the Branson Q-curvature exactly on the round 4-sphere, then extract a discrete Willmore feature on a mesh:
 
 ```python
 # Symbolic: exact Q₄ on S⁴ via SageMath
 from conformal_toolkit import ConformalStructure
 cs = ConformalStructure(g_sphere4)
-cs.q_curvature(order=4)  # → 6  (exact, matches 2(n-1)! / ((n/2-1)!)² for n=4)
+cs.q_curvature(order=4)  # → 6   (exact; Branson's Q_n(Sⁿ) = (n-1)! = 3! = 6)
 
-# Discrete: Q₄ on an icosphere mesh via PyTorch
+# Discrete: a 4th-order surface feature on an icosphere mesh via PyTorch.
+# NOTE: this returns the Willmore integrand H² − K, NOT the 4D GJMS Q₄
+# (a 2-surface quantity cannot reproduce the 4-manifold value 6 — see ERRATA M15/M16).
 from conformal_features.discrete.q_curvature import discrete_q_curvature
 Q4 = discrete_q_curvature(vertices, faces, order=4)
-Q4.mean()  # → 5.94  (converges to 6 as mesh refines)
+Q4.mean()  # → 0   (H² − K vanishes on a round sphere of any radius; → 0 under refinement)
 ```
 
-One toolkit, two representations of the same invariant — symbolic formulas ground-truth the discrete features.
+The symbolic formulas ground-truth the geometry; the discrete features are
+mesh-domain analogues — *not* always the same number, and the docs now say which
+is which.
 
 ---
 
@@ -153,7 +159,7 @@ The 10 features per vertex:
 | 0 | Conformal factor | Conformal | From discrete Yamabe flow |
 | 1 | Willmore density | Conformal | H² (distance from minimality) |
 | 2 | Q₂ | Conformal | Discrete scalar curvature 2K |
-| 3 | Q₄ | Conformal | Higher-order curvature |
+| 3 | Willmore density H²−K | Conformal (integral) | 4th-order surface feature — the Willmore integrand, *not* the 4D GJMS Q₄ ([ERRATA M15](ERRATA.md)) |
 | 4 | Bach norm | Conformal | Bi-Laplacian proxy for non-conformal-flatness |
 | 5–6 | Cross-ratio stats | Möbius | Edge cross-ratio mean and variance |
 | 7 | Gaussian curvature | Isometric | Intrinsic curvature K |
@@ -240,7 +246,7 @@ conformal-toolkit/
 │   ├── discrete/               # Curvature, Q, Bach, Willmore, cross-ratios, Yamabe, spectral
 │   ├── features/               # mesh_conformal_features() pipeline
 │   └── benchmarks/             # ShapeNet, SHREC, FAUST evaluation (WIP)
-├── tests/                      # 157 tests across both packages
+├── tests/                      # 160 tests across both packages
 ├── examples/                   # 6 Jupyter notebooks
 └── paper.md                    # JOSS paper draft
 ```
@@ -273,7 +279,7 @@ conformal-toolkit/
 | `export` | `conformal_feature_vector(cs)` | Dict of all invariants at a point |
 | | `tensor_to_numpy(T)` | SageMath tensor → NumPy array |
 
-†P₆ computes the leading term (−Δ³) only; exact on conformally flat metrics. Obstruction at n=6 is a leading-order approximation.
+†P₆ computes the leading term (+Δ³) only; exact on conformally flat metrics. Obstruction at n=6 is a leading-order approximation (the Graham–Hirachi normalization constant is not applied — see [ERRATA m3](ERRATA.md)).
 
 ### conformal_features (PyTorch)
 
@@ -344,6 +350,66 @@ pytest tests/test_discrete/ tests/test_features/ -v
 
 ---
 
+## Verify it yourself
+
+Don't take our word for any formula — the whole point of a symbolic toolkit is
+that you can check it. Every correction in [ERRATA.md](ERRATA.md) was caught by
+evaluating a claim on a geometry where the answer is known in closed form. Here
+are the anchors we use; copy them into a Sage session and confirm:
+
+```python
+from sage.all import Manifold, sin
+from conformal_toolkit import ConformalStructure
+
+# --- Anchor 1: the round 4-sphere, where Branson's Q_n(Sⁿ) = (n-1)! ---
+# On Sⁿ (sectional curvature 1):  Ric = (n-1)g,  R = n(n-1),
+#   Schouten P = ½g,  J = tr P = n/2,  so  Q₄ = -ΔJ - 2|P|² + (n/2)J² = 6.
+cs = ConformalStructure(g_sphere4)
+assert cs.q_curvature(order=4) == 6          # = (4-1)! = 3!  (ERRATA M1)
+
+# --- Anchor 2: the conformal Laplacian carries a curvature term ---
+# P₂ f = Δf - (n-2)/(4(n-1)) R f.  The R-term is NOT optional for n > 2.
+# On S⁴ its coefficient is n(n-2)/4 = 2, never zero.  (ERRATA M2)
+
+# --- Anchor 3: Fefferman-Graham on the hyperbolic filling of Sⁿ ---
+# g_ρ = (1 - ρ²/4)² g₀  ⟹  g₂ = -½ g₀  and  g₄ = 1/16 g₀ exactly.  (ERRATA C1/M11)
+
+# --- Anchor 4 (discrete, PyTorch only): validate on a NON-constant field ---
+# The cotangent Laplacian L is a *stiffness* matrix; for f = x² on a flat mesh,
+# (L f)_i = -2·A_i  while  (M⁻¹ L f)_i = -2  recovers the pointwise Laplacian.
+# "It vanishes on a sphere" is a false positive — constants are annihilated by
+# any linear operator.  (ERRATA M17)
+```
+
+The method generalizes: **reduce a tensor claim to a scalar on a known geometry,
+and check conformal weights as a free checksum.** That single discipline caught
+most of the errata.
+
+---
+
+## Contributing corrections
+
+We would rather be corrected than be wrong, and this repository is built to make
+that easy. **Finding an error we missed is the system working — please send it.**
+
+1. Open an issue titled `Errata: <one-line claim>`.
+2. Give the counter-evidence the way we give ours: a *concrete geometry* (a
+   sphere radius, a flat patch, an explicit metric) on which the claim returns
+   the wrong number, or a *conformal-weight* argument that the terms can't match.
+   A failing check on a named anchor metric is the gold standard.
+3. Propose the corrected formula, stating your normalization convention (Branson
+   vs. analyst signs differ — half of conformal geometry's "errors" are
+   convention clashes), with a reference if you have one.
+4. If you can, add a regression test under `tests/` pinning the right value on
+   the anchor. *Verified-on-an-anchor beats argued-in-prose.*
+
+Open problems where we explicitly want help are listed at the end of
+[ERRATA.md](ERRATA.md) — the complete weight-4 hypersurface invariant basis
+(`M5`), the Fialkow/Weyl terms in `L₂` (`M7`), the full extrinsic `Q₄` (`M8`),
+and the FG `g₄` Bach differential terms (`C1`).
+
+---
+
 ## Citation
 
 ```bibtex

conformal_features/discrete/bach.py

@@ -1,42 +1,63 @@
 """Discrete approximation of Bach tensor norm.
 
-|B| ~ |L^2(K)| as a 4th-order curvature feature.
-This is not the true Bach tensor but a computationally tractable proxy
-that captures 4th-order curvature variation.
-
-The bi-Laplacian L^2 K (cotangent Laplacian applied twice to the discrete
-Gaussian curvature) acts as a 4th-order differential operator. Since
-discrete_gaussian_curvature already normalizes K = angle_defect / area, the
-raw bi-Laplacian (no additional area weighting) gives values that are small on
-a smooth, nearly-uniform sphere and converge to zero under mesh refinement.
+This is NOT the true Bach tensor (the 4D Fefferman-Graham obstruction) but a
+computationally tractable scalar PROXY that captures 4th-order curvature
+variation on a surface, used as one channel of the per-vertex feature vector.
+
+Discretization honesty note (see ERRATA M17). The cotangent matrix L is the
+FEM *stiffness* matrix, i.e. an INTEGRATED operator: (L f)_i ~ A_i (Delta f)_i,
+not (Delta f)_i. So the quantity computed here, |L L K|, is a *stiffness-
+weighted* (area-integrated) bi-Laplacian, NOT the pointwise bi-Laplacian. The
+true pointwise operator is
+
+    Delta^2 = (M^{-1} L)(M^{-1} L) = M^{-1} L M^{-1} L,   M = diag(vertex areas).
+
+We deliberately keep the integrated form |L L K| as the feature, because it is
+numerically well-behaved across mesh scales; but note two caveats:
+  * its *small absolute magnitude on a sphere is partly a scaling artifact*
+    (the A_i^2 weighting suppresses it), not proof of correctness -- the
+    earlier docstring's claim "K is already area-normalized, so no extra
+    weighting is needed" was the wrong reason for the right-ish behavior;
+  * the mass-lumped pointwise Delta^2 above is the mathematically correct
+    operator but amplifies coarse-mesh curvature noise (e.g. the 12 pentagonal
+    defects of an icosphere), so it is offered as `pointwise=True` rather than
+    the default. Always validate a discretization on a NON-constant field such
+    as f = x^2, where (M^{-1} L)(x^2) = -2 recovers the pointwise Laplacian
+    while (L)(x^2) = -2 A_i does not.
 """
 from __future__ import annotations
 import torch
 from conformal_features.discrete.curvature import discrete_gaussian_curvature
-from conformal_features.discrete.mesh_utils import cotangent_laplacian
+from conformal_features.discrete.mesh_utils import cotangent_laplacian, vertex_areas
 
 
-def discrete_bach_norm(vertices: torch.Tensor, faces: torch.Tensor) -> torch.Tensor:
-    """Per-vertex discrete Bach tensor norm approximation.
-
-    Computes |L @ L @ K| where K is the discrete Gaussian curvature and L is
-    the cotangent Laplacian. This bi-Laplacian captures 4th-order curvature
-    variation and is small on a sphere (where K is nearly constant).
+def discrete_bach_norm(
+    vertices: torch.Tensor, faces: torch.Tensor, pointwise: bool = False
+) -> torch.Tensor:
+    """Per-vertex discrete Bach tensor norm proxy.
 
     Args:
         vertices: (V, 3) Tensor
         faces: (F, 3) LongTensor
+        pointwise: if False (default), return the stiffness-weighted
+            (area-integrated) bi-Laplacian |L L K| -- scale-stable, used as the
+            feature channel. If True, return the mathematically-correct
+            pointwise bi-Laplacian |M^{-1} L M^{-1} L K| (noisier on coarse
+            meshes). See ERRATA M17.
 
     Returns:
         bach_norm: (V,) Tensor
     """
     K = discrete_gaussian_curvature(vertices, faces)
     L = cotangent_laplacian(vertices, faces)
 
-    # Apply the cotangent Laplacian twice. K is already area-normalized
-    # (angle_defect / area), so no additional area weighting is needed —
-    # the raw bi-Laplacian L@L@K converges to zero on a sphere where K is constant.
-    LK = torch.sparse.mm(L, K.unsqueeze(1)).squeeze(1)
-    LLK = torch.sparse.mm(L, LK.unsqueeze(1)).squeeze(1)
+    def stiffness(field: torch.Tensor) -> torch.Tensor:
+        return torch.sparse.mm(L, field.unsqueeze(1)).squeeze(1)
+
+    if pointwise:
+        inv_area = 1.0 / vertex_areas(vertices, faces).clamp_min(1e-12)
+        lap = lambda f: inv_area * stiffness(f)   # pointwise Delta = M^{-1} L
+        return lap(lap(K)).abs()
 
-    return LLK.abs()
+    # Integrated proxy |L L K| (the default feature channel).
+    return stiffness(stiffness(K)).abs()

conformal_features/discrete/q_curvature.py

@@ -9,13 +9,26 @@
 
 
 def discrete_q_curvature(vertices: torch.Tensor, faces: torch.Tensor, order: int = 2) -> torch.Tensor:
-    """Discrete Q-curvature of the given order.
-
-    Q_2 = R = 2K (scalar curvature = 2 * Gaussian curvature for surfaces)
-    Q_4 approximation: H^2 - K  (conformally invariant Willmore integrand).
-        On S^2: H = K = 1, so Q_4 = 0.
-        More precisely this approximates -Delta(J) - 2|P|^2 + J^2 in 2D via
-        the identity Q_4 = 2K^2 - 2KH^2 (rescaled), which equals 0 on S^2.
+    """Discrete curvature surface features (orders 2 and 4).
+
+    order=2:  Q_2 = R = 2K  (scalar curvature = 2 * Gaussian curvature
+              for a surface).
+
+    order=4:  returns the *Willmore integrand* H^2 - K, NOT the 4D GJMS /
+              Branson Q_4. Two honesty caveats (see ERRATA M14/M15/M16):
+        * Naming. The intrinsic 4th-order GJMS Q_4 is a 4-manifold object
+          with Q_4 = (n-1)! = 6 on the round S^4. A 2-surface quantity such
+          as H^2 - K cannot reproduce that value; it is a Willmore-type
+          surface feature, not the GJMS Q_4.
+        * Invariance. H^2 - K is conformally invariant only *under the
+          integral*: ∫(H^2 - K) dA is controlled by Gauss-Bonnet
+          (∫K dA = 2*pi*chi) plus Moebius-invariance of ∫H^2 dA. It is NOT
+          pointwise conformally invariant.
+        * Convergence. On a round S^2 of ANY radius R, H^2 = 1/R^2 = K, so
+          H^2 - K = 0; thus this feature -> 0 under mesh refinement on a
+          sphere (it does NOT converge to 6 -- an earlier README claim).
+          (The earlier "identity" Q_4 = 2K^2 - 2KH^2 was false: that equals
+          -2K (H^2 - K), not a rescaling of H^2 - K.)
 
     Args:
         vertices: (V, 3) Tensor
@@ -32,7 +45,7 @@ def discrete_q_curvature(vertices: torch.Tensor, faces: torch.Tensor, order: int
 
     if order == 4:
         H = discrete_mean_curvature(vertices, faces)
-        # Q_4 ~ H^2 - K: conformally invariant, vanishes on S^2 (H=K=1)
+        # Willmore integrand H^2 - K; -> 0 on a round sphere of any radius.
         return H ** 2 - K
 
     raise ValueError(f"Order {order} not supported (only 2 and 4)")

conformal_toolkit/carroll/bms.py

@@ -58,6 +58,17 @@ def is_bms_symmetry(carroll_struct, vector_field):
     A vector field ξ is a **strict Carroll symmetry** if £_ξ h = 0, and a
     **conformal Carroll symmetry** if £_ξ h = λ h for some scalar λ.
 
+    INCOMPLETE PREDICATE -- necessary but NOT sufficient (see ERRATA M10).
+    A Carroll structure is the PAIR (v, h), so a symmetry must preserve BOTH:
+        strict:    £_ξ h = 0      and  £_ξ v = 0
+        conformal: £_ξ h = 2λ h   and  £_ξ v = -λ v   (Duval-Gibbons-Horvathy)
+    Because h is degenerate (h(v, ·) = 0), the £_ξ h condition places NO
+    constraint on the v-direction, so this routine -- which only tests
+    £_ξ h -- can wrongly accept a ξ that moves v out of ker(h). A correct
+    predicate must ALSO verify £_ξ v ∝ v (strict: = 0; conformal: = -λ v).
+    (The supertranslation generator ξ = f·v built above does satisfy
+    £_ξ v = 0, so it passes; the gap bites for general user input.)
+
     The Lie derivative is computed component-wise:
 
         (£_ξ h)_{ij} = ξ^k ∂_k h_{ij} + h_{kj} ∂_i ξ^k + h_{ik} ∂_j ξ^k

conformal_toolkit/carroll/carroll_connection.py

@@ -1,15 +1,23 @@
 """Carroll connection: Christoffel-like symbols adapted to a Carroll structure.
 
-For a Carroll manifold with adapted coordinates (t, x^i) where v = ∂_t, the
-Carroll connection has the block structure:
+For a Carroll manifold with adapted coordinates (t, x^i) where v = ∂_t, this
+module uses the block structure:
 
     Γ^t_{μν} = 0   (all)
     Γ^i_{tt} = 0
     Γ^i_{tj} = 0
     Γ^i_{jk} = spatial Christoffel symbols from h_{jk}
 
-The spatial Christoffel symbols are the Levi-Civita symbols of the restriction
-of h to the spatial leaves (where h is non-degenerate).
+SCOPE CAVEAT (see ERRATA M9). Carrollian connections are NOT unique (h is
+degenerate), and the choice Γ^i_{tj} = 0 is correct ONLY when the Carroll
+"electric field" E_{ij} = (1/2)(£_v h)_{ij} vanishes -- i.e. when the spatial
+metric is time-independent (∂_t h_{ij} = 0). When ∂_t h_{ij} ≠ 0, metric
+compatibility ∇h = 0 forces Γ^k_{t(i} h_{j)k} = (1/2) ∂_t h_{ij} ≠ 0, so the
+correct symmetric part is Γ^i_{(tj)} = -(1/2) h^{ik}(£_v h)_{kj}. This module
+otherwise computes a generically NONZERO carroll_electric_field, so presenting
+Γ^i_{tj} = 0 as THE Carroll connection is internally inconsistent in the
+time-dependent case. Treat the symbols below as the preserved-h
+(£_v h = 0) representative.
 """
 
 from sage.all import SR

conformal_toolkit/carroll/carroll_structure.py

@@ -58,12 +58,21 @@ def is_valid(self):
         """Check the Carroll compatibility condition h(v, ·) = 0.
 
         Contracts h with v on its first slot and verifies that all components
-        of the resulting 1-form vanish symbolically.
+        of the resulting 1-form vanish symbolically. Also checks that the
+        spatial metric h has rank n (a 1-dimensional kernel spanned by v),
+        as required by the definition of a Carroll structure.
+
+        INCOMPLETENESS NOTE (ERRATA m5): the original check verified only
+        h(v, ·) = 0 and NOT rank(h) = n, so a fully-zero h or any
+        rank-deficient h (kernel dimension > 1) passed validation despite
+        the "rank-n degenerate" definition. The rank check below closes that
+        gap (downstream spatial_christoffel also implicitly needs the spatial
+        block invertible).
 
         Returns
         -------
         bool
-            True if v lies in the kernel of h, False otherwise.
+            True if v lies in the kernel of h and rank(h) = n, else False.
         """
         h = self._h
         v = self._v
@@ -85,6 +94,17 @@ def is_valid(self):
             except Exception:
                 if not bool(val == 0):
                     return False
+
+        # Rank check: h must be degenerate of rank exactly n = dim - 1
+        # (kernel = span(v), one-dimensional). See ERRATA m5.
+        from sage.all import matrix, SR
+        h_mat = matrix(SR, self._dim, self._dim)
+        for a in range(self._dim):
+            for b in range(self._dim):
+                comp = h[frame, a, b]
+                h_mat[a, b] = comp.expr() if hasattr(comp, 'expr') else SR(comp)
+        if h_mat.rank() != self._dim - 1:
+            return False
         return True
 
     # ------------------------------------------------------------------