diff --git a/docs/upstream-bugs.md b/docs/upstream-bugs.md index 0a3842f..70065f2 100644 --- a/docs/upstream-bugs.md +++ b/docs/upstream-bugs.md @@ -1,6 +1,7 @@ # Bugs found in the upstream Python `DiffusionGeometry` -Defects discovered in the reference implementation while porting it to Julia. +Defects discovered in the reference implementation while porting it to Julia — and, +in §5, in its *test suite*, where the shipped code is right but nothing checks it. - **Upstream:** https://github.com/Iolo-Jones/DiffusionGeometry - **Revision audited:** `f45b39f` ("update mdg notebook"), clean tree. @@ -10,7 +11,9 @@ Where a bug has a numerical parity target, the Julia port implements the *correc behaviour and `pyparity/gen_fixtures.py` stores a **corrected** reference, so the parity gate compares against the right answer rather than pinning the bug. Each correction is guarded: if upstream fixes the bug, fixture generation fails loudly -rather than silently double-correcting. See `PORTING_PLAN.md` §8. +rather than silently double-correcting. See `PORTING_PLAN.md` §8. A test defect has +no parity target and so no such guard — §5 is a coverage hole, closed on the Julia +side only. | # | Location | Severity | Status in the port | |---|---|---|---| @@ -18,6 +21,7 @@ rather than silently double-correcting. See `PORTING_PLAN.md` §8. | 2 | `tensors/vector_fields/vector_field.py::VectorField.from_reconstruction` | **Dead code** (`AttributeError`) | Reimplemented; round-trip gated | | 3 | `methods/geodesics.py::geodesic_distances_function` | **Dead code** (`AttributeError`) ×2 | Fixed; fixture corrected | | 4 | `utils/basis_utils.py::form_to_ambient_polyvector` (`k == 0` branch) | Cosmetic (unreachable) | N/A | +| 5 | `tests/test_src/test_hessian.py::test_hessian_02_sym_weak_matrix` | Test defect (coverage hole) | Corrected expectation; check enabled | --- @@ -207,3 +211,43 @@ a `FormSpace(dg, 0)` directly. **In the port.** Not reproduced. `to_ambient(::ScalarFunction)` is defined separately as `to_pointwise_basis(f)`, matching the reachable Python behaviour, and Julia's `FormSpace` constructor asserts `1 ≤ degree ≤ dim`. + + +--- + +## 5. The symmetric Hessian's value check is disabled because the *expectation* is wrong + +**Where:** `tests/test_src/test_hessian.py`, `test_hessian_02_sym_weak_matrix`. + +**What's wrong.** The test assembles the weak symmetric Hessian by hand and then does +not compare against it: + +```python +# Manual computation: H_sym_weak[i, s, I] = ∑_p μ[p] * φ_i[p] * H[p, j1(s), j2(s), I] +... +# NOTE: Disabling strict value check. +assert hessian_sym_02_weak_computed.shape == (n1 * d_sym, n0) +# assert np.allclose(hess_comp_sub, hess_sel_flat) +``` + +The expectation is missing the off-diagonal multiplicity factor. `hessian_02_sym_weak` +doubles the `j₁ ≠ j₂` entries, and it is right to: the symmetric basis element for +`j₁ ≠ j₂` is `dx_{j₁} ⊗ dx_{j₂} + dx_{j₂} ⊗ dx_{j₁}`, so pairing the (symmetric) +Hessian against it picks up both entries. The same convention is baked into +`gamma_02_sym`, whose off-diagonal rows and columns carry the matching ×2. Comparing +the builder against an expectation that drops the factor fails on exactly the +off-diagonal rows — so the assertion was commented out rather than the expectation +fixed, and `hessian_02_sym_weak` was left with no value check on either side of the +port. + +**Severity: test defect.** The shipped builder is correct; the coverage was not. But +it was the *only* weak-operator builder whose values nothing verified — the fixture +parity gate compares Julia against Python's output, which proves nothing if Python is +wrong. + +**In the port.** The check is enabled, with the factor restored, in +`test/test_pysrc.jl` → `hessian` → `hessian_02_sym_weak matrix`, and it passes at +`d = 1..4`. Because that test asserts the very convention in question, a second test +pins it convention-free: expanding the strong symmetric coefficients into the full +`(0,2)` basis and pairing them with the full Gram reproduces `hessian_02_weak`, whose +`⟨H(φ_I), φ_i dx_j ⊗ dx_k⟩` entries involve no symmetric-basis choice at all. diff --git a/test/test_pysrc.jl b/test/test_pysrc.jl index 7c9c0bc..ca26cbb 100644 --- a/test/test_pysrc.jl +++ b/test/test_pysrc.jl @@ -8,6 +8,10 @@ # # `d` is the package's exterior-derivative accessor, so the ambient dimension is # spelled `dim` throughout to avoid shadowing it. +# +# One test here has no upstream counterpart: `hessian` → `hessian_02_sym_weak agrees +# with the full (0,2) Hessian`, which corroborates the value check upstream disabled. +# It is marked as port-added where it appears. using Combinatorics: combinations # `coeffs` and `np_reshape` are internal; the tests need them to poke at the @@ -164,13 +168,52 @@ end # ── test_hessian.py ─────────────────────────────────────────────────────────── @testset "hessian" begin for_each_geom() do dg, (dim, n, n0, n1) - d_sym = size(get_symmetric_basis_indices(dim), 1) + sym_idx = get_symmetric_basis_indices(dim) + d_sym = size(sym_idx, 1) + u, mu = function_basis(dg), measure(dg) + H = hessian_functions(dg.cache) # Upstream asserts only the shape here — its value check is commented out - # (`NOTE: Disabling strict value check`). Ported as-is; see the note in the - # summary rather than silently strengthening it. - @testset "hessian_02_sym_weak shape" begin - @test size(weak(hessian(dg))) == (n1 * d_sym, n0) + # (`NOTE: Disabling strict value check`). The check is enabled below; the + # expectation upstream compares against is the thing that is wrong, not the + # builder. See docs/upstream-bugs.md §5. + @testset "hessian_02_sym_weak matrix" begin + H_weak = weak(hessian(dg)) + @test size(H_weak) == (n1 * d_sym, n0) + + # H_sym_weak[i, s, I] = w_s · Σ_p μ_p φ_i(p) H[p, j₁(s), j₂(s), I], with + # w_s = 2 off the diagonal: the symmetric basis element for j₁ ≠ j₂ is + # dx_{j₁} ⊗ dx_{j₂} + dx_{j₂} ⊗ dx_{j₁}, so it picks up both entries of + # the (symmetric) Hessian. Upstream's expectation drops that factor. + i_sel, s_sel, I_sel = sample_indices(n1), sample_indices(d_sym), sample_indices(n0) + manual = zeros(length(i_sel) * length(s_sel), length(I_sel)) + for (a, i) in enumerate(i_sel), (b, s) in enumerate(s_sel) + j1, j2 = sym_idx[s, 1], sym_idx[s, 2] + w = j1 == j2 ? 1.0 : 2.0 + for (c, I) in enumerate(I_sel) + manual[(a - 1) * length(s_sel) + b, c] = + w * sum(u[p, i] * H[p, j1, j2, I] * mu[p] for p in 1:n) + end + end + rows = vec([(i - 1) * d_sym + s for s in s_sel, i in i_sel]) + @test isapprox(H_weak[rows, I_sel], manual; atol=1e-6) + end + + # No upstream counterpart — added by the port. The ×2 above is a claim about + # the symmetric basis, so pin it without reference to that basis: expand the + # strong symmetric coefficients into the full (0, 2) basis, pair them with the + # full Gram, and the result must be the full weak Hessian + # ⟨H(φ_I), φ_i dx_j ⊗ dx_k⟩ — which carries no such convention. + @testset "hessian_02_sym_weak agrees with the full (0,2) Hessian" begin + c_sym = gram_inv(tensor02sym_space(dg)) * weak(hessian(dg)) + c_full = reduce(hcat, [expand_symmetric_tensor_coeffs(c_sym[:, I], n1, dim) + for I in 1:n0]) + + lhs = gram(tensor02_space(dg)) * c_full + rhs = hessian_02_weak(u, H, mu, n1) + rows = flat_idx(sample_indices(n1), dim * dim) + cols = sample_indices(n0) + @test isapprox(lhs[rows, cols], rhs[rows, cols]; atol=1e-6) end # Strong/weak relation: G · (G⁻¹ · H_weak) == H_weak.