diff --git a/.gitignore b/.gitignore index 6bc4e27a8..77555341d 100644 --- a/.gitignore +++ b/.gitignore @@ -137,4 +137,6 @@ conda/channeldata.json conda/noarch conda/win-64 .DS_Store -sandbox/* \ No newline at end of file +sandbox/* +# Local validation outputs +/validation/ diff --git a/LoopStructural/modelling/core/_model_feature_factory.py b/LoopStructural/modelling/core/_model_feature_factory.py index 95d75e043..11bafc939 100644 --- a/LoopStructural/modelling/core/_model_feature_factory.py +++ b/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 diff --git a/agent/HANDOFF_p2_fold.md b/agent/HANDOFF_p2_fold.md new file mode 100644 index 000000000..596c15b26 --- /dev/null +++ b/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). diff --git a/agent/README.md b/agent/README.md new file mode 100644 index 000000000..cb8c0def6 --- /dev/null +++ b/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. diff --git a/agent/REVIEW_interpolation_hardening.md b/agent/REVIEW_interpolation_hardening.md new file mode 100644 index 000000000..2b6a50351 --- /dev/null +++ b/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: = . +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. diff --git a/agent/p2_fold/README.md b/agent/p2_fold/README.md new file mode 100644 index 000000000..7c3100584 --- /dev/null +++ b/agent/p2_fold/README.md @@ -0,0 +1,49 @@ +# P2 fold interpolation: work plan + +Sources: `agent/HANDOFF_p2_fold.md`, `agent/REVIEW_interpolation_hardening.md`, `packages/loop_interpolation/tests/golden/README.md`. + +## Decisions (settled) + +1. Build `P2FoldInterpolator` on the assembled-matrix path. Do not wait for the matrix-free operator. +2. Add a new type string `"P2Fold"`. Plain `"P2"` does not change. +3. The double `buffer` fix is not part of this work. It is a separate `fix:` commit after Phase 4 (see `phase5_buffer_followup`). + +## Phases + +| Phase | Folder | Goal | +|---|---|---| +| 0 | `phase0_prereqs` | Fix the P2 and shared bugs that the fold class would inherit | +| 1 | `phase1_class` | Write `P2FoldInterpolator` and the quadrature rule | +| 2 | `phase2_wiring` | Register `"P2Fold"` and stop forcing `"DFI"` | +| 3 | `phase3_tests` | Unit, golden and accuracy tests | +| 4 | `phase4_validation` | Re-run Bermagui, write docs, release | +| 5 | `phase5_buffer_followup` | Single-point buffer fix (after Phase 4) | + +Each folder has `plan.md` (why and how) and `tasks.md` (checklist). + +## Rules for all phases + +- One fix per commit. Write the failing test first. +- Use `fix:` and `feat:` prefixes. Release-please sets versions. +- Do not add a Co-Authored-By line to commits (standing user preference). +- Run `python -m pytest tests/golden -q` from `packages/loop_interpolation` before each commit. Review every `--update-golden` diff by hand. The golden files are the Rust parity oracle. +- Single runs vary by about 3 deg. Compare numbers only inside one run, with fixed seeds. + +## Open issues outside this work + +- F2: `TetMesh` returns NaN on the lower x, y, z faces. +- F3: Fold amplitude is damped (gradient ratio 0.35 to 0.66). Decide the intended energy in the Rust operator specification. +- The rotation-angle profile fit is not repeatable. `limb_wl` below 5 stays in a poor local minimum. +- Matrix-free P1/P2 operator and Rust port (review section 4, steps 3 to 5). + + +## Status after phase 4 validation (20 October 2026) + +Phases 0-3 are on master. Phase 4 measurements and the decision gate are complete: +S2 gauss4 4.69 deg / 9,711 unknowns versus DFI 6.15 deg / 12,992; S0 gauss4 +11.14 deg / 85,673 versus DFI 20.63 deg / 112,560. See the +[summary and limitations](../HANDOFF_p2_fold.md). +The old prototype weighting and "nelements is tets" statements are superseded. +The API note and bundled-data example are in [PR #349](https://github.com/Loop3D/LoopStructural/pull/349). Post-merge +release scoping and Phase 5 remain pending; this phase does not change buffer +behaviour or manually bump versions. diff --git a/agent/p2_fold/phase0_prereqs/plan.md b/agent/p2_fold/phase0_prereqs/plan.md new file mode 100644 index 000000000..f61a081bc --- /dev/null +++ b/agent/p2_fold/phase0_prereqs/plan.md @@ -0,0 +1,36 @@ +# Phase 0: prerequisites + +## Goal +Remove the P2 and shared bugs that `P2FoldInterpolator` would inherit. After this phase, P2 weights work, P2 solves converge, and fold weights survive `update()`. + +## Why first +- Bug 4: P2 ignores `gpw`, `npw`, `tpw`. A P2 fold class would ignore its fold weights in the same way. +- Bug 3: The default CG solve on P2 can be wrong by up to 0.64 in a unit box and still return `True`. +- Bug 11: `update()` drops `fold_weights`. Fold classes are the main users. +- F1: The P2 support has no `origin`, `step_vector` or `nsteps_cells`. Registration in Phase 2 needs them. + +## Items +All paths are in `packages/loop_interpolation/src/loop_interpolation/`. + +| Item | Where | Fix | +|---|---|---| +| F1 | `P2UnstructuredTetMesh` constructor | Store `origin`, `step_vector`, `nsteps_cells`. Remove the workaround in `tests/golden/engines.py`. | +| Bug 4 | `_p2interpolator.py:122,151,166` | Apply weights after row normalisation, as P1 does. Turn `test_known_p2_normal_weight_has_effect` from strict xfail into a normal test. | +| Bug 5 | `_p2interpolator.py:210` | Change `np.all(neighbours > 3)` to `>= 0`. Border neighbours are negative. | +| Bug 3 | `_solver_strategy.py:66-71` | Return `False` when CG hits the iteration limit. Keep the warning. Turn `test_known_default_solve_is_converged` into a normal test. | +| Bug 11 | `_discrete_interpolator.py:1337` | `update()` must pass the stored `fold_weights`, `regularisation` and `directional_regularisation` to `setup_interpolator`. Turn `test_known_update_keeps_fold_weights` into a normal test. | +| Bug 6 | `_discrete_interpolator.py:357-359` and `_discrete_fold_interpolator.py` | Drop NaN rows only, not the whole family. Log how many rows you dropped. Apply the same filter in the shared fold setup. | + +## Order and risk +1. F1 first. Other fixes need a working P2 test harness. +2. Bug 5, then bug 4, then bug 3. Each one changes P2 golden files. Review each diff alone. +3. Bug 11 and bug 6 change DFI and FDFold results. Review those golden diffs with care. +4. Bug 3 may make some existing solves report failure. Check callers of the return value of the solver. + +## Done when +- The three `test_known_*` tests named above pass without xfail. +- Golden diffs are reviewed and accepted, each in its own commit. +- `python -m pytest tests/golden -q` is green. + +## Not in this phase +Bugs 1, 2, 7, 8, 9, 10, 12, 13, 14 and the low items. They do not block P2Fold. Bug 7 has no effect while weights are uniform. Re-check it if Phase 1 uses per-sample weights. diff --git a/agent/p2_fold/phase0_prereqs/tasks.md b/agent/p2_fold/phase0_prereqs/tasks.md new file mode 100644 index 000000000..9f931fd79 --- /dev/null +++ b/agent/p2_fold/phase0_prereqs/tasks.md @@ -0,0 +1,32 @@ +# Phase 0 tasks + +Run tests from `packages/loop_interpolation`. One commit per task. + +- [ ] 0.1 F1: store `origin`, `step_vector`, `nsteps_cells` on `P2UnstructuredTetMesh` + - [ ] Write a failing test that reads the three attributes + - [ ] Fix the constructor + - [ ] Remove the manual attribute setting in `tests/golden/engines.py` + - [ ] Confirm the diagnostics report `inside()` works for P2 +- [ ] 0.2 Bug 5: change `neighbours > 3` to `neighbours >= 0` in `_p2interpolator.py:210` + - [ ] Write a failing test: tets next to elements 0 to 3 are kept + - [ ] Fix, then review the P2 golden diff +- [ ] 0.3 Bug 4: apply P2 weights after row normalisation + - [ ] Change `test_known_p2_normal_weight_has_effect` into a normal test (remove xfail) + - [ ] Fix `_p2interpolator.py:122,151,166` + - [ ] Check that `npw` 1 to 100 now moves the field by more than 1e-4 + - [ ] Review the P2 golden diff +- [ ] 0.4 Bug 3: CG iteration limit returns `False` + - [ ] Change `test_known_default_solve_is_converged` into a normal test + - [ ] Fix `_solver_strategy.py:66-71` + - [ ] Search for callers that use the return value. Check that `up_to_date` stays correct + - [ ] Decide the P2 default tolerance or iteration count. Record it in the commit message +- [ ] 0.5 Bug 11: `update()` keeps `fold_weights`, `regularisation`, `directional_regularisation` + - [ ] Change `test_known_update_keeps_fold_weights` into a normal test + - [ ] Fix `_discrete_interpolator.py:1337` + - [ ] Review DFI and FDFold golden diffs +- [ ] 0.6 Bug 6: filter NaN rows, not the whole family + - [ ] Write a failing test with one NaN row in a constraint family + - [ ] Fix `_discrete_interpolator.py:357-359` + - [ ] Add the same filter to the shared fold setup (`_fold_setup.py`) + - [ ] Log the count of dropped rows +- [ ] 0.7 Run the full golden suite. Update the `Known bugs` table in `tests/golden/README.md` diff --git a/agent/p2_fold/phase1_class/plan.md b/agent/p2_fold/phase1_class/plan.md new file mode 100644 index 000000000..39c1d7695 --- /dev/null +++ b/agent/p2_fold/phase1_class/plan.md @@ -0,0 +1,41 @@ +# Phase 1: `P2FoldInterpolator` + +## Goal +A working `P2FoldInterpolator(P2Interpolator)` with a quadrature sample rule. It is not yet registered in the factory (Phase 2). + +New file: `packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py`. +Templates: `_discrete_fold_interpolator.py`, `_fold_setup.py`, `_fd_fold_interpolator.py`. + +## Corrections to the handoff (from the review, section 1) +- Do not multiply by tet volume. Row normalisation cancels it. +- Do not use weight `w / n_samples`. The weight multiplies the residual, so energy goes with w squared. Use `w * sqrt(q_i)`, where the quadrature weights `q_i` sum to 1. Then the total energy per tet does not depend on the rule. + +## Design + +1. **Quadrature.** Add two rules on the reference tet: + - `centroid`: one point (1/4, 1/4, 1/4, 1/4), weight 1. + - `gauss4`: four points with barycentric coordinates (a, a, a, b) and permutations. a = 0.1381966011250105, b = 0.5854101966249685. Each weight is 1/4. + Check `get_quadrature_points` in the P2 support first. Reuse it if it matches. The review did not read it. +2. **Sample points.** Map the rule to every tet (or to every tet in the active region). Keep the tet index of each point. +3. **Shape derivatives.** `support.evaluate_shape_derivatives(points, elements=idx)`. Shape `(n, 3, 10)`. The P2 gradient is linear in a tet. One point does not control it. This is why `gauss4` beats `centroid`. +4. **Three constraint families**, as in `DiscreteFoldInterpolator`: + - gradient . fold_direction = 0 (default weight 10) + - gradient . fold_axis = 0 (default weight 10) + - gradient . dgz = fold_norm (default weight 1) + Get the shared vectors from `_fold_setup.py`. Get the norm target from `resolve_fold_norm_target`, so the sign matches DFI. +5. **Weights.** Row = `einsum(dN, vector)`. Target = `fold_norm` or 0. Weight = `w * sqrt(q_i)`. Pass the weight through `add_constraints_to_least_squares` so it is applied after normalisation. +6. **Regularisation.** Use the P2 default. Add no fold-specific terms (the handoff found no gain). +7. **Options.** `fold_weights["sample_rule"]` is `"gauss4"` (default) or `"centroid"`. +8. **NaN rows.** Filter with the Phase 0 code from bug 6. +9. **`update()`.** Must keep `fold_weights` (Phase 0, bug 11). +10. **Outside points.** `evaluate_shape_derivatives` returns NaN for points outside the mesh, not `LinAlgError`. Make the P1 and P2 return signatures match (review item 17). + +## Risks +- The handoff numbers used the wrong weight rule. Results will differ. Phase 4 re-runs them. +- Per-point weights meet bug 7 (`np.tile` layout). If a per-sample weight vector goes through `add_constraints_to_least_squares`, check the row order. +- Memory: 4 samples per tet at 526k tets is 2.1M rows per family. Check memory at `nelements=1e5`. + +## Done when +- The class builds, solves and evaluates on the `fold_sparse` golden case when you construct it by hand. +- The quadrature unit test passes. +- Both sample rules work, and the weight option changes the field. diff --git a/agent/p2_fold/phase1_class/tasks.md b/agent/p2_fold/phase1_class/tasks.md new file mode 100644 index 000000000..5a71137ef --- /dev/null +++ b/agent/p2_fold/phase1_class/tasks.md @@ -0,0 +1,27 @@ +# Phase 1 tasks + +Needs Phase 0 done. + +- [ ] 1.1 Read the P2 support: `get_quadrature_points`, `evaluate_shape_derivatives`, node ordering, `element_size` + - [ ] Note what already exists. Write findings in a comment at the top of the new file +- [ ] 1.2 Quadrature + - [ ] Add the `centroid` and `gauss4` rules (reuse the support function if it matches) + - [ ] Test: the rule integrates all quadratics exactly on the reference tet + - [ ] Test: weights sum to 1 +- [ ] 1.3 `evaluate_shape_derivatives` for outside points + - [ ] Write a failing test: a point outside the mesh gives NaN + - [ ] Fix. Remove the `LinAlgError` + - [ ] Make P1 and P2 return the same number of values (review item 17). Update callers +- [ ] 1.4 Create `_p2_fold_interpolator.py` with class `P2FoldInterpolator(P2Interpolator)` + - [ ] Constructor takes the support and keeps `fold` + - [ ] Add `fold_weights` defaults: 10, 10, 1, `sample_rule="gauss4"` +- [ ] 1.5 `setup_interpolator` for the three constraint families + - [ ] Sample points per tet, with the tet index + - [ ] Row = `einsum(dN, vector)`, no volume factor + - [ ] Weight `w * sqrt(q_i)`, applied after normalisation + - [ ] Norm target from `resolve_fold_norm_target` + - [ ] Filter NaN rows +- [ ] 1.6 Keep the settings on `update()` (check against the Phase 0 fix) +- [ ] 1.7 Check order of rows when per-sample weights are used (bug 7). Add a test with unequal weights +- [ ] 1.8 Smoke test by hand on the `fold_sparse` golden case. Compare `centroid` and `gauss4` +- [ ] 1.9 Memory check at about 5e5 tets. Record the peak memory and solve time diff --git a/agent/p2_fold/phase2_wiring/plan.md b/agent/p2_fold/phase2_wiring/plan.md new file mode 100644 index 000000000..155767ba6 --- /dev/null +++ b/agent/p2_fold/phase2_wiring/plan.md @@ -0,0 +1,46 @@ +# Phase 2: registration and wiring + +## Goal +`create_and_add_folded_fold_frame(..., interpolatortype="P2Fold")` builds a model with the new class. `"DFI"` stays the default for folds. + +## Decision +New type string `"P2Fold"`. Plain `"P2"` does not change. + +## Items + +### 1. Registration +File: `packages/loop_interpolation/src/loop_interpolation/__init__.py`, with `_interpolatortype.py`. +- Add `InterpolatorType.PIECEWISE_QUADRATIC_FOLD` (name can change; match the existing naming). +- Add the string `"P2Fold"` to the string map (near line 97). +- Add the class to `interpolator_map` (near line 112). +- Add `P2UnstructuredTetMesh` to `support_interpolator_map` for dimension 3. +- Export `P2FoldInterpolator` in `__all__`. +- Check `_normalise_interpolator_type` and `from_dict`. A saved model with this type must load. + +### 2. Stop forcing DFI +Three places: +- `LoopStructural/modelling/core/_model_feature_factory.py:248-250`: accept `DFI` and `P2Fold`. Warn only for other types. +- `LoopStructural/modelling/core/geological_model.py:784`: the default stays `"DFI"`. Update the docstring. +- `LoopStructural/modelling/features/builders/_folded_feature_builder.py:75`: `from_feature_builder` hard-codes `'DFI'`. Take the type from the source builder if it is a fold type. +- Check `_model_feature_factory.py:346`, which also names `"DFI"`. + +### 3. Fold builder +- Check that `FoldedFeatureBuilder` does not use DFI-only attributes (for example `interpolator.fold` and `fold_weights` pass-through). +- Check that the support type works with the fold frame, the region and the fault handling in `_add_faults`. + +### 4. Element count +- `nelements` means tets for P2. 1e4 gives 56k tets and 86k unknowns. 1e5 gives 526k tets and 765k unknowns (solve about 2.5 min). +- Log the number of tets and unknowns when the support is built. +- Describe this in the docstring of `create_and_add_folded_fold_frame`. + +### 5. Buffer +Do not change the buffer behaviour here (see `phase5_buffer_followup`). In docs and examples, say that the mesh is larger than the model box. + +## Risks +- Serialisation (`to_dict` and `from_dict`) may miss the new type. +- Other code may test `interpolatortype == "DFI"`. Search for it. +- Code that checks `isinstance(interpolator, DiscreteFoldInterpolator)` will miss the new class. Search for it. + +## Done when +- A model script with `interpolatortype="P2Fold"` builds, evaluates and saves/loads. +- `"DFI"` models give the same results as before. diff --git a/agent/p2_fold/phase2_wiring/tasks.md b/agent/p2_fold/phase2_wiring/tasks.md new file mode 100644 index 000000000..0d2c9518d --- /dev/null +++ b/agent/p2_fold/phase2_wiring/tasks.md @@ -0,0 +1,25 @@ +# Phase 2 tasks + +Needs Phase 1 done. + +- [x] 2.1 Search for hard-coded fold types + - [x] `grep -rn "DFI\|DISCRETE_FOLD\|DiscreteFoldInterpolator" LoopStructural packages` + - [x] List every place that must accept the new type +- [x] 2.2 Register the type + - [x] Add the `InterpolatorType` member + - [x] Add the string `"P2Fold"` + - [x] Add `interpolator_map` and `support_interpolator_map` entries + - [x] Export the class + - [x] Test: `InterpolatorFactory.create_interpolator("P2Fold", bb, nelements=...)` returns the class with a P2 support + - [x] Test: `to_dict` and `from_dict` round trip +- [x] 2.3 `_model_feature_factory.py:248-250`: allow `DFI` and `P2Fold`. Warn for other types +- [x] 2.4 `_model_feature_factory.py:346`: check the second `"DFI"` default +- [x] 2.5 `geological_model.py:784`: update the docstring. Keep the default +- [x] 2.6 `_folded_feature_builder.py:75`: take the type from the source builder +- [x] 2.7 Check `FoldedFeatureBuilder` with the new class + - [x] `interpolator.fold` is set + - [x] `fold_weights` pass through + - [x] Faults and regions work (`_add_faults`) +- [x] 2.8 Log tets and unknowns when the P2 support is built. Document `nelements` for P2 +- [x] 2.9 End-to-end script: small synthetic model, `create_and_add_folded_fold_frame` with `P2Fold`, evaluate on a grid, save and load +- [x] 2.10 Run the full test suite for `LoopStructural` and `packages/loop_interpolation`. Confirm no DFI change diff --git a/agent/p2_fold/phase3_tests/plan.md b/agent/p2_fold/phase3_tests/plan.md new file mode 100644 index 000000000..d44363871 --- /dev/null +++ b/agent/p2_fold/phase3_tests/plan.md @@ -0,0 +1,43 @@ +# Phase 3: tests + +## Goal +Pin the behaviour of `P2Fold` so that later work (the matrix-free operator and the Rust port) has an oracle. + +Tests live in `packages/loop_interpolation/tests/` and in the golden suite `tests/golden/`. + +## Tests + +### Unit +1. Quadrature: exact for all quadratics on the reference tet. Weights sum to 1. (Written in Phase 1.) +2. Outside points give NaN. (Written in Phase 1.) +3. Weight rule: for the same fold data, the total energy of the fold family is the same for `centroid` and `gauss4` (within a small tolerance), because of `w * sqrt(q_i)`. +4. Weights have an effect: change `fold_weights` from 10 to 50. The field moves by more than a fixed limit. This catches bug 4 returning. +5. Fold norm sign: the sign matches `DiscreteFoldInterpolator` for the same fold event. Use `resolve_fold_norm_target`. +6. `update()` keeps `fold_weights` and `sample_rule`. +7. NaN rows: a fold event that gives NaN at some points still builds. The count of dropped rows is logged. +8. Dot-product test, only if you build the constraints as a `LinearOperator`: ` = `. + +### Golden +1. Add `P2Fold` to `tests/golden/engines.py`. +2. Run it on `fold_sparse` and `refolded_sparse`. Save the golden files (field and constraint fingerprints). +3. Set bounds in `test_golden.py::BOUNDS` about 1.5 times above what the code gives (see golden README). +4. Add `P2Fold` to the figures in `tests.golden.run`. + +### Accuracy +1. On `fold_sparse`, the median gradient direction error of `P2Fold` is not worse than `DFI` at equal unknowns. Use a fixed seed and an explicit buffer for both. +2. Constant fold angle, no orientation data: the median angle to the expected direction is below a fixed limit. +3. `gauss4` is better than `centroid`. Check this on the golden cases. If it is not, record it and do not make it a hard test. + +### Integration +1. `create_and_add_folded_fold_frame(interpolatortype="P2Fold")` does not fall back to DFI (no warning, `type(interpolator)` is `P2FoldInterpolator`). +2. Save and load keeps the type. + +## Notes +- Mesh sizes in tests must be small. The P2 golden mesh is 8^3 cells (4.4k unknowns). Keep the whole suite near 20 s plus the new cases. +- Do not compare numbers between separate runs without a fixed seed. +- The refolded case has a 2 deg floor from the fold event model (see golden README). Do not set a bound below it. + +## Done when +- All new tests pass on a clean checkout. +- Golden diff is reviewed. +- The golden README lists P2 fold as covered. diff --git a/agent/p2_fold/phase3_tests/tasks.md b/agent/p2_fold/phase3_tests/tasks.md new file mode 100644 index 000000000..d15fa1457 --- /dev/null +++ b/agent/p2_fold/phase3_tests/tasks.md @@ -0,0 +1,32 @@ +# Phase 3 tasks + +Needs Phase 2 done. Run from `packages/loop_interpolation`. + +## Unit +- [ ] 3.1 Weight rule: fold-family energy is the same for `centroid` and `gauss4` +- [ ] 3.2 Weights have an effect: `fold_weights` 10 to 50 moves the field +- [ ] 3.3 Fold norm sign matches `DiscreteFoldInterpolator` +- [ ] 3.4 `update()` keeps `fold_weights` and `sample_rule` +- [ ] 3.5 NaN fold rows do not stop the build. Dropped rows are logged +- [ ] 3.6 (Optional) Dot-product test if constraints become a `LinearOperator` + +## Golden +- [ ] 3.7 Add `P2Fold` to `tests/golden/engines.py` +- [ ] 3.8 Run `python -m tests.golden.run --case fold_sparse --engine P2Fold`. Read the figure +- [ ] 3.9 Run the same for `refolded_sparse`. Read the figure +- [ ] 3.10 Set `BOUNDS` in `test_golden.py` (about 1.5 times the observed values) +- [ ] 3.11 `python -m pytest tests/golden --update-golden`. Review the new files only +- [ ] 3.12 Update `tests/golden/README.md`: engine list, "Not covered yet", findings + +## Accuracy +- [ ] 3.13 `P2Fold` is not worse than `DFI` on `fold_sparse` at equal unknowns (fixed seed, same explicit buffer) +- [ ] 3.14 Constant fold angle without orientation data: median angle below a fixed limit +- [ ] 3.15 Compare `gauss4` and `centroid`. Record the result in the golden README + +## Integration +- [ ] 3.16 Factory test: `P2Fold` does not fall back to DFI +- [ ] 3.17 Save and load test keeps the type + +## Finish +- [ ] 3.18 Whole suite runs in a time you accept. Record the time +- [ ] 3.19 Run `python -m pytest` in `LoopStructural` root. No new failures diff --git a/agent/p2_fold/phase4_validation/plan.md b/agent/p2_fold/phase4_validation/plan.md new file mode 100644 index 000000000..b8562fbd4 --- /dev/null +++ b/agent/p2_fold/phase4_validation/plan.md @@ -0,0 +1,45 @@ +# Phase 4: validation, docs, release + +## Goal +Confirm that the benefit seen in the prototype is real after the weight corrections. Then document and release. + +## Why re-run +The handoff table used two wrong rules: a useless volume factor, and weight `w / n_samples` (energy `w^2 / n`). Gauss4 had 1/4 of the centroid strength. The balance between fold constraints and regularisation was different. The numbers (4 to 8 deg for gauss4, 6 to 12 deg for linear DFI) may change. + +## Validation +Data: Bermagui S2 and S0 (`Data/Processed/TrimDataS4.csv`, `S3`, `S2`, `S0_LA2`). Box: origin (238275, 5965000, -50), maximum (238500, 5965250, 50). The scripts are in the author's Honours project (`debug_s0/scripts/`), not in this repo. + +Runs: +1. Linear DFI. +2. `P2Fold` with `centroid`. +3. `P2Fold` with `gauss4`. + +Rules for a fair comparison: +- Same data, same fixed seeds, same explicit `buffer` for all engines. The double buffer issue changes mesh size. +- Compare only inside one run. The model chain varies by about 3 deg. +- Use S0 with `limb_wl=None` (auto). Values 2 to 5 give a poor local minimum. +- Metric: angle between the unit gradient and `cross(fold_direction, fold_axis)`. Report median and share above 30 deg. +- Record unknowns, tets, solve time and peak memory. + +## Decision gate +- If `gauss4` is better than DFI at equal or fewer unknowns: continue. +- If it is not: stop. Write the new numbers in the handoff. Ask the user before you document `P2Fold` as the better choice. Keep it as an experimental option. + +## Docs +Short note, with an example: +- When to use `P2Fold`. +- Mesh size: `nelements` is tets. 1e4 gives 86k unknowns. 1e5 gives 765k unknowns and about 2.5 min. +- `sample_rule` option. +- The field scale is arbitrary (the norm constraint is almost never met). Direction is reliable. A wavelength in field units is not portable between runs. +- Open issues: buffer, F2, F3, rotation-profile fit. + +## Release +- Commits use `feat:` for the class and `fix:` for Phase 0 items. Release-please sets the version. +- No Co-Authored-By line. +- Update the handoff file with the corrected table, or mark it as superseded. +- Update `agent/README.md`. + +## Done when +- The new table is in the handoff or in the docs. +- The decision gate passed, or the user accepted the result. +- Docs are merged. diff --git a/agent/p2_fold/phase4_validation/tasks.md b/agent/p2_fold/phase4_validation/tasks.md new file mode 100644 index 000000000..2af927291 --- /dev/null +++ b/agent/p2_fold/phase4_validation/tasks.md @@ -0,0 +1,70 @@ +# Phase 4 tasks + +Needs Phase 3 done. + +## Validation +- [x] 4.1 Get the prototype scripts from the Honours project (`debug_s0/scripts/`). Port the comparison to use the library class, not the engine swap +- [x] 4.2 Run S2 (1e4) for DFI, `P2Fold` `centroid`, `P2Fold` `gauss4`. Same seed, same explicit buffer +- [x] 4.3 Run S0 (1e5, `limb_wl=None`) for the same three engines +- [x] 4.4 Record: median error, share above 30 deg, unknowns, tets, solve time, peak memory +- [x] 4.5 Compare with the old handoff table. Explain the differences +- [x] 4.6 Decision gate: does `gauss4` beat DFI at equal or fewer unknowns? + - [x] Yes: continue + - [ ] No: stop and ask the user (not applicable: both gates passed) + +## Docs +- [x] 4.7 Write the doc note and an example script +- [x] 4.8 Update `agent/HANDOFF_p2_fold.md`: add the new table. Mark the wrong rules as fixed +- [x] 4.9 Update `agent/README.md` and `agent/p2_fold/README.md` + +## Release +- [x] 4.10 Check all commits have correct prefixes and no Co-Authored-By line +- [x] 4.11 Open a PR. Describe the Phase 0 fixes and the new class +- [ ] 4.12 After merge, check the release-please PR bumps only the right packages (see the release-please path scoping memory) +- [ ] 4.13 Start Phase 5 + + +## Execution notes - 20 October 2026 + +Worktree: `.claude/worktrees/p2-phase4`, branch `p2-phase4-validation`, based on +`origin/master` at b1ca2cd7. External project: +`C:/Users/lachl/scratch/bonnie/Honours_model23.09`. + +All ten candidate solves completed with finite inner probes. The +[handoff](../../HANDOFF_p2_fold.md) records the corrected summary and limitations. +Detailed results, input hashes and figures are generated locally by the comparison scripts. +The code uses the library factory; no engine swap or prototype class is used. +Extra smaller P2 cases make the unknown-budget decision gate meaningful. +S0 uses a fixed DFI S2 driver for every candidate; this is documented explicitly. + +The existing release PRs #348 (root), #347 (loop-common), and #344 +(loop-interpolation) were checked before merge and touch only their intended +package/version files plus manifest entries. This does not substitute for task +4.12 after the phase 4 PR merges. No manual version change or Phase 5 buffer fix +is included. Task 4.13 is intentionally conditional on phase 4 documentation +being merged, per the settled phase ordering. + + +Checks: ten external-data candidate solves; 172 interpolation tests passed +(9 skips, 2 expected failures); 9 gate/integration tests passed; bundled example +smoke-tested; Sphinx HTML built with gallery execution disabled. The 286 docs +warnings are outside the new fold page. Existing golden files are unchanged. + + +Release follow-up found during validation: root dependency ranges for loop-common +and loop-interpolation are `<0.1.0`, excluding the pending 0.1.0 releases. They +must be aligned before publishing the new root release. This is separate from +package path scoping, which currently passes. No release/version PR is merged +by this validation task. + +PR: [#349](https://github.com/Loop3D/LoopStructural/pull/349). Merge and tasks 4.12-4.13 remain pending. + + +Detailed validation artifacts are local outputs and are excluded from the PR. +The scripts default to the ignored `.pytest_cache/p2_fold` directory. Reproduce +from the repository root with workspace source packages on `PYTHONPATH`: + +```powershell +python scripts/validation/p2_fold_bermagui.py --project 'C:/Users/lachl/scratch/bonnie/Honours_model23.09' +python scripts/validation/summarise_p2_fold.py .pytest_cache/p2_fold/runtime +``` diff --git a/agent/p2_fold/phase5_buffer_followup/plan.md b/agent/p2_fold/phase5_buffer_followup/plan.md new file mode 100644 index 000000000..a4cf412c3 --- /dev/null +++ b/agent/p2_fold/phase5_buffer_followup/plan.md @@ -0,0 +1,35 @@ +# Phase 5: single-point buffer fix (after Phase 4) + +## Problem +The model box is expanded twice before the mesh is built. Found by reading the code. Not yet run. + +1. `LoopStructural/modelling/core/_model_feature_factory.py:253` passes `model.bounding_box.with_buffer(buffer)` to the builder. Default `buffer` is 0.1. +2. `LoopStructural/modelling/features/builders/_geological_feature_builder.py:54-59` calls `InterpolatorFactory.create_interpolator(..., buffer=kwargs.get("buffer", 0.2))`. +3. `packages/loop_common/src/loop_common/supports/_support_factory.py:48-49` calls `with_buffer` again. + +`buffer` is a named argument of the factory function, so it is not in `kwargs`. The user's value only sets the first expansion. The second is always 0.2. + +`with_buffer` adds `buffer * (longest box side)` to every side (`_bounding_box.py:458-459`). Defaults give 0.1 + 0.2 = 0.3 per side. With the handoff `buffer=0.2`, the total is 0.4 per side. For a flat box, the mesh volume is about 12 times the box volume, and only about 7% of the elements are inside the box. + +The same pattern is at `_model_feature_factory.py:164, 352, 500`. + +## Fix +- Apply the buffer in one place only: the support factory. +- Remove the `with_buffer` calls in `_model_feature_factory.py` (lines 164, 253, 352, 500). +- Forward the user `buffer` to the builder, so it reaches `create_interpolator`. +- Keep the total default near the present total (0.3), or every model changes. Decide the new default with the user. +- `buffer=0.0` gives a field that is valid only inside the box. Boundary points return NaN (related to F2). Document it. + +## Risks +- Every existing model changes mesh size and results. Golden files change. +- Users who pass `buffer` get a different mesh from before. +- Release note needed. This may be a breaking change for numerical results. + +## Why separate +Mixing it into the P2Fold work makes the accuracy comparison hard to read. In Phase 4 all engines use an explicit buffer to avoid the issue. + +## Done when +- One `with_buffer` call is applied per interpolator build. +- A test checks that the mesh extent equals the box plus the stated buffer. +- Golden files are updated and reviewed. +- Release note is written. diff --git a/agent/p2_fold/phase5_buffer_followup/tasks.md b/agent/p2_fold/phase5_buffer_followup/tasks.md new file mode 100644 index 000000000..1f6b76c03 --- /dev/null +++ b/agent/p2_fold/phase5_buffer_followup/tasks.md @@ -0,0 +1,17 @@ +# Phase 5 tasks + +Needs Phase 4 done. + +- [ ] 5.1 Confirm the problem by running code + - [ ] Build a model with `buffer=0.1` and read the support `origin` and `maximum` + - [ ] Compare with the box plus 0.1, 0.2 and 0.3 times the longest side + - [ ] Check the same for the code at `_model_feature_factory.py:164, 352, 500` +- [ ] 5.2 Ask the user for the new default total buffer +- [ ] 5.3 Write a failing test: mesh extent equals box plus the stated buffer +- [ ] 5.4 Remove `with_buffer` in `_model_feature_factory.py` (lines 164, 253, 352, 500) +- [ ] 5.5 Pass the user `buffer` to `GeologicalFeatureBuilder` and on to `create_interpolator` +- [ ] 5.6 Check other callers: `LoopStructural/interpolators/_api.py:24`, the fault builder (`fault_buffer`, `fault_bounding_box_buffer`), `_geological_feature_builder.py:485` +- [ ] 5.7 Review golden diffs (mesh size changes). Update only after reading them +- [ ] 5.8 Document `buffer=0.0` behaviour (NaN on the boundary) +- [ ] 5.9 Write a release note. State that mesh size and results change +- [ ] 5.10 Commit as `fix:` (or `fix!:` if the user wants a breaking-change mark) diff --git a/docs/source/user_guide/fold_modelling.rst b/docs/source/user_guide/fold_modelling.rst index 6f3adce6e..a63db3f06 100644 --- a/docs/source/user_guide/fold_modelling.rst +++ b/docs/source/user_guide/fold_modelling.rst @@ -1,4 +1,110 @@ Fold modelling -=============== +============== -A specific interpolator can be used for modelling folds geological features. The LoopStructural implements \ No newline at end of file +Fold interpolators constrain a scalar-field gradient using a fold frame, a fold +axis and a fitted limb-rotation profile. ``DFI`` uses linear tetrahedra; +``P2Fold`` uses quadratic tetrahedra. Select either when creating a folded +foliation or a folded fold frame. The latter uses the fold interpolator for its +first coordinate; its other two coordinates still use FDI. + +Using P2Fold +------------ + +``P2Fold`` with its default ``gauss4`` sampling is useful when sparse orientation +data leave the fold constraints to control the field between observations. Its +four interior samples constrain the gradient throughout a quadratic tetrahedron. +In the Bermagui validation it reduced median fold-direction error with fewer +unknowns than DFI for both S2 and S0. This comparison measures agreement with a +shared fold event, rather than independent geological truth, and does not imply +that P2Fold is always the better interpolator. + +Given an existing model and fold frame:: + + folded = model.create_and_add_folded_foliation( + "s0", + fold_frame=frame, + interpolatortype="P2Fold", + nelements=1000, + buffer=0.0, + limb_wl=None, + fold_weights={"sample_rule": "gauss4"}, + solver="cg", + solver_kwargs={"maxiter": 20000, "rtol": 1e-8}, + ) + values = folded.evaluate_value(points) # triggers the lazy build + print(folded.builder.interpolator.dof) + +The same ``interpolatortype`` and ``fold_weights`` options can be passed to +``create_and_add_folded_fold_frame``. The bundled synthetic example is +:ref:`sphx_glr__auto_examples_2_fold_plot_3_p2_fold.py`. + +Sampling and mesh size +---------------------- + +``fold_weights["sample_rule"]`` accepts ``"gauss4"`` (default) or ``"centroid"``. +The library weights each sample residual by ``w * sqrt(q_i)``, with quadrature +weights summing to one in each tetrahedron. Adding more samples therefore does +not multiply the total normalised operator weight by the sample count. + +A single centroid sample does not adequately control a varying P2 gradient when +regularisation is disabled. The unregularised Bermagui centroid runs had about +60 degrees median error and unstable field scales, even when CG converged. +Use gauss4 for this configuration; increasing mesh size did not rescue centroid. +This result does not establish the behaviour of every regularised centroid model. + +Start with a small ``nelements`` request and inspect the realised unknown and +tetrahedron counts:: + + interpolator = folded.builder.interpolator + print(interpolator.support.n_elements, interpolator.dof) + +The current support factory treats ``nelements`` as a structured-grid cell +budget, then splits each cell into five tetrahedra. P2 adds nodes on the edges. +It is not a literal tetrahedron count, and equal requests are not equal cost +between DFI and P2Fold. For the Bermagui box with ``buffer=0.0``: + +.. list-table:: Realised P2 meshes + :header-rows: 1 + :widths: 25 35 40 + + * - Requested nelements + - Tetrahedra + - Unknowns + * - 1,000 + - 5,850 + - 9,711 + * - 10,000 + - 56,420 + - 85,673 + * - 100,000 + - 525,690 + - 765,347 + +The largest gauss4 S0 case took about 93 seconds in ``solve_system`` and used +about 6.6 GiB peak resident memory on the validation machine. These are single +measurements, not runtime or memory guarantees. The +`Bermagui handoff `_ +records the corrected comparison summary and limitations. The scripts in +``scripts/validation`` reproduce the detailed results and figures locally. + +Interpreting the result +----------------------- + +At the smaller budget, S2 gauss4 had 4.69 degrees median error with 9,711 unknowns, +compared with DFI's 6.15 degrees at 12,992 unknowns. S0 had 11.14 degrees at 85,673 +unknowns, compared with 20.63 degrees at 112,560 unknowns. S0 still had 19.4 percent +of inner probes above 30 degrees; inspect the error distribution and fitted +profiles rather than relying only on the median. + +The scalar-field scale is arbitrary. The norm constraint need not be met, so the +field should not be interpreted as a distance function. A wavelength expressed +in field units is not portable between independently fitted models. Inspect each +rotation profile with its observations. The Bermagui S0 validation uses +``limb_wl=None`` (automatic fitting); starts between 2 and 5 were poor local +minima in the earlier experiments. Other datasets require their own checks. + +Remaining issues include buffer application at both model-box and support +construction, lower-face NaNs in the linear tetrahedral support (F2), amplitude +damping (F3), and rotation-profile fit instability. Use the same explicit buffer +when comparing engines. ``buffer=0.0`` keeps this validation inside the model box; +boundary queries can return NaN. These issues are separate follow-up work. diff --git a/examples/2_fold/plot_3_p2_fold.py b/examples/2_fold/plot_3_p2_fold.py new file mode 100644 index 000000000..7d32a28a3 --- /dev/null +++ b/examples/2_fold/plot_3_p2_fold.py @@ -0,0 +1,58 @@ +""" +2c. Quadratic fold interpolation +============================== + +``P2Fold`` uses quadratic tetrahedra and samples the fold constraints inside +each tetrahedron. Select it when constructing a folded foliation or the first +coordinate of a folded fold frame. ``gauss4`` is the default sample rule. + +Start with a small mesh and inspect the realised unknown count. The current +factory interprets ``nelements`` as a grid-cell budget, then splits each cell +into five tetrahedra; P2 also adds nodes on the tetrahedron edges. +""" + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from LoopStructural import GeologicalModel +from LoopStructural.datasets import load_noddy_single_fold + +# This bundled synthetic dataset makes the example independent of Bermagui data. +data, bounds = load_noddy_single_fold() +model = GeologicalModel(bounds[0], bounds[1]) +model.set_model_data(pd.concat([data[:100], data[data["feature_name"] == "s1"]])) +frame = model.create_and_add_fold_frame("s1", nelements=1000, buffer=0.2) +folded = model.create_and_add_folded_foliation( + "s0", + fold_frame=frame, + interpolatortype="P2Fold", + nelements=1000, + buffer=0.0, + fold_axis=[-6.51626577e-06, -5.00013645e-01, -8.66017526e-01], + limb_wl=1, + fold_weights={"sample_rule": "gauss4"}, + solver="cg", + solver_kwargs={"maxiter": 20000, "rtol": 1e-8}, +) + +# Builds are lazy: evaluation triggers the fit and solve. +x = np.linspace( + bounds[0, 0] + 0.1 * np.ptp(bounds[:, 0]), bounds[1, 0] - 0.1 * np.ptp(bounds[:, 0]), 60 +) +z = np.linspace( + bounds[0, 2] + 0.1 * np.ptp(bounds[:, 2]), bounds[1, 2] - 0.1 * np.ptp(bounds[:, 2]), 60 +) +xx, zz = np.meshgrid(x, z) +points = np.c_[xx.ravel(), np.full(xx.size, np.mean(bounds[:, 1])), zz.ravel()] +values = folded.evaluate_value(points) +assert np.isfinite(values).all() +interpolator = folded.builder.interpolator +print(f"P2Fold: {interpolator.support.n_elements:,} tetrahedra, {interpolator.dof:,} unknowns") + +fig, ax = plt.subplots(figsize=(7, 5), layout="constrained") +contours = ax.contourf(xx, zz, values.reshape(xx.shape), levels=20, cmap="viridis") +ax.contour(xx, zz, values.reshape(xx.shape), levels=10, colors="white", linewidths=0.5) +ax.set(xlabel="X", ylabel="Z", title="P2Fold with gauss4: section through the synthetic fold") +fig.colorbar(contours, ax=ax, label="Folded scalar field (arbitrary units)") +plt.show() diff --git a/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py index 77a19472a..caedc7ee7 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_p2_fold_interpolator.py @@ -54,7 +54,7 @@ def __init__(self, support, fold: FoldEvent | None = None): self.fold = fold logger.info( f"P2Fold support: {support.elements.shape[0]} tets, {support.nodes.shape[0]} unknowns. " - "nelements counts tets for P2, so cost grows quickly with it" + "factory nelements requests grid cells; inspect actual tet and unknown counts" ) def copy(self): diff --git a/pyproject.toml b/pyproject.toml index 712f690c2..d7e311799 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -97,6 +97,7 @@ LoopStructural = { workspace = true } [tool.pytest.ini_options] addopts = "--import-mode=importlib" testpaths = ["tests"] +pythonpath = ["."] [tool.isort] profile = 'black' diff --git a/scripts/validation/p2_fold_bermagui.py b/scripts/validation/p2_fold_bermagui.py new file mode 100644 index 000000000..40f44c947 --- /dev/null +++ b/scripts/validation/p2_fold_bermagui.py @@ -0,0 +1,335 @@ +"""Validate library P2Fold against Bermagui data held outside this repository. + +Run from the repository root. Each candidate gets a fresh subprocess and the +same saved upstream DFI frames, so peak RSS and profile fits are comparable. +No prototype class, interpolator replacement, or private data is copied here. +""" + +from __future__ import annotations + +import argparse +import csv +import ctypes +import hashlib +import json +import os +from pathlib import Path +import random +import subprocess +import sys +import time + +import numpy as np + +ORIGIN = np.array([238275.0, 5965000.0, -50.0]) +MAXIMUM = np.array([238500.0, 5965250.0, 50.0]) +AXIS = [-0.31882007, -0.46674375, -0.06258398] +SEED = 12345 +SOLVER = {"maxiter": 20000, "rtol": 1e-8} +FILES = {s: f"TrimData{s}.csv" for s in ("S4", "S3", "S2")} +FILES["S0"] = "TrimDataS0_LA2.csv" + + +def reseed(): + from loop_common.math import rng + + random.seed(SEED) + np.random.seed(SEED) + rng.bit_generator.state = np.random.default_rng(SEED).bit_generator.state + + +def peak_rss_mib(): + if os.name != "nt": + import resource + + value = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + return value / (1024**2 if sys.platform == "darwin" else 1024) + + class Counters(ctypes.Structure): + _fields_ = [("cb", ctypes.c_ulong), ("PageFaultCount", ctypes.c_ulong)] + [ + (name, ctypes.c_size_t) + for name in ( + "PeakWorkingSetSize", + "WorkingSetSize", + "QuotaPeakPagedPoolUsage", + "QuotaPagedPoolUsage", + "QuotaPeakNonPagedPoolUsage", + "QuotaNonPagedPoolUsage", + "PagefileUsage", + "PeakPagefileUsage", + ) + ] + + counts = Counters() + counts.cb = ctypes.sizeof(counts) + handle = ctypes.windll.kernel32.GetCurrentProcess() + ctypes.windll.psapi.GetProcessMemoryInfo.argtypes = [ + ctypes.c_void_p, + ctypes.POINTER(Counters), + ctypes.c_ulong, + ] + if not ctypes.windll.psapi.GetProcessMemoryInfo(handle, ctypes.byref(counts), counts.cb): + raise ctypes.WinError() + return counts.PeakWorkingSetSize / 1024**2 + + +def data(root, name): + import pandas as pd + + return pd.read_csv(root / "Data" / "Processed" / FILES[name]) + + +def save_json(path, payload): + path.write_text(json.dumps(payload, indent=2, sort_keys=True, allow_nan=False) + "\n") + + +def prepare(args): + import dill + from LoopStructural import BoundingBox, GeologicalModel + + reseed() + model = GeologicalModel(BoundingBox(origin=ORIGIN, maximum=MAXIMUM)) + model.set_model_data(data(args.project, "S4")) + s4 = model.create_and_add_fold_frame( + "S4", + nelements=10000, + buffer=0.5, + solver="cg", + damp=True, + solver_kwargs=SOLVER, + ) + s3 = model.create_and_add_folded_fold_frame( + "S3", + limb_wl=100, + data=data(args.project, "S3"), + fold_frame=s4, + av_fold_axis=True, + nelements=10000, + buffer=0.3, + lag=0.05, + solver="cg", + solver_kwargs=SOLVER, + interpolatortype="DFI", + ) + s3[0].evaluate_value((ORIGIN + 1)[None, :]) + reseed() + s2 = model.create_and_add_folded_fold_frame( + "S2", + limb_wl=80, + data=data(args.project, "S2"), + fold_frame=s3, + av_fold_axis=True, + nelements=10000, + buffer=args.buffer, + lag=0.05, + solver="cg", + solver_kwargs=SOLVER, + interpolatortype="DFI", + ) + # Build all coordinates used by a downstream fold before serialising. + for frame in (s3, s2): + for feature in frame: + feature.evaluate_value((ORIGIN + 1)[None, :]) + with (args.out / "drivers.pkl").open("wb") as stream: + dill.dump(model, stream) + save_json( + args.out / "inputs.json", + { + "seed": SEED, + "origin": ORIGIN.tolist(), + "maximum": MAXIMUM.tolist(), + "target_buffer": args.buffer, + "solver": "cg", + "solver_kwargs": SOLVER, + "input_sha256": { + name: hashlib.sha256( + (args.project / "Data" / "Processed" / file).read_bytes() + ).hexdigest() + for name, file in FILES.items() + }, + "commit": subprocess.check_output(["git", "rev-parse", "HEAD"], text=True).strip(), + }, + ) + + +def unit(vectors): + with np.errstate(divide="ignore", invalid="ignore"): + return vectors / np.linalg.norm(vectors, axis=1, keepdims=True) + + +def worker(args): + import dill + from loop_interpolation import DiscreteFoldInterpolator, P2FoldInterpolator + + with (args.out / "drivers.pkl").open("rb") as stream: + model = dill.load(stream) + reseed() + driver = model["S3"] if args.target == "S2" else model["S2"] + options = ( + {"limb_wl": 80, "av_fold_axis": True} + if args.target == "S2" + else { + "limb_wl": None, + "fold_axis": AXIS, + } + ) + folded = model.create_and_add_folded_fold_frame( + args.target + "_validation", + data=data(args.project, args.target), + fold_frame=driver, + interpolatortype=args.engine, + nelements=args.nelements, + buffer=args.buffer, + lag=0.05, + solver="cg", + solver_kwargs=dict(SOLVER), + fold_weights={"sample_rule": args.rule} if args.engine == "P2Fold" else {}, + **options, + ) + feature = folded[0] + interp = feature.builder.interpolator + assert type(interp) is ( + P2FoldInterpolator if args.engine == "P2Fold" else DiscreteFoldInterpolator + ) + solve_seconds = [] + original = interp.solve_system + + def timed_solve(*pos, **kw): + start = time.perf_counter() + ok = original(*pos, **kw) + solve_seconds.append(time.perf_counter() - start) + if not ok: + raise RuntimeError("Candidate CG did not converge") + return ok + + interp.solve_system = timed_solve + start = time.perf_counter() + feature.evaluate_value((ORIGIN + 1)[None, :]) + build_seconds = time.perf_counter() - start + assert solve_seconds, "Lazy build did not solve the candidate" + # Same 60 x 60 x 30 grid and inner 70% box as the prototype. + axes = [np.linspace(a, b, n) for a, b, n in zip(ORIGIN, MAXIMUM, (60, 60, 30))] + points = np.array(np.meshgrid(*axes, indexing="ij")).reshape(3, -1).T + fraction = (points - ORIGIN) / (MAXIMUM - ORIGIN) + points = points[np.all((fraction > 0.15) & (fraction < 0.85), axis=1)] + gradient = feature.evaluate_gradient(points) + fdir, axis, _ = folded.fold.get_deformed_orientation(points) + expected = unit(np.cross(unit(fdir), unit(axis))) + angles = np.degrees( + np.arccos(np.clip(np.abs(np.sum(unit(gradient) * expected, axis=1)), 0.0, 1.0)) + ) + finite = np.isfinite(angles) + if not finite.all(): + raise RuntimeError(f"Nonfinite inner probes: {np.sum(~finite)}/{len(points)}") + profile = folded.fold.fold_limb_rotation + xs = np.asarray(profile.fold_frame_coordinate) + ys = np.asarray(profile.rotation_angle) + used = np.isfinite(xs) & np.isfinite(ys) & (np.abs(ys) < 89) + residual = (profile(xs[used]) - ys[used] + 90) % 180 - 90 + row = { + "target": args.target, + "engine": args.engine, + "rule": args.rule, + "requested_nelements": args.nelements, + "buffer": args.buffer, + "unknowns": int(interp.dof), + "tets": int(interp.support.n_elements), + "median_angle_deg": float(np.median(angles)), + "fraction_above_30_deg": float(np.mean(angles > 30)), + "solve_seconds": float(sum(solve_seconds)), + "build_seconds": build_seconds, + "peak_rss_mib": peak_rss_mib(), + "finite_probe_fraction": float(finite.mean()), + "mean_gradient_norm": float(np.linalg.norm(gradient, axis=1).mean()), + "profile_rms_mod180_deg": float(np.sqrt(np.mean(residual**2))), + "expected_normal_sha256": hashlib.sha256(np.round(expected, 8).tobytes()).hexdigest(), + "profile_params": {k: float(v) for k, v in profile.params.items()}, + } + label = f"{args.target}_{args.engine}_{args.rule}_{args.nelements}" + np.savez_compressed( + args.out / (label + ".npz"), + points=points, + angle_deg=angles, + gradient=gradient, + expected=expected, + profile_x=xs, + profile_y=ys, + curve_x=np.linspace(xs.min(), xs.max(), 2000), + curve_y=profile(np.linspace(xs.min(), xs.max(), 2000)), + ) + save_json(args.out / (label + ".json"), row) + print("RESULT", json.dumps(row), flush=True) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--project", type=Path, required=True) + parser.add_argument("--out", type=Path, default=Path(".pytest_cache/p2_fold/runtime")) + parser.add_argument("--buffer", type=float, default=0.0) + parser.add_argument("--mode", choices=("prepare", "worker", "all"), default="all") + parser.add_argument("--target", choices=("S2", "S0"), default="S2") + parser.add_argument("--engine", choices=("DFI", "P2Fold"), default="DFI") + parser.add_argument("--rule", choices=("centroid", "gauss4"), default="gauss4") + parser.add_argument("--nelements", type=int, default=10000) + args = parser.parse_args() + args.project = args.project.resolve() + args.out = args.out.resolve() + args.out.mkdir(parents=True, exist_ok=True) + if args.mode == "prepare": + prepare(args) + elif args.mode == "worker": + worker(args) + else: + base = [ + sys.executable, + str(Path(__file__).resolve()), + "--project", + str(args.project), + "--out", + str(args.out), + "--buffer", + str(args.buffer), + ] + jobs = [("prepare", [])] + for target, nelements in (("S2", 10000), ("S0", 100000)): + configs = [("DFI", "gauss4", nelements)] + configs += [ + ("P2Fold", rule, nel) + for nel in (nelements, nelements // 10) + for rule in ("centroid", "gauss4") + ] + for engine, rule, nel in configs: + jobs.append( + ( + f"{target}_{engine}_{rule}_{nel}", + [ + "--target", + target, + "--engine", + engine, + "--rule", + rule, + "--nelements", + str(nel), + ], + ) + ) + for label, extra in jobs: + print("START", label, flush=True) + with (args.out / (label + ".log")).open("w") as log: + subprocess.run( + base + ["--mode", "prepare" if label == "prepare" else "worker"] + extra, + stdout=log, + stderr=subprocess.STDOUT, + check=True, + ) + print("DONE", label, flush=True) + rows = [json.loads(p.read_text()) for p in sorted(args.out.glob("S*.json"))] + with (args.out / "results.csv").open("w", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + +if __name__ == "__main__": + main() diff --git a/scripts/validation/summarise_p2_fold.py b/scripts/validation/summarise_p2_fold.py new file mode 100644 index 000000000..12c33c90b --- /dev/null +++ b/scripts/validation/summarise_p2_fold.py @@ -0,0 +1,165 @@ +"""Summarise a completed Bermagui comparison and render review figures.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.ticker import FuncFormatter, LogLocator, NullFormatter +import numpy as np + + +def decision(rows): + """Require comparable fold events and a lower median with no more unknowns.""" + result = {} + for target in ("S2", "S0"): + candidates = [r for r in rows if r["target"] == target] + if not candidates: + raise ValueError(f"Missing {target} results") + if len({r["expected_normal_sha256"] for r in candidates}) != 1: + raise ValueError(f"{target} candidates used different fold events") + if len({r["buffer"] for r in candidates}) != 1: + raise ValueError(f"{target} candidates used different buffers") + (reference,) = [r for r in candidates if r["engine"] == "DFI"] + eligible = [ + r + for r in candidates + if r["engine"] == "P2Fold" + and r["rule"] == "gauss4" + and r["unknowns"] <= reference["unknowns"] + ] + if not eligible: + raise ValueError(f"No {target} P2Fold candidate within the DFI nodal budget") + best = min(eligible, key=lambda r: r["median_angle_deg"]) + result[target] = { + "passed": best["median_angle_deg"] < reference["median_angle_deg"], + "dfi_unknowns": reference["unknowns"], + "p2_unknowns": best["unknowns"], + "dfi_median_deg": reference["median_angle_deg"], + "p2_median_deg": best["median_angle_deg"], + } + result["passed"] = all(result[t]["passed"] for t in ("S2", "S0")) + return result + + +def render(runtime, output, rows): + output.mkdir(parents=True, exist_ok=True) + fig, axes = plt.subplots(2, 3, figsize=(13, 7), layout="constrained") + for i, target in enumerate(("S2", "S0")): + for engine, rule, color, marker in ( + ("DFI", "gauss4", "black", "s"), + ("P2Fold", "centroid", "tab:orange", "o"), + ("P2Fold", "gauss4", "tab:blue", "o"), + ): + group = sorted( + ( + r + for r in rows + if r["target"] == target and r["engine"] == engine and r["rule"] == rule + ), + key=lambda r: r["unknowns"], + ) + label = "DFI" if engine == "DFI" else "P2Fold " + rule + for j, (key, title) in enumerate( + ( + ("median_angle_deg", "Median axial direction error (deg)"), + ("fraction_above_30_deg", "Share above 30 deg (%)"), + ("solve_seconds", "Solve time (s)"), + ) + ): + factor = 100 if j == 1 else 1 + axes[i, j].plot( + [r["unknowns"] for r in group], + [r[key] * factor for r in group], + color=color, + marker=marker, + label=label, + ) + axes[i, j].set(xscale="log", xlabel="Unknowns", ylabel=title, title=target) + axes[i, j].xaxis.set_major_locator(LogLocator(base=10, subs=(1, 3))) + axes[i, j].xaxis.set_major_formatter(FuncFormatter(lambda x, _: f"{x / 1000:g}k")) + axes[i, j].xaxis.set_minor_formatter(NullFormatter()) + axes[i, j].grid(alpha=0.2) + axes[i, 0].legend(fontsize=8) + fig.suptitle("Bermagui: corrected library fold weights, fixed upstream frames, buffer 0") + fig.savefig(output / "comparison.png", dpi=150) + plt.close(fig) + + fig, axes = plt.subplots(2, 3, figsize=(13, 7), layout="constrained") + for i, target in enumerate(("S2", "S0")): + baseline = next(r for r in rows if r["target"] == target and r["engine"] == "DFI") + for j, (engine, rule) in enumerate( + (("DFI", "gauss4"), ("P2Fold", "centroid"), ("P2Fold", "gauss4")) + ): + row = next( + r + for r in rows + if r["target"] == target + and r["engine"] == engine + and r["rule"] == rule + and r["requested_nelements"] == baseline["requested_nelements"] + ) + label = f"{target}_{engine}_{rule}_{row['requested_nelements']}" + with np.load(runtime / (label + ".npz")) as z: + p, a = z["points"], z["angle_deg"] + nearest_y = np.unique(p[:, 1])[np.argmin(np.abs(np.unique(p[:, 1]) - 5965125))] + m = p[:, 1] == nearest_y + ax = axes[i, j] + artist = ax.pcolormesh( + np.unique(p[m, 0]) - 238275, + np.unique(p[m, 2]), + a[m].reshape(len(np.unique(p[m, 0])), -1).T, + shading="nearest", + vmin=0, + vmax=60, + cmap="magma_r", + ) + ax.set( + title=f"{target} {engine} {rule if engine != 'DFI' else ''}\n" + f"{row['unknowns']:,} unknowns; median {row['median_angle_deg']:.2f} deg", + xlabel="Easting minus 238275 (m)", + ylabel="Elevation (m)", + ) + fig.colorbar(artist, ax=axes, label="Axial direction error (deg)", shrink=0.8) + fig.savefig(output / "sections.png", dpi=150) + plt.close(fig) + + fig, axes = plt.subplots(1, 2, figsize=(12, 4), layout="constrained") + for ax, target in zip(axes, ("S2", "S0")): + row = next(r for r in rows if r["target"] == target and r["engine"] == "DFI") + label = f"{target}_DFI_gauss4_{row['requested_nelements']}" + with np.load(runtime / (label + ".npz")) as z: + ax.plot(z["curve_x"], z["curve_y"], label="Fitted limb profile") + ax.scatter(z["profile_x"], z["profile_y"], color="tab:red", label="Observations") + ax.set( + title=f"{target}: shared profile; RMS {row['profile_rms_mod180_deg']:.2f} deg", + xlabel="Driver field coordinate (arbitrary scale)", + ylabel="Rotation angle (deg)", + ylim=(-100, 100), + ) + ax.legend(fontsize=8) + fig.savefig(output / "profiles.png", dpi=150) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("runtime", type=Path) + parser.add_argument("--output", type=Path, default=Path(".pytest_cache/p2_fold")) + args = parser.parse_args() + rows = [json.loads(p.read_text()) for p in sorted(args.runtime.glob("S*.json"))] + outcome = decision(rows) + args.output.mkdir(parents=True, exist_ok=True) + (args.output / "results.json").write_text(json.dumps(rows, indent=2, sort_keys=True) + "\n") + (args.output / "decision.json").write_text(json.dumps(outcome, indent=2, sort_keys=True) + "\n") + render(args.runtime, args.output, rows) + print(json.dumps(outcome, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tests/unit/validation/test_p2_fold_validation.py b/tests/unit/validation/test_p2_fold_validation.py new file mode 100644 index 000000000..c510d84c4 --- /dev/null +++ b/tests/unit/validation/test_p2_fold_validation.py @@ -0,0 +1,49 @@ +"""Check that the validation gate cannot endorse an unfair comparison.""" + +import pytest + +from scripts.validation.summarise_p2_fold import decision + + +def rows(): + return [ + { + "target": target, + "engine": engine, + "rule": "gauss4", + "buffer": 0.0, + "unknowns": unknowns, + "median_angle_deg": angle, + "expected_normal_sha256": target, + } + for target in ("S2", "S0") + for engine, unknowns, angle in (("DFI", 12000, 8.0), ("P2Fold", 9000, 5.0)) + ] + + +def test_gate_requires_improvement_for_both_targets(): + results = rows() + assert decision(results)["passed"] + results[-1]["median_angle_deg"] = 9.0 + assert not decision(results)["passed"] + + +def test_gate_rejects_larger_p2_mesh(): + results = rows() + results[1]["unknowns"] = 85000 + with pytest.raises(ValueError, match="nodal budget"): + decision(results) + + +def test_gate_rejects_different_fold_geometry(): + results = rows() + results[1]["expected_normal_sha256"] = "different profile" + with pytest.raises(ValueError, match="different fold events"): + decision(results) + + +def test_gate_rejects_different_buffers(): + results = rows() + results[1]["buffer"] = 0.2 + with pytest.raises(ValueError, match="different buffers"): + decision(results)