LoopStructural 1.7 - #290
Open
lachlangrose wants to merge 68 commits into
Open
Conversation
add_gradient_constraints indexed self.support directly instead of self.support.elements, and evaluate_d2 used `points == np.nan` (always False) instead of np.isnan, so NaN filtering never ran.
The chunking loop sliced points[:npts+npts_step] instead of points[npts:npts+npts_step], so every iteration recomputed from index 0 (quadratic blowup on a hot path used by every value/gradient evaluation), and its `inside` variable name shadowed the pre-allocated output array. Verified against a 25k-point set spanning multiple chunk boundaries against a brute-force barycentric reference.
The nan-check in add_constraints_to_least_squares only warned (the rejecting return was commented out) and never checked for inf at all, letting bad values reach the sparse solver. solve_system's solver_kwargs mutable default dict was also being mutated in place (atol/btol set, x0/linsys_solver/admm_weight popped), leaking state across calls that omit the argument.
faults/regions used mutable list defaults and regions was assigned directly without copying, so calling add_region() on one feature leaked the region into every other feature created without an explicit regions argument.
add_abutting/faulted_relationship inserted a dead string-keyed entry before the real tuple-keyed one, which broke unpacking in get_fault_relationships/get_matrix for any 2+ char fault name. update_from_dict also iterated adjacency.values() instead of .items(), dropping the ABUTTING/FAULTED distinction and iterating the wrong objects, so to_dict/from_dict was not actually round-trippable. Verified with a standalone round-trip check (add both relationship types, serialise, deserialise, compare).
add_value_constraints required points.shape[0] > 1, silently discarding a single value constraint, inconsistent with the finite-difference interpolator's equivalent check (> 0).
v_t was normalised by its own norm with no zero-guard, so a degenerate (zero-gradient) element produced nan/inf that got fed into the least squares system. Now zero-norm elements are skipped and logged instead. Verified with an all-zero-coefficient case (previously nan) and a normal random-coefficient case (unaffected).
…setter Catching bare except and re-raising BaseException directly swallowed KeyboardInterrupt/SystemExit and is an anti-pattern for callers trying to catch this with `except Exception`. Verified the setter now raises ValueError for invalid data input.
…id export add_structured_grid_to_omf just printed a message and returned None, which StructuredGrid.save() treated as a successful save. OMF genuinely can't represent structured grids, so raise NotImplementedError instead of silently doing nothing.
ValuePoints.save/VectorPoints.save imported _write_pointset from gocad.py, but that function is commented out (and incomplete) there, so saving to .vs raised a confusing ImportError. Raise NotImplementedError with a clear message instead until a real writer exists.
__all__ only listed GeologicalModel while the module also exposes StratigraphicColumn, FaultTopology, LoopInterpolator, InterpolatorBuilder, BoundingBox, and several logging helpers, so `from LoopStructural import *` and API doc tooling didn't reflect what's actually public.
Both files had near-zero coverage and were exactly where the earlier
bug-fixing pass found real, previously-undetected defects.
FaultTopology tests cover the full relationship API (add/remove/update/
change, both fault-fault and fault-stratigraphy relationships) plus the
to_dict/from_dict round trip. Writing the round-trip test surfaced a
second bug in update_from_dict: it treated
stratigraphy_fault_relationships as {unit: [fault, ...]} when to_dict
actually produces {(unit, fault): bool}, raising TypeError on any
populated round trip. Fixed alongside the tests.
P2Interpolator tests use a minimal fake support (no real quadratic mesh
needed) to exercise add_gradient_constraints, evaluate_d2's nan
masking, and the single-point value-constraint case in isolation.
Verified all three regression tests actually fail against the pre-fix
code before confirming they pass against the fix.
…it arrays P2UnstructuredTetMesh (and therefore the whole "P2"/piecewise-quadratic interpolator type) could previously only be constructed from explicit nodes/elements/neighbours arrays. But InterpolatorFactory.create_interpolator - the standard way GeologicalModel and InterpolatorBuilder construct any interpolator - always calls support classes with origin/step_vector/nsteps, so building a P2 interpolator via the documented API always raised TypeError. This made the entire P2 interpolator type unreachable through normal use, which lines up with its near-zero test coverage. __init__ now accepts optional origin/step_vector/nsteps and builds the mesh by tessellating the bounding box into linear tets (via TetMesh) and adding a deduplicated midpoint node per edge to elevate each tet to a 10-node quadratic element. Explicit nodes/elements/neighbours construction still works unchanged. The local edge -> node ordering (corners 0-3, then edges (2,3)->4, (0,3)->5, (0,1)->6, (1,2)->7, (1,3)->8, (0,2)->9) was reverse-engineered from evaluate_shape's basis functions and verified by fitting a genuine quadratic scalar field at every mesh node and confirming the interpolant reproduces it to machine precision elsewhere in the domain (with the default smoothing regularisation disabled, since that intentionally perturbs exact interpolation).
…/P2 constraints add_value_constraints, add_norm_constraints and add_gradient_orthogonal_constraints in both interpolators hardcoded points[:, :3], points[:, 3], points[:, 3:6] and a tile count of 3, unconditionally assuming 3D position/vector columns. For a genuinely 2D interpolator (dimensions=2), points[:, :3] silently pulled in the *value* column as if it were a Z coordinate, corrupting every constraint. Since self.dimensions == 3 for all existing 3D usage, this is a no-op there and only changes behaviour for 2D interpolators.
…structured meshes get_element_for_location computed barycentric weights for vertices[:,1]/vertices[:,2]/vertices[:,0] (in that order, standard Ericson technique) but stored them in columns [0,1,2] - i.e. every returned barycentric/shape-function array was rotated one column off from the vertex order used everywhere else (elements[tri,:], evaluate_value, add_value_constraints). This silently corrupted every 2D value/gradient evaluation rather than crashing. Also fixed the same chunking bug found earlier in the 3D tetra mesh: the loop sliced points[:npts+npts_step] instead of points[npts:npts+npts_step], and its local `inside` variable shadowed the pre-allocated output array. In 2D this wasn't just a performance issue - points within the padded AABB box but not actually inside any triangle inherited a stale `inside=True` from the coarse pre-filter, so they were assigned a bogus tri=0 with all-zero shape weights instead of correctly reporting as outside (nan).
…rays Same root cause as the earlier 3D P2 fix: P1Unstructured2d only accepted explicit elements/vertices/neighbours arrays, so building a 2D P1 interpolator via InterpolatorFactory/InterpolatorBuilder always raised TypeError. __init__ now accepts optional origin/step_vector/ nsteps and tessellates the grid into triangles (splitting each cell along the bottom-left/top-right diagonal), with a vectorised neighbour computation that groups shared edges and cross-references the two triangles either side of each interior edge.
…gradient evaluation
Same construction fix as P1Unstructured2d and the 3D P2 tet mesh:
__init__ now accepts origin/step_vector/nsteps and elevates the P1
triangulation to 6-node quadratic triangles by adding a deduplicated
midpoint node per edge (local ordering: corners 0-2, then edges
(1,2)->3, (0,2)->4, (0,1)->5, matching evaluate_shape's basis
functions).
Also added evaluate_value/evaluate_gradient overrides. The inherited
BaseUnstructured2d implementations only use the 3 linear barycentric
weights against a 6-column elements array, which raised a shape
mismatch (200,3) vs (200,6) - P2Unstructured2d never had quadratic
evaluation of its own, only the shape-function/derivative plumbing
used when building constraints.
Verified end-to-end via InterpolatorFactory.create_interpolator('P2',
2d_bounding_box, ...): fitting a genuine quadratic field at every mesh
node and reproducing it to machine precision elsewhere in the domain.
Covers: bbox-construction error handling and topology validity for both P1Unstructured2d and P2Unstructured2d, factory construction for both interpolator types, and end-to-end field reproduction (linear for P1, quadratic for P2) via InterpolatorFactory with a 2D bounding box. Verified all 8 tests fail against the pre-fix code (reverted the 5 touched files, confirmed every test fails with the original TypeError/ wrong-value symptoms, then restored the fixes) before committing.
cp was built with shape (n_edges, self.ncps, 2) - using the element node count (6) instead of the actual number of quadrature points (2) - and the returned weight array inherited that same wrong shape. This only surfaced when minimise_edge_jumps (used by the default constant-gradient regularisation) actually ran on a 2D P2 interpolator, which nothing exercised until the 2D bbox-construction fix made it reachable. Verified the fix against the pre-fix code (confirmed both new regression tests fail with the original ValueError) before committing.
…D interpolation Fits a scattered sample of Franke's function (a standard scattered-data interpolation benchmark) with scipy.interpolate.RBFInterpolator and with LoopStructural's P1/P2 interpolators built on a 2D bounding box, then compares the reconstructed surfaces and RMSE against the known function. P2 is run with add_value_constraints + minimise_edge_jumps directly rather than setup_interpolator(), since the latter's default regularisation also calls minimise_grad_steepness(), which depends on 2D second-derivative shape-function code that was never finished (P2Unstructured2d.evaluate_shape_d2 references a self.hN attribute that's never set). The discussion section explains this honestly rather than tuning P2 to look artificially good - its RMSE is noticeably worse than P1's here because it's missing that curvature regularisation term.
evaluate_shape_d2 referenced self.hN, an attribute that's never set anywhere in the class - it raised AttributeError the moment minimise_grad_steepness (P2's curvature-minimising regularisation term) actually ran, which nothing did until the 2D bbox-construction fix made 2D P2 reachable at all. Rewritten to follow the same approach as P2UnstructuredTetMesh's 3D evaluate_shape_d2: use the reference-space hessian of the quadratic shape functions (self.hessian, already built in __init__ and confirmed by hand to match the shape functions in evaluate_shape) and the inverse Jacobian chain rule to get physical second derivatives. Also fixed evaluate_d2 (the only in-file caller), which unpacked a 3-tuple that evaluate_shape_d2 never returned and called evaluate_shape() expecting 2 return values instead of 3 - both symptoms of code that had never actually run. Verified evaluate_shape_d2 exactly reproduces the analytic second derivatives of a genuine quadratic field, and that with this fix P2's full setup_interpolator() regularisation (edge-jump + curvature) now gives a comparable RMSE to P1 and scipy's RBF on scattered data, instead of the ~15x worse fit it had when only edge-jump smoothing was available (see the interpolation comparison example). Confirmed both new regression tests fail against the pre-fix code before committing.
…_to_geoh5 function
Covers export/{exporters,gocad,omf_wrapper}.py and
interpolators/{_operator,_api,_interpolator_factory}.py, previously untested.
Documents 3 pre-existing bugs surfaced during testing (VTK surface export
AttributeError, GOCAD volume export LoopValueError, OMF pointset API
mismatch) via characterization tests / strict xfail rather than fixing them.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
(cherry picked from commit 2c1f2b962480cf07f1accad56d9497744d9a3eea)
Covers utils/{helper,observer,regions,_transformation,_surface}.py and
modelling/features/{_region,_unconformity_feature,fold/_svariogram}.py,
previously untested. Documents 8 pre-existing bugs surfaced during testing
(RegionEverywhere/RegionFunction unconstructable, EuclideanTransformation
.inverse_transform slicing bug, svariogram lag-cap and dtype issues, observer
edge cases) via clearly-commented characterization tests rather than fixing
them.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
(cherry picked from commit 1f6fa4029da588961406d674c33ee0fa2009451c)
…into API Replaces unfilled `_description_` numpydoc stubs across interpolators/, modelling/, datatypes/, and utils/ with real descriptions, fixes DiscreteInterpolator.evaluate_value/evaluate_gradient docstrings which documented a since-renamed `evaluation_points` parameter as `locations`, and adds LoopStructural.utils to the Sphinx API reference autosummary (previously missing, so its docstrings never reached the built docs). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> (cherry picked from commit 98f237e8f9c23d92212c48f7cf516863d53b9146)
Narrows broad except-Exception/bare-except blocks to the specific exception types each call site can actually raise, so unrelated bugs stop being silently swallowed as "nan slicing" or generic warnings. Guards GeologicalModel.load's pickle.load with a file-existence check and a LoopValueError on failure instead of raising raw pickle/OS errors. Replaces two validation-purpose `assert isinstance(...)` calls (stripped under python -O) in GeologicalModel with explicit `raise TypeError`. utils/observer.py's broad except is left as-is: it deliberately catches arbitrary exceptions from user-supplied observer callbacks so one bad listener can't break notification of the others (see logging.exception call there), which is correct existing behaviour, not a bug. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> (cherry picked from commit ef542c01c2fbc199fb945e1e5d64d31be2ca2ccc)
Covers modelling/core/stratigraphic_column.py (769 lines, previously the largest untested file in the repo) - construction, add/remove/reorder of units and unconformities, lookups, groups/summaries, notifications, and serialization round-trips. Documents 3 pre-existing bugs via characterization tests: clear(basement=True) actually empties the column instead of restoring the basement; StratigraphicUnit.min()/max() are always inf for every real unit in a normally-constructed column because the basement's inf thickness poisons the cumulative thickness walk; and the thickness setter doesn't live-update cached min/max via the observer wiring. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> (cherry picked from commit a3e6aa377f5fde32f188f2d0e3872bcbd098f3cb)
…iogram helpers Replace plain for/range() loops that do per-row numpy work with vectorized numpy operations, in intrusion_support_functions.py and _svariogram.py. Each rewrite was checked against the original loop implementation on synthetic data before committing to the change (see PR description). Vectorized: - findMinDiff: pairwise abs-diff via broadcasting instead of O(n^2) Python loop - array_from_coords: reshape instead of nested loop copy - find_inout_points: argmax-based column/row lookup instead of manual scan - shortest_path (trailing loop): argmax + boolean mask instead of per-column scan - element_neighbour: single offset-array lookup instead of two range(8) loops - index_min: masked argmin (last-tie-wins) instead of dict-based scan - grid_from_array: fully vectorized index construction for all three axes - find_peaks_and_troughs: vectorized sign-change mask instead of per-index loop - find_wavelengths: vectorized pairwise midpoint averaging Left alone (documented reasoning in code/PR): sort_2_arrays (selection sort is not stable, verified duplicates reorder differently under argsort), calc_semivariogram's lag loop (already numpy-vectorized per iteration, capped at 200 iterations, vectorizing across lags would blow up memory for large datasets), and loops with early-exit/break control flow (find_inout_points' predecessor logic superseded, find_wavelengths' wl1/wl2 peak-picking loops, find_peaks_and_troughs logging loop). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> (cherry picked from commit 25db859e6afbe63fb5e75fb3df881ff0630196ba)
… and import checks
…tured grids - Introduced _face_table.py for initializing face relationships in grids. - Added _structured_grid.py for defining a 3D structured grid with properties and methods for visualization and manipulation. - Created _structured_grid_2d.py for handling 2D structured grid geometry. - Implemented _structured_grid_3d.py for 3D structured grid geometry, including rotation and volume calculations. - Developed _unstructured_mesh.py for managing unstructured mesh geometry, including nodes, elements, and neighbors. - Added utility functions in utils.py for logging, bounding box calculations, and geometric transformations. - Updated test_imports.py to verify imports for new geometry and utility modules.
- Moved the Observer, Observable, and Disposable classes from LoopStructural to loop_common, consolidating the observer pattern implementation. - Updated LoopStructural/utils/observer.py to re-export the observer classes from loop_common. - Created a new logging infrastructure in loop_common, including LogSink, StreamSink, FileSink, and SqliteSink classes for flexible logging. - Introduced timing utilities (timed_stage and timed) for structured logging of long-running operations. - Removed obsolete logging and utility classes from LoopStructural/utils.py, ensuring a cleaner codebase. - Updated ROADMAP.md to reflect the changes and improvements made during the refactor.
…g stability checks
- Deleted the _surface.py and _unstructured_mesh.py files from the geometry module. - Updated pyproject.toml to remove references to the deleted files. - Adjusted api_surface_snapshot.json to reflect changes in the Surface class import path. - Modified test__surface.py to import Surface directly from the geometry module.
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
…ntain permissions' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
…ntain permissions' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
…ntain permissions' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No major changes to API except bounding box now uses affine transformations instead of confusing local/global coordinates.