feat(folds): validate P2Fold on Bermagui and document gauss4

.gitignore

@@ -137,4 +137,6 @@ conda/channeldata.json
 conda/noarch
 conda/win-64
 .DS_Store
-sandbox/*
+sandbox/*
+# Local validation outputs
+/validation/

LoopStructural/modelling/core/_model_feature_factory.py

@@ -223,8 +223,9 @@ def build_folded_foliation(
             ``"DFI"`` (default, piecewise linear) or ``"P2Fold"`` (piecewise quadratic).
             Other values are replaced by ``"DFI"`` with a warning.
         nelements : int
-            number of elements in the support. For ``"P2Fold"`` this is the number
-            of tetrahedra, e.g. 1e4 gives about 56k tets and 86k unknowns
+            requested structured-grid cell budget. Each cell splits into five
+            tetrahedra. For ``"P2Fold"``, a request of 1e4 gave about 56k tets
+            and 86k unknowns in the Bermagui validation; inspect the actual counts
         kwargs
             additional kwargs to be passed through to other functions
 
@@ -315,7 +316,8 @@ def build_folded_fold_frame(
         interpolatortype : str
             ``"DFI"`` or ``"P2Fold"`` sets the interpolator of the first (fold) coordinate.
             Any other value gives ``"DFI"``. The other two coordinates always use FDI.
-            For ``"P2Fold"``, ``nelements`` is the number of tetrahedra.
+            ``nelements`` is a structured-grid cell budget, not a literal
+            tetrahedron count. P2 adds edge nodes; inspect the actual dof.
         fold_frame : StructuralFrame, optional
             the fold frame for the fold if not specified uses last feature added
         nelements : int

agent/HANDOFF_p2_fold.md

@@ -0,0 +1,145 @@
+# P2 fold handoff: phase 4 validation
+
+## Current status - 20 October 2026
+
+Phases 0-3 are implemented on `master`. The original prototype narrative below
+is retained as historical context and is **superseded** for current results and
+weighting rules by the corrected measurements in this section.
+The library class and factory registration are in use; an engine swap is no
+longer required.
+
+The phase 4 decision gate passed on both targets with the same fixed upstream
+frames, seed 12345, explicit target buffer 0, pinned CG settings, and S0
+`limb_wl=None`:
+
+| Target | Method | Unknowns | Median axial error (deg) | Above 30 deg (%) |
+|---|---|---:|---:|---:|
+| S2 | DFI | 12,992 | 6.15 | 0.009 |
+| S2 | P2Fold gauss4 | 9,711 | 4.69 | 0.000 |
+| S0 | DFI | 112,560 | 20.63 | 37.282 |
+| S0 | P2Fold gauss4 | 85,673 | 11.14 | 19.379 |
+
+The comparison also covers the requested full P2 meshes (S2 1e4, S0 1e5), both
+sampling rules, tets, solve time, peak resident memory, input hashes and figures.
+The full S0 gauss4 run has 765,347 unknowns, median 8.69 deg, 16.1% above 30 deg,
+92.5 s solve time and about 6.6 GiB peak RSS. The remaining S0 error tail matters.
+Unregularised centroid is unstable (about 60 deg median), even with converged CG.
+
+Corrected rules: use `w*sqrt(q_i)` after row normalisation, not `w/n_samples`.
+The common volume factor cancels and is removed. The gauss4 squared residual
+weight is four times the prototype's at the same named family weight. Actual
+residuals can still differ between sampling rules. `nelements` is currently a
+structured-grid cell budget; realised P2 tets and unknowns exceed that request.
+The older statement that it is a literal tet count is inaccurate.
+
+These results are target-interpolator comparisons on a fixed **DFI** upstream
+chain. The old all-P2 S2/S0 table changed the S0 driver and profile, so its
+6.20 deg / 0.94% S0 result is not a direct before/after weight comparison.
+Field scale remains arbitrary. Buffer, F2, F3 and profile-fit issues remain open.
+Documentation and a bundled-data example are prepared; merge and post-merge
+release checks precede Phase 5. Root dependency ranges currently exclude the
+pending 0.1.0 workspace-package releases and need alignment before publication.
+
+---
+
+## Historical prototype handoff (superseded)
+
+
+## Goal
+Make the quadratic (P2) tetra interpolator usable for fold modelling (`create_and_add_folded_fold_frame`).
+Today P2 has no fold constraints, and the model factory forces `"DFI"` (linear) for folded features.
+
+## Why do it
+Tests on the Bermagui S2 and S0 folded foliations (`debug_s0/`). The metric is the angle between the unit
+gradient of the solved field and the direction expected from the fold constraints, `cross(fold_direction, fold_axis)`.
+
+| Mesh | Unknowns | Median error |
+|---|---|---|
+| Linear DFI | 12k | 6 to 12 deg |
+| P2, centroid constraint only | 9k to 28k | 9 to 17 deg (does not improve with refinement) |
+| P2, 4-point quadrature | 9k to 28k | 4 to 8 deg |
+
+P2 with quadrature is better than linear DFI at the same or fewer unknowns. The hinge-zone error nearly disappears.
+All values are from single runs. The model chain varies by about 3 deg between runs, so compare only inside one run.
+
+## What the prototype does
+File: `debug_s0/scripts/p2_fold.py`, class `P2FoldInterpolator(P2Interpolator)`.
+It uses the loop_interpolation package from `~/dev/LoopStructural` (version 1.8.x).
+
+1. Put sample points inside every tet. The best rule is the degree-2 rule, 4 points with barycentric coordinates
+   (a, a, a, b), a = 0.1381966011250105, b = 0.5854101966249685.
+   The centroid alone is not enough: the P2 gradient is linear inside a tet, so one point does not control it.
+   Using all 10 nodes is slightly worse (shared vertices get too much weight).
+2. Get the shape-function derivatives with `support.evaluate_shape_derivatives(points, elements=idx)`. Shape: `(n, 3, 10)`.
+3. Add three constraint families, as in `DiscreteFoldInterpolator`:
+   - gradient . fold_direction = 0 (w = 10)
+   - gradient . fold_axis = 0 (w = 10)
+   - gradient . dgz = fold_norm (w = 1, default fold_norm = -1 in the new version)
+4. Row = `einsum(dN, vector) * tet_volume`, right side = `target * tet_volume`.
+   Give each sample weight `w / n_samples`, so the total weight per tet does not depend on the rule.
+5. Regularisation: leave it off or use the P2 default. Fold-frame versions did not help (see below).
+
+The prototype is attached to the builder by swapping `builder.interpolator` and `feature.interpolator` after
+`create_and_add_folded_fold_frame`. That is a workaround. The real change should go in the library.
+
+## Suggested library work
+- Add a class (for example `P2FoldInterpolator`) in `packages/loop_interpolation/src/loop_interpolation/`.
+  Use `_discrete_fold_interpolator.py` (shared setup in `_fold_setup.py`) and `_fd_fold_interpolator.py` as templates.
+- Register a new `InterpolatorType`, string key, `interpolator_map` entry and `support_interpolator_map` entry
+  (support = `P2UnstructuredTetMesh`, built from `origin`, `step_vector`, `nsteps_cells`) in `loop_interpolation/__init__.py`.
+- Stop forcing `"DFI"`: `LoopStructural/modelling/core/_model_feature_factory.py` (about line 248),
+  `geological_model.py` (about line 784), `modelling/features/builders/_folded_feature_builder.py` (about line 75).
+- Make `nelements` clear for P2. 1e4 gives 56k tets and 86k unknowns. 1e5 gives 526k tets and 765k unknowns (solve about 2.5 min).
+- Expose the sample rule (`centroid` or `gauss4`) in `fold_weights`. Default `gauss4`.
+- Return NaN, not `LinAlgError ("Singular matrix")`, for points outside the mesh in `evaluate_shape_derivatives`.
+- Add tests (see below) and a short doc note.
+
+## Tests to add
+- The quadrature rule integrates a quadratic exactly on a reference tet.
+- Constant fold angle, no orientation data: median angle to the expected direction below a fixed limit.
+- At equal unknowns, P2 error is not worse than linear DFI on a synthetic fold.
+- Fold norm sign: use `resolve_fold_norm_target` as `DiscreteFoldInterpolator` does.
+
+## Findings to keep in mind
+- Constraint rows are normalised by `add_constraints_to_least_squares`. The size of `fold_direction` does not matter
+  (the "sum, not norm" scaling in old versions is harmless). Weights `w` are applied after the normalisation.
+- The three fold constraints fix the gradient completely. Regularisation weights 0.1 and 0.01 change nothing.
+  Large weights (about 100 to 1000) make the fold fit worse. Fold-frame face-jump and curvature terms gave
+  at most a small gain over the isotropic terms. Do not add them for now.
+- The norm constraint is almost never satisfied, and the field scale is arbitrary (mean gradient 0.1 to 0.5).
+  The direction metric is not affected, but a wavelength in S2 units is not portable between runs.
+- `fold_norm` default changed from +1 (old) to -1 (new). With +1, the old version gave a field with
+  opposite polarity to the data normals (80% of S2 data flipped).
+- `buffer` is applied twice (model box, then support). `buffer=0.2` gives a mesh about 12 times the box volume,
+  and only 7% of the elements are inside the box. `buffer=0.0` roughly halves the error at equal element count.
+  With `buffer=0.0` the field is valid only inside the box and boundary points return NaN.
+- Builds are lazy. Evaluate a feature before you read `fold` attributes.
+- Some steps use random shuffles. Do not compare numbers between separate runs.
+
+## Separate issue: rotation-angle profile fit (not P2)
+- The fit changes between runs. 2 of 8 rebuilds gave a degenerate profile (angle 0 everywhere).
+- `limb_wl` is only the optimiser start. Values 2 to 5 stay in a poor local minimum
+  (misfit 40 to 48 deg). Values 8 to 50 and auto converge to the same fit (misfit 11 deg).
+- Full P2 run (S2 1e4, S0 1e5, 765k unknowns), S0 `limb_wl=5` against `limb_wl=None` (auto):
+  profile misfit 36 deg against 12.7 deg; S0 median error 11.8 against 6.2 deg; share above 30 deg 19.5% against 0.9%;
+  roughness 44 against 30 deg per 10 m; data match 0.77 (2 of 7 flipped) against 0.92 (none flipped).
+  Plots: `debug_s0/images/p2_s2_s0_rotation_profiles.png` and `..._wlauto.png`.
+- S2 observations contradict each other (9 deg and 69 deg only 0.7 S3 units apart), so its profile fits poorly.
+- Always plot the profile with its observations: `debug_s0/scripts/p2_s2_s0.py` does this.
+
+## Reproduce
+- Environment: a venv with numpy, pandas, scipy, scikit-image, scikit-learn, pydantic, pyyaml, pyvista, dill, matplotlib
+  (and trame, trame-vtk, trame-vuetify, trame-pyvista, nest_asyncio2 for HTML export).
+- Set `PYTHONPATH` to `~/dev/LoopStructural;~/dev/LoopStructural/packages/loop_common/src;~/dev/LoopStructural/packages/loop_interpolation/src`.
+- Scripts in `debug_s0/scripts/`:
+  - `p2_fold.py` is the prototype class and the comparison runs.
+  - `s0_p2.py` compares linear DFI and P2 on S0.
+  - `p2_s2_s0.py` is S2 (P2, 1e4) and S0 (P2, 1e5) with profile plots (`--wl=auto` for the automatic wavelength).
+  - `s0_wavelength_sweep.py` is the `limb_wl` sweep.
+  - `compat_chain.py` runs the same chain in the old (1.6.27) and new (1.8.x) versions.
+- Results: `debug_s0/newversion/*.csv` and `*.log`, images in `debug_s0/images/`, interactive pages in `debug_s0/html/`.
+- Data: `Data/Processed/TrimDataS4.csv`, `S3`, `S2`, `S0_LA2`. The model box is origin (238275, 5965000, -50), maximum (238500, 5965250, 50).
+
+## Caveats
+- Single runs. The metric measures agreement with the fold constraints, not geological quality.
+- The prototype has no anisotropic fold regularisation, and it was tested only on this one dataset (5 S2 and 7 S0 points).

agent/README.md

@@ -0,0 +1,18 @@
+# agent/
+
+Working notes for `loop_interpolation`. Plans and reviews, not user documentation.
+
+| File | What it is |
+|---|---|
+| `REVIEW_interpolation_hardening.md` | Code review of the shared interpolation path, P1, P2, FD, DFI and fold setup. Bugs, and the order of work for a matrix-free Rust port. Written from reading the code. The bugs it lists are now pinned by `packages/loop_interpolation/tests/golden` (`test_known_*`). |
+| `HANDOFF_p2_fold.md` | Current P2Fold handoff with the corrected phase 4 gate table; the original prototype narrative is marked historical. The prototype scripts are in the author's Honours project (`debug_s0/scripts/p2_fold.py`), not in this repo. The historical volume-factor and per-sample-weight rules are explicitly superseded by the current section. |
+
+Golden tests: `packages/loop_interpolation/tests/golden/README.md`.
+
+
+Phase 4: [corrected Bermagui validation](HANDOFF_p2_fold.md),
+[execution checklist](p2_fold/phase4_validation/tasks.md), and
+[usage documentation](../docs/source/user_guide/fold_modelling.rst).
+Both gauss4 gates passed at fewer unknowns. S0 retains a material error tail;
+unregularised centroid was unstable. The comparison uses a shared DFI driver,
+not an end-to-end all-P2 chain. Release checks after merge remain pending.

agent/REVIEW_interpolation_hardening.md

@@ -0,0 +1,72 @@
+# Review: interpolation code hardening and matrix-free Rust port
+
+Status: review only. No code was changed. No test was run. All findings come from reading the code.
+
+Scope: shared path (`_discrete_interpolator.py`, `_solver_strategy.py`, `_solver_pipeline.py`), P1, P2, FD, DFI, fold setup.
+Package: `packages/loop_interpolation/src/loop_interpolation/`
+Related: `HANDOFF_p2_fold.md` (the P2 fold interpolator plan).
+Update: step 1 of section 4 (golden tests) is done, see `packages/loop_interpolation/tests/golden/`. It confirmed bugs 3, 4 and 11 by running code, and found three new issues (F1 to F3), listed in its README. The "No test was run" statement above describes the original review only.
+
+## 1. Corrections to the first review of the P2 fold handoff
+
+1. The volume factor does nothing. `add_constraints_to_least_squares` divides each row of A and B by the row norm (`_discrete_interpolator.py:361-366`). A volume factor on both A and B cancels exactly. The prototype and DFI both multiply by volume. It has no effect.
+2. The rule "weight `w/n_samples` per sample" is wrong. The weight multiplies the residual (`:697-698`). The energy is therefore proportional to w². With n samples of weight w/n, the energy is w²/n. Gauss4 has 1/4 of the centroid strength. The correct split is `w*sqrt(q_i)`, where the quadrature weights q_i sum to 1. Re-run the handoff table with this fix. The balance between fold constraints and regularisation changes.
+
+## 2. Bugs found
+
+### High (wrong results, no error)
+
+| # | Where | Problem |
+|---|---|---|
+| 1 | `_discrete_interpolator.py:835` | `add_equality_block` is never called. Equality constraints have no effect. Users: `_fault_builder.py:693`, `_geological_feature_builder.py:336`. |
+| 2 | `_solver_strategy.py` | Inequality constraints work only with ADMM. CG and LSMR build `Q` and ignore it. No warning. |
+| 3 | `_solver_strategy.py:66-71` | CG returns `True` after an iteration-limit warning. The model is marked up to date. |
+| 4 | `_p2interpolator.py:122,151,166` | P2 multiplies weights into A and B before row normalisation. `gpw`, `npw`, `tpw` have no effect in P2. |
+| 5 | `_p2interpolator.py:210` | `np.all(neighbours > 3)` is wrong. Border neighbours are negative. It should be `>= 0`. Tets next to elements 0 to 3 are dropped by mistake. |
+| 6 | `_discrete_interpolator.py:357-359` | One NaN row removes the whole constraint family (warning only). DFI does not filter NaN rows. The prototype does. Check S0 and S2 logs for "Constraints contain nan". |
+| 7 | `_discrete_interpolator.py:349` | `np.tile(w, n)` gives a component-major layout. Rows are point-major. Per-point weights go to the wrong rows. No effect today because the weights are uniform. |
+
+### Medium
+
+| # | Where | Problem |
+|---|---|---|
+| 8 | `_discrete_fold_interpolator.py:192,200` | `fold_normalisation` scales B and is also the row weight. Any value other than 1 changes the target. |
+| 9 | P1, P2, FD | Weight conventions differ. FD value constraints use the per-point weight column; P1 and P2 do not. FD splits `npw` by dimension; P1 and P2 do not. Default `npw`: FD 10, P1 1. The weight is a residual multiplier, not an energy weight, and no document says so. |
+| 10 | P1 `:287`, FD `:379,513` | The caller's `vectors` and `points[:, 3:6]` are changed in place. FD hard-codes columns 3:6. |
+| 11 | `_discrete_interpolator.py:1337` | `update()` calls `setup_interpolator()` with no arguments. `fold_weights`, `regularisation`, `directional_regularisation` are lost. A re-solve can differ from the first solve. |
+| 12 | `_pre_solve` | `c` is set to NaN. A failed solve leaves it NaN with `up_to_date=False`. Every evaluation retries the solve. |
+| 13 | `_post_solve` | Commented out. All constraint matrices stay in memory after the solve. This matters at 500k+ tets. |
+| 14 | `:376-380` | The duplicate-name code calls `int()` on the second `_` token. A name such as `gradient_orthogonal` crashes on its second use. |
+
+### Low
+
+15. P2 `wtfunc` replaces `w` and receives element indices, not points (`_p2interpolator.py:274-275`).
+16. P2 steepness: off-diagonal Hessian terms count once. A Frobenius norm counts them twice.
+17. `evaluate_shape_derivatives` returns 3 values in P1 and 2 in P2.
+18. `ridge_factor=1e-8` becomes 1e-16 in the normal equations. This is below float precision. The term has no effect.
+
+## 3. Problems for a matrix-free Rust port
+
+1. Row normalisation defines the weights. Each regularisation row is divided by its own norm. The result is not a mesh integral. It does not change with element size. The value-constraint weight changes by up to 2x across one element (the shape-function norm is between 0.5 and 1). Decide the intended energy first, for example the integral over each element. Implement that in Rust. Do not copy the current behaviour.
+2. FD already has a partial matrix-free path (`_assemble_operator`, `_build_fused_cg_regularisation_operator`). It covers regularisation only. Data constraints stay as a sparse matrix (`AtA_data`). P1 and P2 have none. The fused FD kernel is not exact near the boundary, so the code patches the 2-cell shell. Reuse the design (stencil kernels, Jacobi preconditioner). Do not reuse the patch.
+3. Both paths solve the normal equations `AtA x = Atb`. This squares the condition number. Use one of:
+   - CG on the normal equations with a Jacobi preconditioner. Column scaling does this today. The diagonal of AtA is easy to compute matrix-free.
+   - LSQR or LSMR on the operator (no squaring).
+4. Tet meshes are the hard part. The operator is a sum over elements and faces. Edge jumps need face-neighbour tables. Define each operator as `apply(x)` and `apply_transpose(y)`. Test every pair with the dot-product check: <Ax, y> = <x, ATy>.
+5. Point location and sampling are in Python. `evaluate_shape_derivatives` raises `LinAlgError` for points outside the mesh. Rust must return NaN or a mask. For inside points, it must give the same answer as Python.
+6. DFI uses `rng.shuffle` and `[::step]`. Runs are not repeatable. Rust needs a deterministic rule.
+
+## 4. Proposed order of work
+
+1. Pin the current behaviour. Write golden tests (small meshes, fixed seeds) for P1, P2, FD and DFI. Save the matrix for each constraint family. These tests become the Rust parity tests.
+2. Fix bugs 1 to 7 and 14. Fix each in its own commit. Write a failing test first. Keep the golden files to show exactly what changes.
+3. Write the operator specification: one page per family. State the energy, the weight convention (w or w squared), and the row-normalisation rule. Settle items 9 and 10 here.
+4. Add a matrix-free `LinearOperator` for P1 and P2 in Python, with the transpose test. Check it against the assembled matrix to about 1e-10.
+5. Port to Rust behind the same operator interface. Use the assembled matrix as the oracle.
+6. Build the P2 fold interpolator after step 4. It then needs only a quadrature loop. This avoids writing it twice.
+
+## 5. Not checked
+
+- Support code was not read in detail: shape functions, `element_size`, `get_quadrature_points`, P2 node ordering.
+- `constraints.py` and the `get_*_constraints` code were not read.
+- No test was run. Items 5, 6 and 7 come from reading the code. Confirm item 6 from the logs.