From fdf83c428c52eb66a08dfa655a82d2bb7de3e146 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 09:11:45 +1030 Subject: [PATCH 01/10] docs: make the documentation build complete - fix homepage and user-guide examples that used bounding_box.structured_grid().properties (that returns the interpolation-support grid, which has no properties dict) - fix dtm_creator importing Union from ctypes, which made the module unimportable - make docs/requirements.txt install the docs extra; add trame packages to it - start gallery examples from an empty stratigraphic column so the unused Basement group no longer logs 'Model does not contain Group_N' errors - fix reST errors (toctree, code-block directives, literal blocks, title underline) - add preparing_data and structural_frames to the user guide toctree Co-Authored-By: Claude Sonnet 5.5 --- .../modelling/core/geological_model.py | 7 ++++--- LoopStructural/utils/dtm_creator.py | 2 +- docs/requirements.txt | 9 +++------ docs/source/getting_started/about.rst | 2 +- .../getting_started/contributors_guide.rst | 3 ++- docs/source/getting_started/index.rst | 1 - docs/source/index.rst | 16 +++++++-------- docs/source/user_guide/debugging.rst | 8 ++++---- docs/source/user_guide/index.rst | 4 +++- .../user_guide/interpolation_options.rst | 6 +++--- .../user_guide/what_is_a_geological_model.rst | 20 ++++++++++++++++--- examples/1_basic/plot_2_surface_modelling.py | 2 ++ .../1_basic/plot_3_model_visualisation.py | 2 ++ .../plot_5_using_stratigraphic_column.py | 2 ++ examples/4_advanced/plot_2_using_logging.py | 2 ++ .../_finite_difference_interpolator.py | 20 ++++++++++--------- pyproject.toml | 3 +++ 17 files changed, 68 insertions(+), 41 deletions(-) diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index cabaff0cf..d8f2ece4e 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -636,12 +636,13 @@ def set_stratigraphic_column(self, stratigraphic_column, cmap="tab20"): Notes ----- - stratigraphic_column is a nested dictionary with the format - {'group': + stratigraphic_column is a nested dictionary with the format:: + + {'group': {'series1': {'min':0., 'max':10.,'id':0,'colour':} } - } + } """ warnings.warn( diff --git a/LoopStructural/utils/dtm_creator.py b/LoopStructural/utils/dtm_creator.py index 54f93663d..801980d7b 100644 --- a/LoopStructural/utils/dtm_creator.py +++ b/LoopStructural/utils/dtm_creator.py @@ -1,5 +1,5 @@ -from ctypes import Union from pathlib import Path +from typing import Union from .logging import getLogger diff --git a/docs/requirements.txt b/docs/requirements.txt index 8222675cb..0a9c23e62 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,6 +1,3 @@ -sphinx -sphinx_gallery -sphinxcontrib-bibtex -pydata_sphinx_theme -myst-parser -sphinx-design \ No newline at end of file +# Build the docs with the project and its docs extra: +# pip install -r docs/requirements.txt +.[docs] diff --git a/docs/source/getting_started/about.rst b/docs/source/getting_started/about.rst index a53811168..10897683c 100644 --- a/docs/source/getting_started/about.rst +++ b/docs/source/getting_started/about.rst @@ -11,7 +11,7 @@ The interpolation algorithms behind LoopStructural are implementations of resear Within LoopStructural we have implemented three discrete interpolation algorithms for implicit 3D modelling: * A Piecewise Linear Interpolation algorithm where the interpolation is performed on a tetraheral mesh minimising the second derivative of the implicit funtion between neighbouring tetrahedron. -This interpolation algorithm is an implementation of the discrete smooth interpolation algorithm included in Gocad-Skua, which is heavily based on the work of Frank et al., 2007 and Caumon et al., 2007. + This interpolation algorithm is an implementation of the discrete smooth interpolation algorithm included in Gocad-Skua, which is heavily based on the work of Frank et al., 2007 and Caumon et al., 2007. * This interpolation algorithm can also be framed by minimising the second derivatives using finite different on a regular cartesian grid which has been presented by Irakarama et al., diff --git a/docs/source/getting_started/contributors_guide.rst b/docs/source/getting_started/contributors_guide.rst index 4d723ac46..3fe7cf5ab 100644 --- a/docs/source/getting_started/contributors_guide.rst +++ b/docs/source/getting_started/contributors_guide.rst @@ -31,7 +31,7 @@ If you would like to contribute code to LoopStructural please open an `issue [optional scope]: [optional body] diff --git a/docs/source/getting_started/index.rst b/docs/source/getting_started/index.rst index 0a1fa2a93..e7fcdc81a 100644 --- a/docs/source/getting_started/index.rst +++ b/docs/source/getting_started/index.rst @@ -9,4 +9,3 @@ Getting Started about contributors_guide loopstructural_design - CHANGELOG diff --git a/docs/source/index.rst b/docs/source/index.rst index 9be018a78..5b4656fd2 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -52,7 +52,7 @@ evaluate the scalar field and gradient of the interpolator at some random locati model = GeologicalModel(bb[0,:],bb[1,:]) model.data = data # nelements specifies the number of discrete interpolation elements - # 'stratí' is the feature name in the data dataframe + # 'strati' is the feature name in the data dataframe model.create_and_add_foliation('strati',nelements=1e5) model.update() # get the value of the interpolator at some random locations @@ -84,6 +84,7 @@ constraints. :force_static: from LoopStructural import BoundingBox, InterpolatorBuilder + from LoopStructural.geometry import StructuredGrid from LoopStructural.utils import EuclideanTransformation from LoopStructural.datasets import load_claudius @@ -105,8 +106,13 @@ constraints. .build() ) # Set the number of elements in the bounding box to 10000 and create a structured grid + # to hold the results bounding_box.nelements = 10000 - mesh = bounding_box.structured_grid() + mesh = StructuredGrid( + origin=bounding_box.origin, + step_vector=bounding_box.step_vector, + nsteps=bounding_box.nsteps, + ) # add the interpolated values to the mesh at the nodes mesh.properties["v"] = interpolator.evaluate_value(mesh.nodes) @@ -128,12 +134,6 @@ constraints. getting_started/index _auto_examples/index user_guide/index - - - .. toctree:: - :caption: LoopStructural API - :hidden: - API diff --git a/docs/source/user_guide/debugging.rst b/docs/source/user_guide/debugging.rst index f7c20ae12..1432325b5 100644 --- a/docs/source/user_guide/debugging.rst +++ b/docs/source/user_guide/debugging.rst @@ -43,7 +43,7 @@ Interpolator is taking a long time - Check how many elements you are using, :code:`1e5` should take under 5 minutes on a modern laptop - How many faults are associated with the feature? :code:`len(feature.faults)`, you can remove the faults by setting this to an empty list -.. code-block: +.. code-block:: faults = feature.faults feature.faults = [] @@ -74,7 +74,7 @@ No displacement - Check that the displacement is large enough that it should be visible. E.g. compare displacement to model bounding box - Has the fault frame interpolated correctly? Coordinates 0 and 1 should both have values. -.. code-block: +.. code-block:: view.add_scalar_field(fault[0]) view.add_scalar_field(fault[1]) @@ -84,13 +84,13 @@ No displacement - Is the fault parallel to the feature being faulted? - Has the fault been added to the feature being faulted? -.. code-block: +.. code-block:: print([f.name for f in faulted_feature.faults]) - Is the fault displacement vector correct? Add the vector field to the visualisation -.. code-block: +.. code-block:: view = LavaVuModelViewer(model) view.add_vector_field(model['fault'][1],locations=model.regular_grid()[::200]) #random 200 locations diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst index c754775b3..fb6f2e141 100644 --- a/docs/source/user_guide/index.rst +++ b/docs/source/user_guide/index.rst @@ -11,10 +11,12 @@ for a list of topics covered in this section. what_is_a_geological_model input_data + preparing_data geological_model + structural_frames interpolation_options native_solver fold_modelling fault_modelling visualisation - debugging \ No newline at end of file + debugging diff --git a/docs/source/user_guide/interpolation_options.rst b/docs/source/user_guide/interpolation_options.rst index a8370a384..786eedf26 100644 --- a/docs/source/user_guide/interpolation_options.rst +++ b/docs/source/user_guide/interpolation_options.rst @@ -168,7 +168,7 @@ Iterative solvers: * algorithmic multigrid solver - using pyamg * lsqr - this solver uses the rectangular matrix directly, therefore does not require computing A.T A and A.T B. Uses scipy -.. code-block: +.. code-block:: model.create_and_add_foliation('my_foliation',solver='lu') model.create_and_add_foliation('my_foliation',solver='chol') @@ -185,7 +185,7 @@ You can also pass a function that solves to LoopStructural if you want to use another solver by using the `external` keyword . -.. code-block: +.. code-block:: def mysolver(A,B): from scipy.sparse.linalg import gmres @@ -195,7 +195,7 @@ to LoopStructural if you want to use another solver by using the `external` keyw The solution to the least squares problem will be stored in the interpolator object and can be easily accessed: -.. code-block: +.. code-block:: model.create_and_add_foliation('my_foliation',solver='pyamg') pyamg_solution = model['my_foliation'].interpolator.c[model['my_foliation'].interpolator.region] diff --git a/docs/source/user_guide/what_is_a_geological_model.rst b/docs/source/user_guide/what_is_a_geological_model.rst index a8e86f611..f33d7efc1 100644 --- a/docs/source/user_guide/what_is_a_geological_model.rst +++ b/docs/source/user_guide/what_is_a_geological_model.rst @@ -49,6 +49,7 @@ The interpolator can be evaluated on a grid to visualise the geometry of the geo import numpy as np import pandas as pd from LoopStructural import InterpolatorBuilder, BoundingBox + from LoopStructural.geometry import StructuredGrid from LoopStructural.visualisation import Loop3DView bounding_box = BoundingBox(np.array([-2, -2, -2]), np.array([2, 2, 2])) @@ -60,7 +61,11 @@ The interpolator can be evaluated on a grid to visualise the geometry of the geo .add_normal_constraints(vector.values).setup_interpolator() .build() ) - mesh = bounding_box.structured_grid() + mesh = StructuredGrid( + origin=bounding_box.origin, + step_vector=bounding_box.step_vector, + nsteps=bounding_box.nsteps, + ) mesh.properties['val'] = interpolator.evaluate_value(mesh.nodes) view = Loop3DView() view.add_mesh(mesh.vtk()) @@ -118,7 +123,11 @@ If we were to change the length of the vector we would change the distance betwe .add_normal_constraints(vector.values).setup_interpolator() .build() ) - mesh2 = bounding_box.structured_grid() + mesh2 = StructuredGrid( + origin=bounding_box.origin, + step_vector=bounding_box.step_vector, + nsteps=bounding_box.nsteps, + ) mesh2.properties['val'] = interpolator.evaluate_value(mesh2.nodes) view = Loop3DView() view.add_mesh(mesh2.vtk().contour([-1, 0, 1]),color='red') @@ -151,6 +160,7 @@ Following the example above we will use two lines of points with a value of 0 an import numpy as np import pandas as pd from LoopStructural import InterpolatorBuilder, BoundingBox + from LoopStructural.geometry import StructuredGrid from LoopStructural.visualisation import Loop3DView x = np.linspace(0, 1, 10) @@ -175,7 +185,11 @@ Following the example above we will use two lines of points with a value of 0 an .setup_interpolator() .build() ) - mesh = bounding_box.structured_grid() + mesh = StructuredGrid( + origin=bounding_box.origin, + step_vector=bounding_box.step_vector, + nsteps=bounding_box.nsteps, + ) mesh.properties["val"] = interpolator.evaluate_value(mesh.nodes) view = Loop3DView() view.add_mesh(mesh.vtk().contour([0, 0.5]), opacity=0.4) diff --git a/examples/1_basic/plot_2_surface_modelling.py b/examples/1_basic/plot_2_surface_modelling.py index b4c4c6ea5..c998d7766 100644 --- a/examples/1_basic/plot_2_surface_modelling.py +++ b/examples/1_basic/plot_2_surface_modelling.py @@ -98,6 +98,8 @@ # add a foliation using the `create_and_add_foliation` method. # Define stratigraphic column with scalar field ranges for each unit +# the default column ends in an empty 'Basement' group; start from an empty column +model.stratigraphic_column.clear(basement=False) vals = [0, 60, 250, 330, 600] for i in range(len(vals) - 1): model.stratigraphic_column.add_unit( diff --git a/examples/1_basic/plot_3_model_visualisation.py b/examples/1_basic/plot_3_model_visualisation.py index d27ab7c51..9eb5318a1 100644 --- a/examples/1_basic/plot_3_model_visualisation.py +++ b/examples/1_basic/plot_3_model_visualisation.py @@ -28,6 +28,8 @@ model = GeologicalModel(bb[0, :], bb[1, :]) model.data = data strati = model.create_and_add_foliation("strati",nelements=1e4) +# the default column ends in an empty 'Basement' group; start from an empty column +model.stratigraphic_column.clear(basement=False) vals = [0, 60, 250, 330, 600] for i in range(len(vals) - 1): model.stratigraphic_column.add_unit( diff --git a/examples/1_basic/plot_5_using_stratigraphic_column.py b/examples/1_basic/plot_5_using_stratigraphic_column.py index 3a6c46459..2ec2d8a6e 100644 --- a/examples/1_basic/plot_5_using_stratigraphic_column.py +++ b/examples/1_basic/plot_5_using_stratigraphic_column.py @@ -60,6 +60,8 @@ # "strati" (oldest group) - four units from shallowest to deepest, the # last of which extends to infinite thickness +# the default column ends in an empty 'Basement' group; start from an empty column +model.stratigraphic_column.clear(basement=False) model.stratigraphic_column.add_unit("unit2", thickness=60, id=1) model.stratigraphic_column.add_unit("unit3", thickness=190, id=2) model.stratigraphic_column.add_unit("unit4", thickness=80, id=3) diff --git a/examples/4_advanced/plot_2_using_logging.py b/examples/4_advanced/plot_2_using_logging.py index cb710372a..79a01a080 100644 --- a/examples/4_advanced/plot_2_using_logging.py +++ b/examples/4_advanced/plot_2_using_logging.py @@ -35,6 +35,8 @@ def build_claudius_model(): model = GeologicalModel(bb[0, :], bb[1, :]) model.data = data + # the default column ends in an empty 'Basement' group; start from an empty column + model.stratigraphic_column.clear(basement=False) vals = [0, 60, 250, 330, 600] for i in range(len(vals) - 1): model.stratigraphic_column.add_unit( diff --git a/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py index cedad6746..03ee60c39 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_finite_difference_interpolator.py @@ -142,15 +142,17 @@ def setup_interpolator(self, **kwargs): Default masks are the second derivative in x,y,z direction and the second derivative of x wrt y and y wrt z and z wrt x. Custom masks can be used by specifying the operator as a 3d numpy array - e.g. [ [ [ 0 0 0 ] - [ 0 1 0 ] - [ 0 0 0 ] ] - [ [ 1 1 1 ] - [ 1 1 1 ] - [ 1 1 1 ] ] - [ [ 0 0 0 ] - [ 0 1 0 ] - [ 0 0 0 ] ] + e.g.:: + + [ [ [ 0 0 0 ] + [ 0 1 0 ] + [ 0 0 0 ] ] + [ [ 1 1 1 ] + [ 1 1 1 ] + [ 1 1 1 ] ] + [ [ 0 0 0 ] + [ 0 1 0 ] + [ 0 0 0 ] ] ] Returns ------- diff --git a/pyproject.toml b/pyproject.toml index 8e5a7b76b..ba13e124c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,6 +66,9 @@ docs = [ "sphinxcontrib-bibtex", "myst-parser", "sphinx-design", + "trame", + "trame-vtk", + "trame-vuetify", ] [project.urls] Documentation = 'https://Loop3d.org/LoopStructural/' From 26dd57ab232bc2339e86da4b43a4670624e6aa40 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:05:55 +1030 Subject: [PATCH 02/10] docs: rewrite the visualisation guide, fix the API reference, add a docs check - visualisation.rst and debugging.rst described the removed LavaVuModelViewer; rewrite them around Loop3DView. Every snippet was run. Lambda features need an explicit bounding_box - map2loop.rst documented GeologicalModel.from_map2loop_directory, which does not exist; describe Map2LoopProcessor + GeologicalModel.from_processor (checked against the code) and add the page to the User Guide. Fix the stale docstring that still named the old method - API reference: index each class once, where it is defined. Re-exports are shown on each page but marked :no-index:; document loop_common and loop_interpolation, where most of the re-exported classes are defined; leave third-party members (pydantic, sklearn, pyvista, abc) out of the class and module pages - fix the remaining reST errors in docs, examples and two loop_interpolation docstrings; exclude the internal to-do page; drop the LavaVu notes from the contributors guide - add .github/workflows/docs-check.yml: build the docs on pull requests that touch them; a failing gallery example fails the build Sphinx warnings: 98 -> 5. Co-Authored-By: Claude Sonnet 5.5 --- .github/workflows/docs-check.yml | 69 ++++++ .../modelling/input/process_data.py | 7 +- docs/source/API.rst | 6 + docs/source/_static/overview.rst | 2 +- .../_templates/custom-class-template.rst | 5 +- .../_templates/custom-module-template.rst | 18 +- docs/source/conf.py | 68 +++++- .../getting_started/contributors_guide.rst | 11 +- docs/source/user_guide/debugging.rst | 35 +-- docs/source/user_guide/index.rst | 1 + .../user_guide/interpolation_options.rst | 2 +- docs/source/user_guide/map2loop.rst | 40 +++- docs/source/user_guide/preparing_data.rst | 5 +- docs/source/user_guide/visualisation.rst | 199 ++++++++++-------- examples/1_basic/plot_1_data_preparation.py | 1 + examples/2_fold/plot_3_p2_fold.py | 2 +- .../_fd_fold_interpolator.py | 1 + .../src/loop_interpolation/_fold_event.py | 1 + 18 files changed, 343 insertions(+), 130 deletions(-) create mode 100644 .github/workflows/docs-check.yml diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml new file mode 100644 index 000000000..4f4a23e67 --- /dev/null +++ b/.github/workflows/docs-check.yml @@ -0,0 +1,69 @@ +name: "📖 Docs check" + +# Builds the documentation on pull requests that touch it. The build fails if Sphinx errors or if +# any gallery example raises (sphinx-gallery reports all failing examples at the end). Unlike documentation.yml this does not deploy anything. +on: + pull_request: + branches: + - master + paths: + - "docs/**" + - "examples/**" + - "LoopStructural/**" + - "packages/*/src/**" + - "pyproject.toml" + - ".github/workflows/docs-check.yml" + workflow_dispatch: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + docs: + name: Build documentation + runs-on: ubuntu-latest + timeout-minutes: 45 + env: + PYVISTA_OFF_SCREEN: "True" + MPLBACKEND: Agg + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + cache: "pip" + + - name: Install system libraries for off-screen rendering + run: | + sudo apt-get update + sudo apt-get install -y libosmesa6-dev libgl1-mesa-dev python3-tk pandoc + + - name: Install the packages and the docs extra + run: | + pip install -e packages/loop_common -e packages/loop_interpolation \ + -e "packages/loopstructural_visualisation[all]" + pip install -e ".[docs]" + + - name: Use the OSMesa build of VTK + run: | + pip uninstall vtk -y + pip install vtk-osmesa==9.3.0 --index-url https://gitlab.kitware.com/api/v4/projects/13/packages/pypi/simple + + # pyvista's viewer directive writes _static as a file on a fresh output directory + - name: Build + run: | + mkdir -p docs/build/html/_static + python -m sphinx -b html docs/source docs/build/html -w sphinx_warnings.txt --keep-going + + - name: Show the warnings + if: always() + run: if [ -e sphinx_warnings.txt ]; then cat sphinx_warnings.txt; fi + + - name: Upload HTML + if: always() + uses: actions/upload-artifact@v4 + with: + name: docs-html + path: docs/build/html/ diff --git a/LoopStructural/modelling/input/process_data.py b/LoopStructural/modelling/input/process_data.py index c0d19466b..d501536dc 100644 --- a/LoopStructural/modelling/input/process_data.py +++ b/LoopStructural/modelling/input/process_data.py @@ -217,10 +217,11 @@ def stratigraphy_cmap(self, supergroup="supergroup_0"): -------- Add a isosurfaces to the model coloured by the map colours - >>> model, processor = GeologicalModel.from_map2loop_directory(directory) + >>> processor = Map2LoopProcessor(directory) + >>> model = GeologicalModel.from_processor(processor) >>> view = Loop3DView(model) - >>> view.add_isosurface(model['supergroup_0'],**processor.stratigraphy_cmap('supergroup_0')) - >>> view.interactive() + >>> view.plot_surface(model['supergroup_0'], **processor.stratigraphy_cmap('supergroup_0')) + >>> view.display() Raises diff --git a/docs/source/API.rst b/docs/source/API.rst index d6f7d2e54..ecf530eee 100644 --- a/docs/source/API.rst +++ b/docs/source/API.rst @@ -1,6 +1,10 @@ API --- +Classes are documented where they are defined. Many are re-exported by the top-level +``LoopStructural`` package and its subpackages (for example ``LoopStructural.BoundingBox`` is +``loop_common.geometry.BoundingBox``); those entries are shown on each page but indexed once. + .. autosummary:: :caption: API @@ -14,3 +18,5 @@ API LoopStructural.visualisation LoopStructural.geometry LoopStructural.utils + loop_common + loop_interpolation diff --git a/docs/source/_static/overview.rst b/docs/source/_static/overview.rst index b7af02bd5..1a147cd6d 100644 --- a/docs/source/_static/overview.rst +++ b/docs/source/_static/overview.rst @@ -1,7 +1,7 @@ LoopStructural Philosophy ========================== Loopstructural is an open source 3D geological modelling library with minimal dependencies - *numpy, scipy and scikit-learn for model building and evaluation*. -For visualisation *LavaVu and scikit-image* are required. +For visualisation *PyVista (through loopstructuralvisualisation) and scikit-image* are required. The main entry point for using LoopStructural is to create a GeologicalModel. The GeologicalModel manages the dataset, model area, creating of different elements inside the GeologicalModel. diff --git a/docs/source/_templates/custom-class-template.rst b/docs/source/_templates/custom-class-template.rst index 46cae576b..2737e22bd 100644 --- a/docs/source/_templates/custom-class-template.rst +++ b/docs/source/_templates/custom-class-template.rst @@ -3,9 +3,10 @@ .. currentmodule:: {{ module }} .. autoclass:: {{ objname }} + {% if is_reexport(module, objname) %}:no-index:{% endif %} :members: - :show-inheritance: - :inherited-members: + :show-inheritance: + :inherited-members:{% set skip = third_party_bases(module, objname) %}{% if skip %} {{ skip }}{% endif %} {% block methods %} .. automethod:: __init__ diff --git a/docs/source/_templates/custom-module-template.rst b/docs/source/_templates/custom-module-template.rst index 26f263563..ba9e74921 100644 --- a/docs/source/_templates/custom-module-template.rst +++ b/docs/source/_templates/custom-module-template.rst @@ -15,39 +15,43 @@ {% endblock %} {% block functions %} - {% if functions %} + {% set functions_ours = ours(fullname, functions) %} + {% if functions_ours %} .. rubric:: {{ _('Functions') }} .. autosummary:: :toctree: :recursive: - {% for item in functions %} + {% for item in functions_ours %} {{ item }} {%- endfor %} {% endif %} {% endblock %} {% block classes %} - {% if classes %} + {% set classes_ours = ours(fullname, classes) %} + {% if classes_ours %} .. rubric:: {{ _('Classes') }} .. autosummary:: :toctree: :template: custom-class-template.rst :recursive: - {% for item in classes %} + {% for item in classes_ours %} {{ item }} {%- endfor %} {% endif %} {% endblock %} {% block exceptions %} - {% if exceptions %} + {% set exceptions_ours = ours(fullname, exceptions) %} + {% if exceptions_ours %} .. rubric:: {{ _('Exceptions') }} .. autosummary:: - :toctree: - {% for item in exceptions %} + :toctree: + :template: custom-class-template.rst + {% for item in exceptions_ours %} {{ item }} {%- endfor %} {% endif %} diff --git a/docs/source/conf.py b/docs/source/conf.py index 318164eb7..dfcfa56f6 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -44,6 +44,68 @@ napoleon_numpy_docstring = True # False # Force consistency, leave only Google napoleon_use_rtype = False # More legible autosummary_imported_members = True + + +def _is_reexport(modname, objname): + """True when ``modname.objname`` is defined in another module and only imported here.""" + import importlib + + try: + obj = getattr(importlib.import_module(modname), objname) + except Exception: # noqa: BLE001 - anything unimportable is documented as usual + return False + return getattr(obj, "__module__", modname) != modname + + +# Classes re-exported by several packages (LoopStructural, LoopStructural.utils, ...) are +# indexed once, where they are defined. The other pages still show them, marked ``:no-index:``, +# so Sphinx does not warn about duplicate objects or ambiguous cross-references. +_OUR_PACKAGES = ("LoopStructural", "loop_common", "loop_interpolation", "loopstructuralvisualisation") + + +def _ours(modname, names): + """Keep the members of ``modname`` that belong to this project. + + ``autosummary_imported_members`` also lists whatever a module imports (pydantic, sklearn, + abc, ...), whose docstrings are not written for these docs. + """ + import importlib + + module = importlib.import_module(modname) + kept = [] + for name in names: + owner = getattr(getattr(module, name, None), "__module__", None) or modname + if owner.split(".")[0] in _OUR_PACKAGES: + kept.append(name) + return kept + + +def _third_party_bases(modname, objname): + """Comma-separated names of the base classes of ``modname.objname`` that are not ours. + + Used to leave members inherited from e.g. ``pyvista.Plotter`` or ``pydantic.BaseModel`` out of + the class pages: their docstrings link to labels that only exist in those projects' docs. + """ + import importlib + + try: + cls = getattr(importlib.import_module(modname), objname) + mro = cls.__mro__ + except Exception: # noqa: BLE001 - not a class, or not importable + return "" + names = [ + base.__name__ + for base in mro[1:] + if base is not object and base.__module__.split(".")[0] not in _OUR_PACKAGES + ] + return ", ".join(dict.fromkeys(names)) + + +autosummary_context = { + "is_reexport": _is_reexport, + "ours": _ours, + "third_party_bases": _third_party_bases, +} autosummary_ignore_module_all = False # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom @@ -72,7 +134,11 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns = [] +exclude_patterns = ["getting_started/todo.rst"] + +# sphinx-gallery and autosummary_context hold functions, which cannot be pickled into the +# environment cache; this only means those two settings are re-read on every build. +suppress_warnings = ["config.cache"] html_theme_options = { "icon_links": [ diff --git a/docs/source/getting_started/contributors_guide.rst b/docs/source/getting_started/contributors_guide.rst index 3fe7cf5ab..f9e0f90a5 100644 --- a/docs/source/getting_started/contributors_guide.rst +++ b/docs/source/getting_started/contributors_guide.rst @@ -38,15 +38,20 @@ To get started, fork and clone the repository. It's then best to get started wit pip install -r requirements.txt conda list -N.B. On Linux, the LavaVu package requires the installation of extra tools:: - - sudo apt install build-essential libgl1-mesa-dev libx11-dev zlib1g-dev +N.B. The docs and examples render 3D views with PyVista. On a headless Linux machine install the +OpenGL libraries it needs (e.g. ``sudo apt install libgl1-mesa-dev libosmesa6-dev``) and set +``PYVISTA_OFF_SCREEN=true``. For changes to ``./docs``:: pip install -r docs/requirements.txt make -C ./docs html +``docs/requirements.txt`` installs the project with its ``docs`` extra. The build runs every +example in ``examples/`` and takes several minutes; an example that raises an error fails the +check in CI. If the first build stops with ``NotADirectoryError`` for ``_static``, create the +directory first with ``mkdir -p docs/build/html/_static`` and build again. + Building the docs the first time takes a little while. After this, you will see a build directory. From here you can preview the generated HTML in your browser. diff --git a/docs/source/user_guide/debugging.rst b/docs/source/user_guide/debugging.rst index 1432325b5..d90569858 100644 --- a/docs/source/user_guide/debugging.rst +++ b/docs/source/user_guide/debugging.rst @@ -7,8 +7,8 @@ Surfaces are not displaying ~~~~~~~~~~~~~~~~~~~~~~~~~~~ - check to see what the scalar field looks linked - look at the values of the nodes of the implicit function supports `feature.interpolator.c` if all of the nodes are 0 the interpolator has not solved. -- Does the scalar field show different colours? `view.add_scalar_field(feature)` -- If the scalar field shows the right geometry, are the surface values correct? +- Does the scalar field show different colours? `view.plot_scalar_field(feature)` +- If the scalar field shows the right geometry, are the surface values correct? Compare the range of the field with the ``value`` you pass to `view.plot_surface(feature, value=...)` Surface geometry is erroneous @@ -19,12 +19,13 @@ Surface geometry is erroneous - Try using only value constraints - Change gradient norm constraints to gradient direction constraints (rename :code:`nx,ny,nz` to :code:`gx,gy,gz`) - Try removing some of the data points -- Try changing to using interface constraints instead of value, to do this rename column :code:`val` to code:`interface` but only use :code:`FDI`. -Alternatively use surfe with `method=interface`. Using either of these methods the value of the scalar field will change. +- Try changing to using interface constraints instead of value, to do this rename column :code:`val` to code:`interface` but only use :code:`FDI`. + Alternatively use surfe with `method=interface`. Using either of these methods the value of the scalar field will change. - The tolerance of the conjugate gradient solver is scaled to the size of the model assuming that the scalar field is interpolated to have a magnitude of the gradient -norm to be close to 1. The tolerance is saved in `feature.builder.build_arguments['tol']` but can be overwritten using the argument in code:`create_and_add` `tol=val`. -Smaller values mean the solution will be a closer fit to the data + regularisation. + norm to be close to 1. The tolerance is saved in `feature.builder.build_arguments['tol']` but can be overwritten using the argument in code:`create_and_add` `tol=val`. + Smaller values mean the solution will be a closer fit to the data + regularisation. - Try changing the relative weighting of the data points and regularisation. + * :code:`regularisation` - changes the weighting of the constant gradient or the finite difference second derivative. 1.0 is equivalent to 0.1 constant gradient. * :code:`npw=1.0` changes the weighting of the gradient norm constraints * :code:`gpw=1.0` changes the weighting of gradient direction constraints @@ -62,9 +63,9 @@ S-Plot doesn't fit the observations - Is the wavelength obvious in the s-variogram? If it is, check the guessed wavelength is reasonable Check the s - Try specifying the wavelength using :code:`limb_wl` or :code:`axis_wl` - If the data points in the S-Plot do not look correct, check the data points in the model by plotting the axial foliation and the data. The fold geometry -should be visible. + should be visible. - If the fold is non-cylindrical, try building a cylindrical fold by removing some data -- Check the polarity of the datapoints, the fold modelling requires the polarity of the vectors is constant. Check uising vector plots :code:`view.add_data(folded_foliation,vectors=True)` +- Check the polarity of the datapoints, the fold modelling requires the polarity of the vectors is constant. Check using vector plots :code:`view.plot_data(folded_foliation, vector=True)` Faults ------ @@ -76,11 +77,11 @@ No displacement .. code-block:: - view.add_scalar_field(fault[0]) - view.add_scalar_field(fault[1]) - view.add_scalar_field(fault[2]) + view.plot_scalar_field(fault[0]) + view.plot_scalar_field(fault[1]) + view.plot_scalar_field(fault[2]) -- Increase the resolution of the visualisation mesh :code:`view.nelements=1e6` +- Increase the resolution of the visualisation grid :code:`model.bounding_box.nelements = 1e6` - Is the fault parallel to the feature being faulted? - Has the fault been added to the feature being faulted? @@ -92,10 +93,12 @@ No displacement .. code-block:: - view = LavaVuModelViewer(model) - view.add_vector_field(model['fault'][1],locations=model.regular_grid()[::200]) #random 200 locations - view.add_isosurface(model['fault'][0],0) - view.interactive() + from LoopStructural.visualisation import Loop3DView + + view = Loop3DView(model) + view.plot_vector_field(model['fault'][1], scale=50, name='slip direction') + view.plot_surface(model['fault'][0], value=0, name='fault surface') + view.display() diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst index fb6f2e141..84acae3d8 100644 --- a/docs/source/user_guide/index.rst +++ b/docs/source/user_guide/index.rst @@ -20,3 +20,4 @@ for a list of topics covered in this section. fault_modelling visualisation debugging + map2loop diff --git a/docs/source/user_guide/interpolation_options.rst b/docs/source/user_guide/interpolation_options.rst index 786eedf26..6ec43f71e 100644 --- a/docs/source/user_guide/interpolation_options.rst +++ b/docs/source/user_guide/interpolation_options.rst @@ -26,7 +26,7 @@ Geological observations including the location of contacts and orientation of st .. math:: \nabla f(x,y,z) \cdot \textbf{t} = 0 -**To solve the implicit system there needs to be a minimum of two different value observations, or a value observation and a gradient normal constraint. ** +**To solve the implicit system there needs to be a minimum of two different value observations, or a value observation and a gradient normal constraint.** LoopStructural has two types of interpolation algorithms that can be used. 1. Discrete interpolation where the implicit function is approximated by basis functions on a predefined support diff --git a/docs/source/user_guide/map2loop.rst b/docs/source/user_guide/map2loop.rst index 74094cb58..ae7d44745 100644 --- a/docs/source/user_guide/map2loop.rst +++ b/docs/source/user_guide/map2loop.rst @@ -1,11 +1,29 @@ -The class method (GeologicalModel.from_map2loop_directory(m2l_directory, **kwargs)) creates an instance of a GeologicalModel from a root map2loop output directory. -There are two steps to this function the first is a pre-processing step that prepares the extracted datasets for LoopStructural. This involves combining all observations (fault surface location, fault orientation, stratigraphic orientation and stratigraphic contacts) into the one pandas DataFrame. The cumulative thickness is calculated by iterating over the stratigraphic column for each group and starting with the base of the group as 0 and adding the thickness of the lower unit. This function returns a python dictionary which is structured as shown in Table 5. The processed data dictionary can then be passed to a utility function that: -1. initialises a GeologicalModel for the specified bounding box (map area chosen for map2loop and depth specified). -2. associates the pandas DataFrame to the model -3. iterate through the fault-stratigraphy table -• adding stratigraphic packages that are associated with no faults -• add faults in reverse order, fault intersections are manually, and fault geometries are only interpolated within their object aligned bounding box to reduce computational to ensure that -• The base of each group is assigned as an unconformity for any older stratigraphic units or faults that occur below it -4. associate stratigraphic column to geological model - -A parameters dictionary can be passed to the GeologicalModel.from_map2loop_directory function which allows for the interpolator parameters (listed in Table 2) for the model to be specified. This can be passed using two python dictionaries: foliation_params contains the parameters for modelling any stratigraphic groups and fault_params contains parameters for the faults. +Using map2loop output +===================== + +`map2loop `_ extracts the observations needed for a 3D model +(stratigraphic contacts, orientations, fault traces and their relationships) from a geological map. +LoopStructural turns that output into a model in two steps: a *processor* gathers the observations +and the order of the geological features, and the model is built from the processor. + +.. code-block:: python + + from LoopStructural import GeologicalModel + from LoopStructural.modelling.input import Map2LoopProcessor + + processor = Map2LoopProcessor("path/to/map2loop/project") + model = GeologicalModel.from_processor(processor) + model.update() # run the interpolation + +:class:`~LoopStructural.modelling.input.Map2LoopProcessor` reads the directory layout written by +the original map2loop (the ``output`` and ``tmp`` folders). For other sources, or to control the +inputs yourself, build a :class:`~LoopStructural.modelling.input.ProcessInputData` from data frames +and use :meth:`GeologicalModel.from_processor` in the same way; see +:doc:`/_auto_examples/4_advanced/plot_1_model_from_geological_map` for a worked example. + +``from_processor`` creates the model for the extent given by the processor (``origin`` to ``maximum``) +and attaches its data. It then adds the faults of the processor's fault network, including splay and +abutting relationships between faults, and adds each stratigraphic group as a foliation affected by the +faults listed for it in the fault-stratigraphy table. The base of each interpolated group is added as an +unconformity, and the processor's stratigraphic column is attached to the model. The model is returned +without having been solved; call ``model.update()``. diff --git a/docs/source/user_guide/preparing_data.rst b/docs/source/user_guide/preparing_data.rst index 8c7f25379..4b783b88a 100644 --- a/docs/source/user_guide/preparing_data.rst +++ b/docs/source/user_guide/preparing_data.rst @@ -56,8 +56,9 @@ The lists are to be formated so that the first level of the list contains the yo The next level contains the order of the formations within the conformable group, with them order youngest to oldest. In the example below, there are two groups the youngest group contains a,b,c and the older group contains d,e. -.. code-block:: - [[a,b,c],['d','e]] +.. code-block:: python + + [['a', 'b', 'c'], ['d', 'e']] Note: If the stratigraphy can be interpreted as conformable, it is advisable to include all stratigraphy inside a single group. diff --git a/docs/source/user_guide/visualisation.rst b/docs/source/user_guide/visualisation.rst index 56661aac4..fbc92c480 100644 --- a/docs/source/user_guide/visualisation.rst +++ b/docs/source/user_guide/visualisation.rst @@ -1,136 +1,171 @@ Visualisation -================= -LoopStructural uses the powerful lavaVu visualisaton library for visualising the 3D models as surfaces and scalar fields. -LavaVu provides an interative viewer inside jupyter notebooks and can use matplotlib colourmaps as well as export the model to webgl. -The LavaVuModelViewer provides an easy to use interfaces for adding GeologicalFeatures to the viewer, the base lavaVu object can still be used. -Some of the functions in the viewer require a bounding box to be able to draw isosurfaces or scalar fields. -By default the viewer is initialised with a GeologicalModel and the origin and maximum of the model are used to set the viewers extents. -It is also possible to modify this manually by changing the bounding_box attribute. +============= -Visualising GeologicalFeatures ------------------------------- +LoopStructural visualises models with :class:`LoopStructural.visualisation.Loop3DView`, a thin +subclass of the `PyVista `_ ``Plotter`` that knows how to draw +geological features, faults, data and whole models. It comes from the separate +``loopstructuralvisualisation`` package; install it with -There are four main ways a GeologicalFeature can be visualised. -1. Scalar field showing the value of the implicit function on a block representing the model domain +.. code-block:: bash -.. code-block:: + pip install loopstructural[visualisation] - view.plot_scalar_field(geological_feature) +Because ``Loop3DView`` *is* a PyVista plotter, everything PyVista can do (``add_mesh``, +``add_points``, ``add_arrows``, ``camera_position``, ``screenshot``, ...) is available as well. -2. Surfaces extracted from isovalues of the implicit function +.. code-block:: python -.. code-block:: + from LoopStructural.visualisation import Loop3DView - view.plot_surface(geological_feature,value=0) + view = Loop3DView(model) # the model is used for bounds and defaults + view.plot_surface(feature, value=0) + view.display() # opens the viewer (inline in a notebook) -3. Vector field representing the gradient of the implicit function +Some methods need a bounding box to know where to evaluate the feature. By default this is +the bounding box of the model passed to ``Loop3DView``; every ``plot_*`` method also accepts a +``bounding_box=`` argument to override it. -.. code-block:: +Visualising geological features +------------------------------- - view.plot_surface(geological_feature) +There are four main ways to look at a geological feature. -4. Cross section showing the value of the implicit function at a particular location inside the model +1. A scalar field: the value of the implicit function on a volume. -.. code-block:: + .. code-block:: python - view.add_section(geological_feature,axis='x',value=0) + view.plot_scalar_field(feature, cmap="viridis") -When adding to the viewer there are a number of common arguments that can be specified to change the behaviour: +2. Surfaces: isosurfaces of the implicit function. ``value`` is one isovalue, a list of + isovalues, or an integer number of regularly spaced surfaces. - .. list-table:: Viewer kwargs - :widths: 25 75 - :header-rows: 1 + .. code-block:: python - * - Parameter - - Options - * - vmin - - minimum value for colourmap - * - vmax - - maximum value for colourmap - * - cmap - - matplotlib colourmap or string with name of matplotlib colourmap e.g. 'rainbow' - * - pyvista_kwargs - - any arguments to be passed to the :meth:`pyvista.Plotter.add_mesh` method - -The default behaviour of the viewer can also be modified. -To change the resolution of the scalar field or isosurfacing, either set the number of steps in x,y,z or the total number of elements of tht bounding box + view.plot_surface(feature, value=[0, 100, 300], cmap="rainbow") -.. code-block:: - - model.bounding_box.nsteps = np.array([nx,ny,nz]) - print(model.bounding_box.nsteps) - model.bounding_box.nelements = 1e6 - print(model.bounding_box.nsteps) +3. A vector field: the gradient of the implicit function, drawn as arrows. + .. code-block:: python - -Visualising data ----------------- + view.plot_vector_field(feature, scale=50) + +4. Sliders: pass ``slicer=True`` to ``plot_scalar_field`` or ``plot_surface`` to add an + interactive slicing widget. -The input data for the geological model can also be visualised within the viewer. -To add all of the data attached to a GeologicalFeature +The methods share a few arguments: -.. code-block:: +.. list-table:: Common arguments + :widths: 25 75 + :header-rows: 1 - view.plot_data(geological_feature,value=True,vector=True,scale=10,pyvista_kwargs={}) + * - Parameter + - Meaning + * - ``vmin``, ``vmax`` + - limits of the colour map + * - ``cmap`` + - a matplotlib colour map, or its name e.g. ``"rainbow"`` + * - ``opacity`` + - transparency of the object, between 0 and 1 + * - ``name`` + - name of the object in the viewer. Adding two objects with the same name replaces the first + * - ``pyvista_kwargs`` + - a dictionary passed on to :meth:`pyvista.Plotter.add_mesh` + * - ``bounding_box`` + - the region to evaluate the feature in, by default the model's bounding box -Alternatively, points can be added into the model using the: +Resolution +~~~~~~~~~~ -.. code-block:: +The scalar field and isosurfaces are evaluated on a regular grid over the bounding box. +To change its resolution either set the number of points in each direction or the approximate +number of cells: - view.add_points(xyz) +.. code-block:: python -or points with an associated value: + model.bounding_box.nsteps = np.array([nx, ny, nz]) + print(model.bounding_box.nsteps) + + model.bounding_box.nelements = 1e6 # updates nsteps for roughly cubic cells + print(model.bounding_box.nsteps) -.. code-block:: - view.add_points(xyz,scalars=np.array(values)) +Visualising data +---------------- -or points with an associated vector: +The data attached to a feature can be drawn with the feature: -.. code-block:: +.. code-block:: python - view.add_arrows(xyz,direction) + view.plot_data(feature, value=True, vector=True, scale=50) -**Warning, you cannot have two objects in the viewer with the same name** +``value`` draws the value constraints as points and ``vector`` draws the orientation +constraints as arrows. Arbitrary points and arrows can be added with the PyVista methods: -Visualising geological model ----------------------------- +.. code-block:: python -Where the GeologicalModel has an assocaited stratigraphic column including a list of faults in the model, this can be used for automatically visualising the model. -To add a block model showing the stratigraphic unit distribution, where the value is linked to the *id* in the stratigraphic column. + view.add_points(xyz, name="my points") + view.add_points(xyz, scalars=values, name="my values") + view.add_arrows(xyz, directions, name="my vectors") -.. code-block:: +.. warning:: + + Pass a distinct ``name`` for each object. Reusing a name replaces the earlier object. + +Visualising the geological model +-------------------------------- + +When a model has a stratigraphic column (and, optionally, faults), the whole model can be +drawn at once. A block model colours each cell by the ``id`` of its stratigraphic unit: + +.. code-block:: python view.plot_block_model() -alternatively, the surfaces can be added to the visualisation. +or draw a surface for each stratigraphic unit and, optionally, each fault: + +.. code-block:: python -.. code-block:: + view.plot_model_surfaces(strati=True, faults=True) - view.plot_model_surfaces(faults=True,strati=True) +Faults can also be drawn individually, with the slip vector and displacement volume: -Where the faults flag allows you to disable visualising the fault surface. -This function will add a surface for every stratigraphic unit in the stratigraphic column. +.. code-block:: python + + view.plot_fault(model["fault_name"], surface=True, slip_vector=True) Custom visualisation -------------------- -It is also possible to pass custom functions to the viewer using the LambdaGeologicalFeature. -The LambdaGeologicalFeature is a class allowing for a function describing the value at locations or vector field. +A :class:`LoopStructural.modelling.features.LambdaGeologicalFeature` wraps any function of +location, so arbitrary quantities can be drawn with the same methods. -.. code-block:: +.. code-block:: python from LoopStructural.modelling.features import LambdaGeologicalFeature + def x_function(xyz): - return xyz[:,0] - custom_feature = LambdaGeologicalFeature(x_function,name='x feature') - view.plot_surface(custom_feature,value) - view.plot_scalar_field(custom_feature) + return xyz[:, 0] + + custom_feature = LambdaGeologicalFeature(x_function, name="x feature") + + # a lambda feature has no model, so say where to evaluate it + bb = model.bounding_box + view.plot_surface(custom_feature, value=0.5 * (bb.origin[0] + bb.maximum[0]), bounding_box=bb) + view.plot_scalar_field(custom_feature, bounding_box=bb) -The function passed to the LambdaGeologicalFeature can be as simple or complicated as required. -It will be evaluated for the locations within the model that the visualisation requires, usually between the origin and maximum of the geoloical model. -By default the min and max values are 0, however these can be overwritten by setting the attribute min and max to the required values. +The function is evaluated at the locations the visualisation needs, normally between the +origin and maximum of the bounding box. Set the ``min`` and ``max`` attributes of the feature +if the colour range should differ from the default. +Using the grid directly +----------------------- +``feature.scalar_field()`` evaluates a feature on the bounding box's regular grid and returns a +:class:`LoopStructural.geometry.StructuredGrid`. Its ``properties`` hold the values at the nodes and +``cell_properties`` hold values at the cell centres. ``vtk()`` returns a PyVista +grid and ``save()`` writes it to VTK, GeoH5, OMF or GOCAD files. +.. code-block:: python + grid = feature.scalar_field() + grid.plot() # quick look with PyVista + grid.save("scalar_field.vtk") diff --git a/examples/1_basic/plot_1_data_preparation.py b/examples/1_basic/plot_1_data_preparation.py index fcbe84bbd..17fbf0dcc 100644 --- a/examples/1_basic/plot_1_data_preparation.py +++ b/examples/1_basic/plot_1_data_preparation.py @@ -88,6 +88,7 @@ # * :code:`gx,gy,gz` represent a constraint where the interpolated scalar field is parallel to this vector # * :code:`nx,ny,nz` represent a constraint which set the partial derivatives of the function. # * :code:`feature_name` assigns which geologicalfeature the observations control +# # **Note** for the interpolator to solve there needs to be two unique values or a norm constraint # for the interpolator to be able to find a solution. diff --git a/examples/2_fold/plot_3_p2_fold.py b/examples/2_fold/plot_3_p2_fold.py index 7d32a28a3..e0c19f5c1 100644 --- a/examples/2_fold/plot_3_p2_fold.py +++ b/examples/2_fold/plot_3_p2_fold.py @@ -1,6 +1,6 @@ """ 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 diff --git a/packages/loop_interpolation/src/loop_interpolation/_fd_fold_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_fd_fold_interpolator.py index 003616a3f..2757814cf 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_fd_fold_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_fd_fold_interpolator.py @@ -73,6 +73,7 @@ def setup_interpolator(self, **kwargs): :meth:`add_fold_constraints`. Keys are the same as the parameters of that method (e.g. ``fold_orientation``, ``fold_regularisation``, …). + All other kwargs are forwarded to the parent :meth:`FiniteDifferenceInterpolator.setup_interpolator`. """ diff --git a/packages/loop_interpolation/src/loop_interpolation/_fold_event.py b/packages/loop_interpolation/src/loop_interpolation/_fold_event.py index f7a8b9691..b8c4d9e5e 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_fold_event.py +++ b/packages/loop_interpolation/src/loop_interpolation/_fold_event.py @@ -16,6 +16,7 @@ class FoldEvent: """Describes fold geometry via a fold frame and two rotation-angle functions. The fold frame provides two coordinate fields: + - ``foldframe.features[0]``: axial-surface scalar field (gx). Its gradient (``evaluate_gradient``) gives the across-fold direction and its value (``evaluate_value``) parameterises the fold limb rotation angle. From 230cfa3a491c7728cb1a45c542137e3fa4207759 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:18:12 +1030 Subject: [PATCH 03/10] docs: add an Interpolation gallery section built on the golden-test models Three examples that measure interpolator behaviour against synthetic models with a known scalar field, the idea behind packages/loop_interpolation/tests/golden: - 5a choosing an interpolator: FDI, PLI and P2 on a gaussian dome (error, crest error, data misfit, gradient angle, time) - 5b regularisation: sweep the weight on exact and noisy data; the best weight moves from the smallest to about the default as noise increases - 5c fold constraints: three plain interpolators against FDFold, DFI and P2Fold with two contacts and one orientation on a cylindrical fold The prose in each was written from the printed results and the rendered figures. Co-Authored-By: Claude Sonnet 5.5 --- examples/5_interpolation/README.rst | 7 + .../plot_1_choosing_an_interpolator.py | 208 +++++++++++++++ .../5_interpolation/plot_2_regularisation.py | 209 ++++++++++++++++ .../plot_3_fold_constraints.py | 236 ++++++++++++++++++ examples/README.rst | 3 + 5 files changed, 663 insertions(+) create mode 100644 examples/5_interpolation/README.rst create mode 100644 examples/5_interpolation/plot_1_choosing_an_interpolator.py create mode 100644 examples/5_interpolation/plot_2_regularisation.py create mode 100644 examples/5_interpolation/plot_3_fold_constraints.py diff --git a/examples/5_interpolation/README.rst b/examples/5_interpolation/README.rst new file mode 100644 index 000000000..62abebd2d --- /dev/null +++ b/examples/5_interpolation/README.rst @@ -0,0 +1,7 @@ +Interpolation +------------- + +These examples use synthetic geological models whose scalar field is known exactly, so every +result can be measured against the truth. They show how the interpolators behave: which to +choose, what the regularisation weight does, and what fold constraints add when data are sparse. +They use the interpolators directly rather than through ``GeologicalModel``. diff --git a/examples/5_interpolation/plot_1_choosing_an_interpolator.py b/examples/5_interpolation/plot_1_choosing_an_interpolator.py new file mode 100644 index 000000000..419251a8b --- /dev/null +++ b/examples/5_interpolation/plot_1_choosing_an_interpolator.py @@ -0,0 +1,208 @@ +""" +5a. Choosing an interpolator +============================ + +LoopStructural can solve for the same implicit function on three kinds of support: + +* ``FiniteDifferenceInterpolator`` (FDI) on a regular grid, +* ``P1Interpolator`` (PLI) on a mesh of linear tetrahedra, +* ``P2Interpolator`` on a mesh of quadratic tetrahedra. + +They minimise the same kind of misfit (fit the data, keep the second derivative small) but +discretise it differently, so the same data give slightly different surfaces. To see how much +that matters we need to know the right answer, so this example builds a synthetic geological +model with a known scalar field: a gaussian dome. The data are sampled from it, then each +interpolator is asked to recover it, and we measure the error. + +The same comparison, on several synthetic models, is what the interpolation regression tests +(``packages/loop_interpolation/tests/golden``) run. +""" + +import time + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy.optimize import brentq + +from LoopStructural.interpolators import ( + FiniteDifferenceInterpolator, + P1Interpolator, + P2Interpolator, + P2UnstructuredTetMesh, + StructuredGrid, + TetMesh, +) + +# %% +# A synthetic model with a known answer +# ------------------------------------- +# The stratigraphic scalar field is the height above a gaussian dome, ``s = z - dome(x, y)``. +# Horizons are the surfaces ``s = 0.3, 0.5, 0.7``. Everything lives in the unit cube. + +RELIEF, CENTRE, WIDTH = 0.2, 0.5, 0.2 +HORIZONS = (0.3, 0.5, 0.7) + + +def truth(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + return xyz[:, 2] - RELIEF * np.exp(-r2 / (2 * WIDTH**2)) + + +def truth_gradient(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + e = RELIEF * np.exp(-r2 / (2 * WIDTH**2)) / WIDTH**2 + return np.c_[e * (xyz[:, 0] - CENTRE), e * (xyz[:, 1] - CENTRE), np.ones(len(xyz))] + + +# %% +# Sample data from it the way a geologist would: contact points where nine boreholes cross +# each horizon, and five orientation measurements. + + +def borehole_contacts(xy): + points = [] + for x, y in xy: + for level in HORIZONS: + + def misfit(z, x=x, y=y, level=level): + return truth(np.array([[x, y, z]]))[0] - level + + if misfit(0.0) * misfit(1.0) < 0: + points.append([x, y, brentq(misfit, 0.0, 1.0), level]) + return np.array(points) + + +rng = np.random.default_rng(4) +borehole_xy = np.vstack( + [[0.5, 0.5], [0.35, 0.5], [0.65, 0.5], [0.5, 0.35], [0.5, 0.65], rng.uniform(0.1, 0.9, (4, 2))] +) +values = borehole_contacts(borehole_xy) # x, y, z, value +orientation_xyz = np.array( + [[0.4, 0.4, 0.5], [0.6, 0.6, 0.5], [0.3, 0.7, 0.6], [0.8, 0.3, 0.4], [0.5, 0.5, 0.5]] +) +normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] # x, y, z, nx, ny, nz +print(f"{len(values)} contact points and {len(normals)} orientations") + +# %% +# Three supports +# -------------- +# FDI and PLI use a grid of 16 cells (17 nodes) per side. The quadratic mesh has far more +# unknowns per cell, so 8 cells per side gives a similar number of unknowns. + +CELLS = 16 +origin = np.zeros(3) + + +def make_interpolator(kind): + if kind == "FDI": + grid = StructuredGrid( + origin=origin, nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + return FiniteDifferenceInterpolator(grid) + if kind == "PLI": + mesh = TetMesh( + origin=origin, nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + return P1Interpolator(mesh) + mesh = P2UnstructuredTetMesh( + origin=origin, nsteps_cells=np.full(3, 8), step_vector=np.full(3, 1 / 8) + ) + return P2Interpolator(mesh) + + +def interpolate(kind): + interpolator = make_interpolator(kind) + interpolator.set_value_constraints(values) + interpolator.set_normal_constraints(normals) + start = time.perf_counter() + interpolator.setup_interpolator() + interpolator.solve_system("cg") + return interpolator, time.perf_counter() - start + + +results = {kind: interpolate(kind) for kind in ("FDI", "PLI", "P2")} +interpolators = {kind: interpolator for kind, (interpolator, _) in results.items()} + +# %% +# How good is each one? +# --------------------- +# Evaluate on a regular set of points well inside the cube and compare with the truth: +# +# * ``rms error`` of the scalar field over the whole interior, +# * ``crest error``, the same near the top of the dome, where smoothing hurts most, +# * ``data misfit``, how closely the surface passes through the contact points, +# * ``median gradient angle``, how far the gradient direction is from the true one, +# * ``time``, to build and solve the system (indicative only, it depends on the machine). + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T +crest = points[(np.abs(points[:, 0] - CENTRE) < 0.15) & (np.abs(points[:, 1] - CENTRE) < 0.15)] + + +def angle_degrees(a, b): + cosine = np.sum(a * b, axis=1) / np.linalg.norm(a, axis=1) / np.linalg.norm(b, axis=1) + return np.rad2deg(np.arccos(np.clip(cosine, -1, 1))) + + +rows = {} +for kind, interpolator in interpolators.items(): + error = interpolator.evaluate_value(points) - truth(points) + crest_error = interpolator.evaluate_value(crest) - truth(crest) + misfit = interpolator.evaluate_value(values[:, :3]) - values[:, 3] + angle = angle_degrees(interpolator.evaluate_gradient(points), truth_gradient(points)) + rows[kind] = { + "unknowns": interpolator.dof, + "rms error": np.sqrt(np.mean(error**2)), + "crest error": np.sqrt(np.mean(crest_error**2)), + "data misfit": np.sqrt(np.mean(misfit**2)), + "median gradient angle": np.median(angle), + "time (s)": results[kind][1], + } +print(pd.DataFrame(rows).T.astype({"unknowns": int}).round(4).to_string()) + +# %% +# Looking at a section +# -------------------- +# The section through the middle of the dome (``y = 0.5``). The dashed lines are the true +# horizons and the solid lines are each interpolator's. + +x = np.linspace(0, 1, 120) +z = np.linspace(0, 1, 120) +xx, zz = np.meshgrid(x, z) +section = np.c_[xx.ravel(), np.full(xx.size, 0.5), zz.ravel()] +true_field = truth(section).reshape(xx.shape) +near_section = np.abs(values[:, 1] - 0.5) < 0.2 + +fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True) +for ax, (kind, interpolator) in zip(axes, interpolators.items()): + field = interpolator.evaluate_value(section).reshape(xx.shape) + ax.contour(xx, zz, true_field, levels=HORIZONS, colors="k", linestyles="--", linewidths=1) + ax.contour(xx, zz, field, levels=HORIZONS, colors="tab:red", linewidths=1.5) + ax.plot(values[near_section, 0], values[near_section, 2], "ko", ms=4) + ax.set_title(f"{kind}: rms error {rows[kind]['rms error']:.3f}") + ax.set_xlabel("x") + ax.set_aspect("equal") +axes[0].set_ylabel("z") +fig.suptitle("True horizons (dashed) and interpolated horizons (red), section y = 0.5") +fig.tight_layout() +plt.show() + +# %% +# What to take from this +# ---------------------- +# No interpolator wins on every measure here: +# +# * Over the whole interior P2 has the lowest error and the lowest gradient-angle error, then +# PLI, then FDI. +# * Near the crest the order changes: FDI has the lowest error and PLI the highest. +# * P2 passes through the contact points (misfit about zero) and FDI has the largest misfit, +# so FDI is the one that trades the data for smoothness. +# * The quadratic mesh is the slowest to solve, FDI the fastest. FDI is also the default +# used by ``GeologicalModel`` (``interpolatortype="FDI"``). +# +# These results come from one synthetic model, one random draw of data and one resolution, +# so treat the ranking as a property of this example, not a rule. Errors of a few percent of +# the box height are normal for sparse data: the interpolator is inventing the surface between +# the data. More data, or a finer support, reduce them. The next example looks at the other +# lever, the regularisation weight. diff --git a/examples/5_interpolation/plot_2_regularisation.py b/examples/5_interpolation/plot_2_regularisation.py new file mode 100644 index 000000000..7190c62a0 --- /dev/null +++ b/examples/5_interpolation/plot_2_regularisation.py @@ -0,0 +1,209 @@ +""" +5b. Regularisation: smoothing versus fitting +============================================ + +An implicit surface has to do two things at once: pass through the data, and behave sensibly +where there are no data. The interpolators do both by minimising a weighted sum of terms, + +* data terms, which penalise the distance from the contact points and the angle to the + orientation measurements, and +* a regularisation term, which penalises the second derivative of the scalar field, so that + the surface is as smooth as the data allow. + +The ``regularisation`` weight sets how much the smoothness term counts. Too small and the +surface chases every error in the data; too large and it flattens real structure. This example +measures that trade-off on a synthetic dome whose true shape is known, with clean and with noisy +data. For noise-free data the best weight is the smallest one that still solves, while for noisy +data the best weight is larger. +""" + +import matplotlib.pyplot as plt +import numpy as np +from scipy.optimize import brentq + +from LoopStructural.interpolators import FiniteDifferenceInterpolator, StructuredGrid + +# %% +# The synthetic model and its data +# -------------------------------- +# A gaussian dome in the unit cube, ``s = z - dome(x, y)``, with horizons at ``s = 0.3, 0.5, 0.7`` +# (the same model as the previous example). + +RELIEF, CENTRE, WIDTH = 0.2, 0.5, 0.2 +HORIZONS = (0.3, 0.5, 0.7) + + +def truth(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + return xyz[:, 2] - RELIEF * np.exp(-r2 / (2 * WIDTH**2)) + + +def truth_gradient(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + e = RELIEF * np.exp(-r2 / (2 * WIDTH**2)) / WIDTH**2 + return np.c_[e * (xyz[:, 0] - CENTRE), e * (xyz[:, 1] - CENTRE), np.ones(len(xyz))] + + +def borehole_contacts(xy): + points = [] + for x, y in xy: + for level in HORIZONS: + + def misfit(z, x=x, y=y, level=level): + return truth(np.array([[x, y, z]]))[0] - level + + if misfit(0.0) * misfit(1.0) < 0: + points.append([x, y, brentq(misfit, 0.0, 1.0), level]) + return np.array(points) + + +rng = np.random.default_rng(4) +borehole_xy = np.vstack( + [[0.5, 0.5], [0.35, 0.5], [0.65, 0.5], [0.5, 0.35], [0.5, 0.65], rng.uniform(0.1, 0.9, (4, 2))] +) +exact_values = borehole_contacts(borehole_xy) +orientation_xyz = np.array( + [[0.4, 0.4, 0.5], [0.6, 0.6, 0.5], [0.3, 0.7, 0.6], [0.8, 0.3, 0.4], [0.5, 0.5, 0.5]] +) +exact_normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] + + +def add_noise(noise, seed): + """Perturb the contacts' height by ``noise`` and the orientation vectors by ``2.5 * noise``.""" + generator = np.random.default_rng(seed) + values = exact_values.copy() + values[:, 2] += generator.normal(0, noise, len(values)) + normals = exact_normals.copy() + normals[:, 3:] += generator.normal(0, 2.5 * noise, (len(normals), 3)) + return values, normals + + +# %% +# Solving with a given regularisation weight +# ------------------------------------------ +# ``regularisation`` is passed to ``setup_interpolator``. Its default is 1, which +# gives each second-derivative term a weight of 1 against data weights of 1 for contacts and 10 +# for orientation vectors. We use a grid of 16 cells per side. + + +def solve(values, normals, regularisation): + grid = StructuredGrid( + origin=np.zeros(3), nsteps_cells=np.full(3, 16), step_vector=np.full(3, 1 / 16) + ) + interpolator = FiniteDifferenceInterpolator(grid) + interpolator.set_value_constraints(values) + interpolator.set_normal_constraints(normals) + interpolator.setup_interpolator(regularisation=regularisation) + solved = interpolator.solve_system("cg") + return interpolator, bool(solved) + + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T + + +def errors(interpolator, values): + """RMS error against the truth over the interior, and RMS distance from the contact points. + + Noise can move a contact outside the model, where there is no value: those are skipped. + """ + field_error = interpolator.evaluate_value(points) - truth(points) + data_error = interpolator.evaluate_value(values[:, :3]) - values[:, 3] + return np.sqrt(np.mean(field_error**2)), np.sqrt(np.nanmean(data_error**2)) + + +# %% +# Sweeping the weight +# ------------------- +# For each weight, solve with exact data and with noisy data (height noise 0.1, a tenth of the +# box, averaged over five random draws so that the curve does not depend on one lucky sample), +# and compare the result with the truth. + +WEIGHTS = (0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30) +NOISE = 0.1 +DRAWS = 5 + +curves = {"exact data": [], "noisy data": []} +misfits = {"exact data": [], "noisy data": []} +not_solved = {"exact data": [], "noisy data": []} +for weight in WEIGHTS: + interpolator, solved = solve(exact_values, exact_normals, weight) + error, misfit = errors(interpolator, exact_values) + curves["exact data"].append(error) + misfits["exact data"].append(misfit) + not_solved["exact data"].append(not solved) + + draw_errors, draw_misfits, draw_unsolved = [], [], 0 + for draw in range(DRAWS): + values, normals = add_noise(NOISE, seed=100 + draw) + interpolator, solved = solve(values, normals, weight) + error, misfit = errors(interpolator, values) + draw_errors.append(error) + draw_misfits.append(misfit) + draw_unsolved += not solved + curves["noisy data"].append(np.mean(draw_errors)) + misfits["noisy data"].append(np.mean(draw_misfits)) + not_solved["noisy data"].append(draw_unsolved > 0) + +print(f"{'weight':>8} | {'exact: error':>12} {'misfit':>8} | {'noisy: error':>12} {'misfit':>8}") +for i, weight in enumerate(WEIGHTS): + flag = " <- solver did not converge" if not_solved["exact data"][i] or not_solved["noisy data"][i] else "" + print( + f"{weight:>8} | {curves['exact data'][i]:12.4f} {misfits['exact data'][i]:8.4f} |" + f" {curves['noisy data'][i]:12.4f} {misfits['noisy data'][i]:8.4f}{flag}" + ) + +# %% +# The error against the truth, and the sections +# --------------------------------------------- +# Left: error against the true field for each weight. The default weight (1) is marked. Right: +# the middle section through the dome for one noisy data set at a weak, a moderate and a strong +# weight, against the true horizons (dashed). + +fig = plt.figure(figsize=(13, 4.2)) +ax = fig.add_subplot(1, 4, 1) +for label, colour in (("exact data", "tab:blue"), ("noisy data", "tab:red")): + ax.plot(WEIGHTS, curves[label], "o-", color=colour, label=label) +ax.axvline(1, color="grey", linestyle=":") +ax.set_xscale("log") +ax.set_xlabel("regularisation weight") +ax.set_ylabel("rms error against the truth") +ax.legend() + +noisy_values, noisy_normals = add_noise(NOISE, seed=100) +x = np.linspace(0, 1, 100) +xx, zz = np.meshgrid(x, x) +section = np.c_[xx.ravel(), np.full(xx.size, 0.5), zz.ravel()] +for i, weight in enumerate((0.03, 1, 30), start=2): + interpolator, _ = solve(noisy_values, noisy_normals, weight) + sub = fig.add_subplot(1, 4, i) + sub.contour(xx, zz, truth(section).reshape(xx.shape), levels=HORIZONS, colors="k", linestyles="--") + sub.contour( + xx, zz, interpolator.evaluate_value(section).reshape(xx.shape), levels=HORIZONS, colors="tab:red" + ) + near = np.abs(noisy_values[:, 1] - 0.5) < 0.2 + sub.plot(noisy_values[near, 0], noisy_values[near, 2], "ko", ms=3) + sub.set_title(f"weight {weight}") + sub.set_aspect("equal") +fig.tight_layout() +plt.show() + +# %% +# Reading the results +# ------------------- +# * With exact data the error rises with the weight, because any smoothing moves the surface +# away from the truth, which the data describe. The weakest weights fit the data almost +# exactly (misfit near zero). +# * With noisy data the picture changes. A very small weight fits the noise: the error is +# several times larger than at the default, and the horizons in the left section wiggle. A +# very large weight flattens the dome: the error levels off near 0.12, which is roughly the +# error of a flat horizon. The lowest error is in between, near the default weight of 1. +# * The data misfit always grows with the weight. It tells you how far the surface has been +# pulled from the data, but a small misfit is not the same as a small error: the noisy data +# are matched best by the weights that are furthest from the truth. +# +# The right weight depends on how much you trust the data, which is a geological judgement. +# A practical approach, as in the debugging guide, is to lower the weight until the surface +# starts to follow noise, then raise it until the structure starts to flatten, and pick a value +# between the two. The data weights (``cpw`` for contacts, ``npw`` for orientation vectors, +# ``gpw`` and ``tpw`` for the other constraint types) are the same trade-off from the other side. diff --git a/examples/5_interpolation/plot_3_fold_constraints.py b/examples/5_interpolation/plot_3_fold_constraints.py new file mode 100644 index 000000000..cff178fa5 --- /dev/null +++ b/examples/5_interpolation/plot_3_fold_constraints.py @@ -0,0 +1,236 @@ +""" +5c. Fold constraints with sparse data +===================================== + +A fold can be very hard to model from a few measurements: two contact points and one +orientation say little about a surface that folds back on itself. A plain interpolator, +which only asks for a smooth surface, will return something close to a plane. + +The fold interpolators add the knowledge that the surface is folded. Given a *fold frame* (a +coordinate system with one axis across the fold and one along it) and a description of how the +limb dip changes across the fold, they constrain the *direction* of the gradient everywhere, +not only where there are data. + +This example measures what that adds. The truth is a cylindrical fold, and each of six +interpolators is given the same three data: plain ``FDI``, ``PLI`` and ``P2``, and their +fold-constrained versions ``FDFold``, ``DFI`` (the fold version of PLI) and ``P2Fold``. + +Real models get the fold frame from ``GeologicalModel.create_and_add_fold_frame``; see the +fold examples. Here we write it down analytically so that the answer is exactly known. +""" + +import time + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from loop_interpolation import ( + DiscreteFoldInterpolator, + FDFoldInterpolator, + FiniteDifferenceInterpolator, + FoldEvent, + P1Interpolator, + P2FoldInterpolator, + P2Interpolator, + P2UnstructuredTetMesh, + StructuredGrid, + TetMesh, +) +from scipy.optimize import brentq + +# %% +# The folded surface +# ------------------ +# A cylindrical fold with its axis along ``y``: the horizons are ``s = z - h(x)`` with +# ``h(x) = A sin(2 pi x / wavelength)``, in the unit cube. + +AMPLITUDE, WAVELENGTH = 0.12, 0.6 +K = 2 * np.pi / WAVELENGTH + + +def height(x): + return AMPLITUDE * np.sin(K * x) + + +def height_slope(x): + return AMPLITUDE * K * np.cos(K * x) + + +def truth(xyz): + return xyz[:, 2] - height(xyz[:, 0]) + + +def truth_gradient(xyz): + return np.c_[-height_slope(xyz[:, 0]), np.zeros(len(xyz)), np.ones(len(xyz))] + + +# %% +# Three data: two contact points on the 0.5 horizon, far apart, and one orientation in the middle. + + +def contact(x, y, level=0.5): + z = brentq(lambda z: truth(np.array([[x, y, z]]))[0] - level, 0.0, 1.0) + return [x, y, z, level] + + +values = np.array([contact(0.2, 0.3), contact(0.8, 0.7)]) +orientation_xyz = np.array([[0.5, 0.5, 0.5]]) +normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] + +# %% +# The fold frame and the fold event +# --------------------------------- +# The fold frame's first coordinate ``gx`` runs across the fold and the second, ``gy``, along +# the axis. For this fold they are just ``x`` and ``y``. A ``FoldEvent`` combines the frame with +# the fold axis and a function giving the limb rotation angle (the dip of the limb) in terms of +# ``gx``. For a cylindrical fold the axis is constant. The signs follow the library's +# convention for the direction of the frame's third axis. + + +class Coordinate: + """Scalar field equal to one coordinate; its gradient is the unit vector.""" + + def __init__(self, axis): + self.axis = axis + + def evaluate_value(self, xyz): + return np.asarray(xyz)[:, self.axis].astype(float) + + def evaluate_gradient(self, xyz): + gradient = np.zeros((len(xyz), 3)) + gradient[:, self.axis] = 1.0 + return gradient + + +class FoldFrame: + features = (Coordinate(0), Coordinate(1)) + + +def limb_rotation(gx): + return -np.rad2deg(np.arctan(height_slope(gx))) + + +def make_fold_event(): + return FoldEvent( + FoldFrame(), fold_limb_rotation=limb_rotation, fold_axis=np.array([0.0, -1.0, 0.0]) + ) + + +# %% +# Six interpolators +# ----------------- +# FDI/FDFold use a grid of 16 cells per side, PLI/DFI a tetrahedral mesh of the same grid, and +# P2/P2Fold a mesh of quadratic tetrahedra with 8 cells per side. + +CELLS = 16 + + +def grid(): + return StructuredGrid( + origin=np.zeros(3), nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + + +def tetrahedra(): + return TetMesh( + origin=np.zeros(3), nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + + +def quadratic(): + return P2UnstructuredTetMesh( + origin=np.zeros(3), nsteps_cells=np.full(3, 8), step_vector=np.full(3, 1 / 8) + ) + + +FACTORIES = { + "FDI": lambda: FiniteDifferenceInterpolator(grid()), + "PLI": lambda: P1Interpolator(tetrahedra()), + "P2": lambda: P2Interpolator(quadratic()), + "FDFold": lambda: FDFoldInterpolator(grid(), fold=make_fold_event()), + "DFI": lambda: DiscreteFoldInterpolator(tetrahedra(), fold=make_fold_event()), + "P2Fold": lambda: P2FoldInterpolator(quadratic(), fold=make_fold_event()), +} + + +def interpolate(kind): + interpolator = FACTORIES[kind]() + interpolator.set_value_constraints(values) + interpolator.set_normal_constraints(normals) + start = time.perf_counter() + interpolator.setup_interpolator() + interpolator.solve_system("cg") + return interpolator, time.perf_counter() - start + + +results = {kind: interpolate(kind) for kind in FACTORIES} + +# %% +# Measuring the result +# -------------------- +# On a regular set of points inside the cube compare each interpolator with the truth: +# +# * ``rms error``, the error in the scalar field (the position of the horizons), +# * ``median gradient angle``, how far the gradient direction (the dip of the surface) is from +# the true one, in degrees. + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T + + +def angle_degrees(a, b): + cosine = np.sum(a * b, axis=1) / np.linalg.norm(a, axis=1) / np.linalg.norm(b, axis=1) + return np.rad2deg(np.arccos(np.clip(cosine, -1, 1))) + + +rows = {} +for kind, (interpolator, elapsed) in results.items(): + error = interpolator.evaluate_value(points) - truth(points) + angle = angle_degrees(interpolator.evaluate_gradient(points), truth_gradient(points)) + rows[kind] = { + "unknowns": interpolator.dof, + "rms error": np.sqrt(np.mean(error**2)), + "median gradient angle": np.median(angle), + "time (s)": elapsed, + } +print(pd.DataFrame(rows).T.astype({"unknowns": int}).round(3).to_string()) + +# %% +# The section through the middle of the fold +# ------------------------------------------ +# Dashed lines are the true horizons, red lines are each interpolator's. Contact points near +# the section are black dots. + +x = np.linspace(0, 1, 100) +xx, zz = np.meshgrid(x, x) +section = np.c_[xx.ravel(), np.full(xx.size, 0.5), zz.ravel()] +true_field = truth(section).reshape(xx.shape) + +fig, axes = plt.subplots(2, 3, figsize=(12, 7), sharex=True, sharey=True) +for ax, (kind, (interpolator, _)) in zip(axes.ravel(), results.items()): + field = interpolator.evaluate_value(section).reshape(xx.shape) + ax.contour(xx, zz, true_field, levels=[0.3, 0.5, 0.7], colors="k", linestyles="--", linewidths=1) + ax.contour(xx, zz, field, levels=[0.3, 0.5, 0.7], colors="tab:red", linewidths=1.5) + ax.plot(values[:, 0], values[:, 2], "ko", ms=5) + ax.set_title(f"{kind}: gradient angle {rows[kind]['median gradient angle']:.0f} deg") + ax.set_aspect("equal") +fig.tight_layout() +plt.show() + +# %% +# What to take from this +# ---------------------- +# * The three plain interpolators return straight, gently tilted horizons, about 45 degrees +# out in the direction of the surface. With three data nothing in the problem says the surface +# is folded. +# * All three fold interpolators bring the direction to within about 5 degrees. This is what +# the fold constraint adds: the dip of the limbs everywhere, not only near the data. +# * Getting the direction right does not by itself put the horizons in the right place. +# ``FDFold`` and ``DFI`` recover the folded shape but exaggerate its amplitude, and their +# scalar-field error (0.11 and 0.12) is only modestly better than the plain interpolators +# (about 0.16). ``P2Fold`` is the only one close to the truth in both respects. This example +# does not investigate why; it is one set of data on one fold. +# +# The fold event here is exact because we wrote the geometry down. In a real model it is +# fitted to the structural data, so the quality of the fold constraint is only as good as that +# fit: see the fold examples, and the S-plot checks in the debugging guide. diff --git a/examples/README.rst b/examples/README.rst index c4ece1f82..46a57ac34 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -12,6 +12,9 @@ learn them: and updating fault geometry after a model has been built. 4. **Advanced use** - building models directly from geological map data, logging, controlling data weighting, and comparing interpolators. +5. **Interpolation** - how the interpolators behave, measured against synthetic + models with a known answer: choosing an interpolator, the regularisation weight, and + fold constraints with sparse data. Each example is a standalone, runnable Python script. Most examples load one of the sample datasets bundled in :code:`LoopStructural.datasets`, so From fc1da01cc97bc3dc548ab3da9a23d1aa3371f727 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 14:34:08 +1030 Subject: [PATCH 04/10] docs: add refolded-fold and orientation-data examples to the Interpolation section - 5d refolded folds: a dome-and-basin interference pattern. A fold event with a rotating axis recovers it; a constant axis leaves the median direction error at 18-19 degrees and the horizon collapses to stripes, even though the scalar-field error looks as good - 5e kinds of orientation data with InterpolatorBuilder: with six contacts, error falls from contacts only through one tangent, two tangents or a gradient direction, to a full normal Co-Authored-By: Claude Sonnet 5.5 --- .../5_interpolation/plot_4_refolded_folds.py | 273 ++++++++++++++++++ .../plot_5_kinds_of_orientation_data.py | 161 +++++++++++ examples/README.rst | 2 +- 3 files changed, 435 insertions(+), 1 deletion(-) create mode 100644 examples/5_interpolation/plot_4_refolded_folds.py create mode 100644 examples/5_interpolation/plot_5_kinds_of_orientation_data.py diff --git a/examples/5_interpolation/plot_4_refolded_folds.py b/examples/5_interpolation/plot_4_refolded_folds.py new file mode 100644 index 000000000..4cec24b41 --- /dev/null +++ b/examples/5_interpolation/plot_4_refolded_folds.py @@ -0,0 +1,273 @@ +""" +5d. Refolded folds: when the fold axis changes +============================================== + +In a refolded fold the axis of the first fold is itself folded by a later event, so the fold +axis is not constant: it plunges one way here and the other way a little further along. A fold +constraint that assumes a straight axis then pushes the surface the wrong way wherever the axis +has rotated. + +This example builds the simplest refolded geometry, a dome-and-basin interference pattern, and +asks what the fold interpolators need to recover it from three contacts and six orientations: + +* the early fold has its axis along ``y`` and shape ``sin(x)``, +* a later fold with axis along ``x`` and shape ``sin(y)`` tilts that axis, so its plunge + changes along ``y``. + +The scalar field is ``s = z - A2 sin(k2 x) - B3 sin(k3 y)``. The fold event describes this with +two functions: the limb rotation as a function of the coordinate across the fold, and the +*axis rotation* as a function of the coordinate along it. We compare an event that includes the +axis rotation with one that assumes a constant axis. For the sparse-data fold with a constant +axis see the previous example. +""" + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from loop_interpolation import ( + DiscreteFoldInterpolator, + FDFoldInterpolator, + FiniteDifferenceInterpolator, + FoldEvent, + P1Interpolator, + P2FoldInterpolator, + P2Interpolator, + P2UnstructuredTetMesh, + StructuredGrid, + TetMesh, +) +from scipy.optimize import brentq + +# %% +# The refolded surface +# -------------------- + +A2, WAVELENGTH2 = 0.10, 0.5 # early fold: amplitude and wavelength (axis along y) +B3, WAVELENGTH3 = 0.08, 0.8 # later fold: amplitude and wavelength (axis along x) +K2, K3 = 2 * np.pi / WAVELENGTH2, 2 * np.pi / WAVELENGTH3 + + +def truth(xyz): + return xyz[:, 2] - A2 * np.sin(K2 * xyz[:, 0]) - B3 * np.sin(K3 * xyz[:, 1]) + + +def truth_gradient(xyz): + return np.c_[ + -A2 * K2 * np.cos(K2 * xyz[:, 0]), + -B3 * K3 * np.cos(K3 * xyz[:, 1]), + np.ones(len(xyz)), + ] + + +# %% +# Data: three contact points on the 0.5 horizon and six orientations at random places. + + +def contact(x, y, level=0.5): + z = brentq(lambda z: truth(np.array([[x, y, z]]))[0] - level, 0.0, 1.0) + return [x, y, z, level] + + +values = np.array([contact(0.25, 0.25), contact(0.75, 0.5), contact(0.4, 0.8)]) +orientation_xyz = np.random.default_rng(5).uniform(0.15, 0.85, (6, 3)) +normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] + +# %% +# The fold frame and two fold events +# ---------------------------------- +# As in the previous example the fold frame is just ``gx = x`` (across the early fold) and +# ``gy = y`` (along it). The limb rotation is the dip of the early fold's limbs. +# +# The axis rotation is the plunge of the early fold's axis, which the later fold varies with +# ``gy``. The ``+ 180`` is a sign convention of the library: the frame's third axis must point +# against the orientation vectors. + + +class Coordinate: + def __init__(self, axis): + self.axis = axis + + def evaluate_value(self, xyz): + return np.asarray(xyz)[:, self.axis].astype(float) + + def evaluate_gradient(self, xyz): + gradient = np.zeros((len(xyz), 3)) + gradient[:, self.axis] = 1.0 + return gradient + + +class FoldFrame: + features = (Coordinate(0), Coordinate(1)) + + +def limb_rotation(gx): + return -np.rad2deg(np.arctan(A2 * K2 * np.cos(K2 * gx))) + + +def axis_rotation(gy): + return np.rad2deg(np.arctan(B3 * K3 * np.cos(K3 * gy))) + 180.0 + + +def fold_with_rotating_axis(): + return FoldEvent(FoldFrame(), fold_limb_rotation=limb_rotation, fold_axis_rotation=axis_rotation) + + +def fold_with_constant_axis(): + return FoldEvent( + FoldFrame(), fold_limb_rotation=limb_rotation, fold_axis=np.array([0.0, -1.0, 0.0]) + ) + + +# %% +# Interpolators +# ------------- +# The three plain interpolators, the three fold interpolators with the rotating axis, and the +# same two of them (``DFI`` and ``P2Fold``) with a constant axis. + +CELLS = 16 + + +def grid(): + return StructuredGrid( + origin=np.zeros(3), nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + + +def tetrahedra(): + return TetMesh( + origin=np.zeros(3), nsteps_cells=np.full(3, CELLS), step_vector=np.full(3, 1 / CELLS) + ) + + +def quadratic(): + return P2UnstructuredTetMesh( + origin=np.zeros(3), nsteps_cells=np.full(3, 8), step_vector=np.full(3, 1 / 8) + ) + + +FACTORIES = { + "FDI": lambda: FiniteDifferenceInterpolator(grid()), + "PLI": lambda: P1Interpolator(tetrahedra()), + "P2": lambda: P2Interpolator(quadratic()), + "FDFold": lambda: FDFoldInterpolator(grid(), fold=fold_with_rotating_axis()), + "DFI": lambda: DiscreteFoldInterpolator(tetrahedra(), fold=fold_with_rotating_axis()), + "P2Fold": lambda: P2FoldInterpolator(quadratic(), fold=fold_with_rotating_axis()), + "DFI (constant axis)": lambda: DiscreteFoldInterpolator( + tetrahedra(), fold=fold_with_constant_axis() + ), + "P2Fold (constant axis)": lambda: P2FoldInterpolator( + quadratic(), fold=fold_with_constant_axis() + ), +} + + +def interpolate(kind): + interpolator = FACTORIES[kind]() + interpolator.set_value_constraints(values) + interpolator.set_normal_constraints(normals) + interpolator.setup_interpolator() + interpolator.solve_system("cg") + return interpolator + + +interpolators = {kind: interpolate(kind) for kind in FACTORIES} + +# %% +# Results +# ------- +# * ``rms error``: error in the scalar field over the interior, +# * ``median`` and ``90th percentile gradient angle``: how far the dip direction is from the +# truth (degrees), +# * ``worst angle at the data``: the largest angle between the solution's gradient and one of +# the six orientation measurements. A solution that respects the data should be near zero. + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T + + +def angle_degrees(a, b): + cosine = np.sum(a * b, axis=1) / np.linalg.norm(a, axis=1) / np.linalg.norm(b, axis=1) + return np.rad2deg(np.arccos(np.clip(cosine, -1, 1))) + + +rows = {} +for kind, interpolator in interpolators.items(): + error = interpolator.evaluate_value(points) - truth(points) + angle = angle_degrees(interpolator.evaluate_gradient(points), truth_gradient(points)) + at_data = angle_degrees(interpolator.evaluate_gradient(normals[:, :3]), normals[:, 3:]) + rows[kind] = { + "rms error": np.sqrt(np.mean(error**2)), + "median angle": np.median(angle), + "p90 angle": np.percentile(angle, 90), + "worst angle at the data": np.max(at_data), + } +print(pd.DataFrame(rows).T.round(3).to_string()) + +# %% +# The height of the 0.5 horizon +# ----------------------------- +# A map of the height ``z`` of the horizon ``s = 0.5`` for the truth, a plain interpolator, and +# the fold interpolator with a constant and with a rotating axis. + +x = np.linspace(0.05, 0.95, 40) +xx, yy = np.meshgrid(x, x) +z_levels = np.linspace(0.0, 1.0, 80) + + +def horizon_height(field): + """Height z at which ``field(x, y, z) = 0.5``, found by linear interpolation along z.""" + heights = np.full(xx.shape, np.nan) + lattice = np.array(np.meshgrid(x, x, z_levels, indexing="ij")).reshape(3, -1).T + values_ = np.asarray(field(lattice)).reshape(len(x), len(x), len(z_levels)) - 0.5 + for i in range(len(x)): + for j in range(len(x)): + column = values_[i, j] + crossing = np.where(np.sign(column[:-1]) * np.sign(column[1:]) < 0)[0] + if len(crossing): + k = crossing[0] + heights[j, i] = z_levels[k] + (z_levels[k + 1] - z_levels[k]) * column[k] / ( + column[k] - column[k + 1] + ) + return heights + + +panels = { + "truth": truth, + "FDI (no fold constraint)": interpolators["FDI"].evaluate_value, + "P2Fold, constant axis": interpolators["P2Fold (constant axis)"].evaluate_value, + "P2Fold, rotating axis": interpolators["P2Fold"].evaluate_value, +} +fig, axes = plt.subplots(1, 4, figsize=(15, 3.8), sharey=True) +for ax, (title, field) in zip(axes, panels.items()): + image = ax.contourf(xx, yy, horizon_height(field), levels=np.linspace(0.35, 0.65, 13), extend="both") + ax.plot(values[:, 0], values[:, 1], "ko", ms=4) + ax.set_title(title) + ax.set_xlabel("x") + ax.set_aspect("equal") +axes[0].set_ylabel("y") +fig.colorbar(image, ax=axes, label="height of the 0.5 horizon", shrink=0.8) +plt.show() + +# %% +# What to take from this +# ---------------------- +# * The plain interpolators miss the dome-and-basin pattern: the gradient direction is about 40 +# degrees out and the horizon is a smooth ramp. +# * With the axis rotation included, all three fold interpolators get the median direction to +# within 2 to 10 degrees, and ``P2Fold`` reproduces the pattern in the map. +# * With a constant axis the median direction error is about 18 to 19 degrees for both ``DFI`` +# and ``P2Fold``, and the map shows what has gone wrong: the horizon is straight stripes +# along ``y``. The later fold, which varies the axis, has been lost entirely. +# * A single error number can hide this. ``P2Fold`` with a constant axis has a scalar-field +# error (0.065) as small as with the rotating axis (0.062). Look at the gradient direction +# and at maps or sections, not only at the error. +# * The fold constraint can override the data. ``FDFold`` ends 40 degrees from one of the six +# orientation measurements and ``DFI`` 14 degrees, while ``P2Fold`` stays within 6 degrees. +# This example does not examine why. When a fold constraint is used it is worth checking the +# solution against the measurements. +# * A fold constraint improves the direction without necessarily improving the scalar-field +# error: ``DFI`` has a larger error (0.16) than the plain interpolators (about 0.10). +# +# The fold event can describe the limb rotation only as a function of the coordinate across the +# fold, so for this surface the constraint direction is itself a few degrees off in places. +# Treat these numbers as one example, not a benchmark. diff --git a/examples/5_interpolation/plot_5_kinds_of_orientation_data.py b/examples/5_interpolation/plot_5_kinds_of_orientation_data.py new file mode 100644 index 000000000..43b134041 --- /dev/null +++ b/examples/5_interpolation/plot_5_kinds_of_orientation_data.py @@ -0,0 +1,161 @@ +""" +5e. Kinds of orientation data, with InterpolatorBuilder +======================================================= + +Contact points tell the interpolator where a surface is. Orientation data tell it which way +the surface faces, but there are several ways to say that, and they carry different amounts of +information: + +* a **tangent** is one direction that lies *in* the surface (for example a strike line or a + lineation). One tangent leaves the surface free to rotate about it, +* **two tangents** that are not parallel fix the plane, and so the direction of the gradient, +* a **gradient** constraint gives the direction of the gradient (the surface normal) but not its + length, +* a **normal** constraint gives the full gradient vector, direction and length. Its length + says how quickly the scalar field changes across the surface, so it also fixes the scale of + the field. + +This example uses ``InterpolatorBuilder``, which assembles an interpolator from a bounding box +and arrays of constraints, and measures what each kind of orientation data adds when there are +few contact points. The truth is the same gaussian dome as in the earlier examples. +""" + +import numpy as np +import pandas as pd +from scipy.optimize import brentq + +from LoopStructural import BoundingBox, InterpolatorBuilder + +# %% +# The model and its data +# ---------------------- + +RELIEF, CENTRE, WIDTH = 0.2, 0.5, 0.2 +HORIZONS = (0.3, 0.5, 0.7) + + +def truth(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + return xyz[:, 2] - RELIEF * np.exp(-r2 / (2 * WIDTH**2)) + + +def truth_gradient(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + e = RELIEF * np.exp(-r2 / (2 * WIDTH**2)) / WIDTH**2 + return np.c_[e * (xyz[:, 0] - CENTRE), e * (xyz[:, 1] - CENTRE), np.ones(len(xyz))] + + +def borehole_contacts(xy): + points = [] + for x, y in xy: + for level in HORIZONS: + + def misfit(z, x=x, y=y, level=level): + return truth(np.array([[x, y, z]]))[0] - level + + if misfit(0.0) * misfit(1.0) < 0: + points.append([x, y, brentq(misfit, 0.0, 1.0), level]) + return np.array(points) + + +# %% +# Few contacts: two boreholes, so six contact points. Five orientation measurements, each +# written four ways. + +contacts = borehole_contacts(np.array([[0.5, 0.5], [0.35, 0.5]])) +print(f"{len(contacts)} contact points") + +orientation_xyz = np.array( + [[0.4, 0.4, 0.5], [0.6, 0.6, 0.5], [0.3, 0.7, 0.6], [0.8, 0.3, 0.4], [0.5, 0.5, 0.5]] +) +gradient = truth_gradient(orientation_xyz) +unit_normal = gradient / np.linalg.norm(gradient, axis=1, keepdims=True) + +# two directions in the surface at each point: ``first`` is horizontal where possible +reference = np.where( + np.abs(unit_normal[:, [2]]) > 0.9, np.array([[1.0, 0, 0]]), np.array([[0, 0, 1.0]]) +) +first = np.cross(unit_normal, reference) +first /= np.linalg.norm(first, axis=1, keepdims=True) +second = np.cross(unit_normal, first) + +constraints = { + "contacts only": {}, + "one tangent per point": {"tangent": np.c_[orientation_xyz, first]}, + "two tangents per point": { + "tangent": np.vstack([np.c_[orientation_xyz, first], np.c_[orientation_xyz, second]]) + }, + "gradient direction": {"gradient": np.c_[orientation_xyz, unit_normal]}, + "normal (direction and length)": {"normal": np.c_[orientation_xyz, gradient]}, +} + +# %% +# Building each interpolator +# -------------------------- +# ``InterpolatorBuilder`` takes the interpolator type, the bounding box of the model and the +# approximate number of elements. Each ``add_*_constraints`` call returns the builder, so the +# calls chain, and ``build()`` returns the interpolator. The arrays are: +# +# =============== ========================================== +# method columns +# =============== ========================================== +# value x, y, z, value +# normal x, y, z, nx, ny, nz +# gradient x, y, z, gx, gy, gz (direction only) +# tangent x, y, z, tx, ty, tz (lies in the surface) +# =============== ========================================== + +bounding_box = BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[1.0, 1.0, 1.0]) + + +def build(extra): + builder = InterpolatorBuilder( + "FDI", bounding_box=bounding_box, nelements=4000, buffer=0.1 + ).add_value_constraints(contacts) + for kind, array in extra.items(): + getattr(builder, f"add_{kind}_constraints")(array) + return builder.build() + + +interpolators = {label: build(extra) for label, extra in constraints.items()} + +# %% +# What each kind of data adds +# --------------------------- +# Compare with the truth on a regular set of interior points. + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T + + +def angle_degrees(a, b): + cosine = np.sum(a * b, axis=1) / np.linalg.norm(a, axis=1) / np.linalg.norm(b, axis=1) + return np.rad2deg(np.arccos(np.clip(cosine, -1, 1))) + + +rows = {} +for label, interpolator in interpolators.items(): + error = interpolator.evaluate_value(points) - truth(points) + angle = angle_degrees(interpolator.evaluate_gradient(points), truth_gradient(points)) + rows[label] = { + "rms error": np.sqrt(np.nanmean(error**2)), + "median gradient angle": np.nanmedian(angle), + } +print(pd.DataFrame(rows).T.round(3).to_string()) + +# %% +# What to take from this +# ---------------------- +# * Each step adds information and the error falls in order: contacts only, one tangent, two +# tangents or a gradient direction, then a full normal. +# * One tangent helps only a little. It removes one of the two directions the surface could +# rotate in, and leaves the other free. +# * Two non-parallel tangents are nearly equivalent to a gradient direction: both fix the +# plane, and the errors are within a few percent of each other. +# * A normal is the most informative because its length constrains the scale of the field as +# well as its direction. That is also why it is the usual choice in ``GeologicalModel`` +# (the ``nx, ny, nz`` columns); use a gradient direction when you have no reliable measure of +# how fast the field changes, and tangents when the data really are lines in the surface. +# +# These are errors from one set of six contacts and five orientations. With more contact points +# the differences shrink, because the contacts alone already constrain the surface well. diff --git a/examples/README.rst b/examples/README.rst index 46a57ac34..69348d2f5 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -14,7 +14,7 @@ learn them: logging, controlling data weighting, and comparing interpolators. 5. **Interpolation** - how the interpolators behave, measured against synthetic models with a known answer: choosing an interpolator, the regularisation weight, and - fold constraints with sparse data. + fold constraints with sparse data, refolded folds, and the kinds of orientation data. Each example is a standalone, runnable Python script. Most examples load one of the sample datasets bundled in :code:`LoopStructural.datasets`, so From d6685b0649cb25b0345f9dda3876ba163813e915 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 17:21:26 +1030 Subject: [PATCH 05/10] docs: add a native solver example to the Interpolation section 5f solves the dome with about 55 000 unknowns with the Python CG solver and the native solver: the fields agree to the solver tolerance (9e-9) and the native solve is faster (2.2 s against 3.9 s), the Python-side peak memory is about 1 MB against 102 MB when the solver is chosen before setup (28 MB when chosen after), and warm_start cuts the iterations of a re-solve by a third. It also shows the thing that catches people out: the native solver stops at its default maxiter=5000 without converging, logs a warning, leaves up_to_date False, and returns a field that is 0.9 off, although the residual in the warning is close to the tolerance. The user guide gets a short section on checking convergence. The example imports the compiled core, which has wheels for Linux, macOS and Windows on PyPI, so the docs requirements and both docs workflows install it. Co-Authored-By: Claude Sonnet 5.5 --- .github/workflows/docs-check.yml | 2 +- .github/workflows/documentation.yml | 1 + docs/requirements.txt | 2 + docs/source/user_guide/native_solver.rst | 9 + .../5_interpolation/plot_6_native_solver.py | 229 ++++++++++++++++++ examples/README.rst | 3 +- 6 files changed, 244 insertions(+), 2 deletions(-) create mode 100644 examples/5_interpolation/plot_6_native_solver.py diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml index 4f4a23e67..f05a03c86 100644 --- a/.github/workflows/docs-check.yml +++ b/.github/workflows/docs-check.yml @@ -42,7 +42,7 @@ jobs: - name: Install the packages and the docs extra run: | - pip install -e packages/loop_common -e packages/loop_interpolation \ + pip install -e packages/loop_common -e "packages/loop_interpolation[native]" \ -e "packages/loopstructural_visualisation[all]" pip install -e ".[docs]" diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 5763251b6..b212008a8 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -42,6 +42,7 @@ jobs: - name: Install PyVista and dependencies run: | pip install -e .[docs] + pip install loop-interpolation-core - name: Install custom OSMesa VTK variant run: | diff --git a/docs/requirements.txt b/docs/requirements.txt index 0a9c23e62..92c8f02d5 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,5 @@ # Build the docs with the project and its docs extra: # pip install -r docs/requirements.txt .[docs] +# the native solver example needs the compiled core (wheels for Linux, macOS and Windows) +loop-interpolation-core diff --git a/docs/source/user_guide/native_solver.rst b/docs/source/user_guide/native_solver.rst index b6cbe752d..7abf9a115 100644 --- a/docs/source/user_guide/native_solver.rst +++ b/docs/source/user_guide/native_solver.rst @@ -64,6 +64,15 @@ A start value needs 8 more bytes for each node until the solve starts. The memor Use ``warm_start`` when you solve again after a small change, for example a new weight. The solve then needs fewer iterations. Any other keyword is rejected with a ``NotImplementedError``. The message lists the accepted keywords. +Checking that the solve converged +--------------------------------- +The solver stops after ``maxiter`` iterations even when it has not reached ``tol``. When this happens it logs a warning +(``The native solver did not converge``), sets ``interpolator.up_to_date`` to ``False``, and leaves the unfinished field in the interpolator. +Nothing is raised. On a test problem with about 55 000 unknowns the default limit of 5000 iterations was not enough: the field was far +from the Python solution, although the residual in the warning was close to ``tol``. Raising ``maxiter`` fixed it. +Check ``up_to_date`` after a native solve, and raise ``maxiter`` if it is ``False``. +The gallery example :doc:`/_auto_examples/5_interpolation/plot_6_native_solver` shows this, and the memory saving, on a worked problem. + Unsupported settings -------------------- The native solver does not support these settings. When one of them is set, the solver raises ``NotImplementedError``. The message contains the name of the setting. diff --git a/examples/5_interpolation/plot_6_native_solver.py b/examples/5_interpolation/plot_6_native_solver.py new file mode 100644 index 000000000..a0ecbe01f --- /dev/null +++ b/examples/5_interpolation/plot_6_native_solver.py @@ -0,0 +1,229 @@ +""" +5f. The native solver +===================== + +The native solver is a Rust implementation of the same conjugate gradient solve that the Python +interpolator does. It solves the same discrete system, uses several threads, and, when it is chosen +*before* the interpolator is set up, the Python side never builds the sparse matrix, which is where +most of the memory goes. + +This example checks the claims on a problem big enough to matter (about 55 000 unknowns) and looks +at the one thing that catches people out: the native solver stops at a default iteration limit, and +what it returns when it stops early is not the answer. The native solver needs the compiled +package ``loop_interpolation_core``; install it with ``pip install loop-interpolation[native]``. +It works on uniform 3D grids, which is what ``FDI`` uses (see the user guide for the full list of +limits). +""" + +import importlib.util +import time +import tracemalloc + +import numpy as np +import pandas as pd +from loop_common.geometry import BoundingBox +from scipy.optimize import brentq + +from LoopStructural import InterpolatorBuilder + +if importlib.util.find_spec("loop_interpolation_core") is None: + raise ImportError( + "This example needs the compiled native core: pip install loop-interpolation[native]" + ) + +# %% +# The problem +# ----------- +# The gaussian dome used in the earlier interpolation examples, in the unit cube, sampled with nine +# boreholes and five orientation measurements. + +RELIEF, CENTRE, WIDTH = 0.2, 0.5, 0.2 +HORIZONS = (0.3, 0.5, 0.7) + + +def truth(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + return xyz[:, 2] - RELIEF * np.exp(-r2 / (2 * WIDTH**2)) + + +def truth_gradient(xyz): + r2 = (xyz[:, 0] - CENTRE) ** 2 + (xyz[:, 1] - CENTRE) ** 2 + e = RELIEF * np.exp(-r2 / (2 * WIDTH**2)) / WIDTH**2 + return np.c_[e * (xyz[:, 0] - CENTRE), e * (xyz[:, 1] - CENTRE), np.ones(len(xyz))] + + +def borehole_contacts(xy): + points = [] + for x, y in xy: + for level in HORIZONS: + + def misfit(z, x=x, y=y, level=level): + return truth(np.array([[x, y, z]]))[0] - level + + if misfit(0.0) * misfit(1.0) < 0: + points.append([x, y, brentq(misfit, 0.0, 1.0), level]) + return np.array(points) + + +rng = np.random.default_rng(4) +borehole_xy = np.vstack( + [[0.5, 0.5], [0.35, 0.5], [0.65, 0.5], [0.5, 0.35], [0.5, 0.65], rng.uniform(0.1, 0.9, (4, 2))] +) +values = borehole_contacts(borehole_xy) +orientation_xyz = np.array( + [[0.4, 0.4, 0.5], [0.6, 0.6, 0.5], [0.3, 0.7, 0.6], [0.8, 0.3, 0.4], [0.5, 0.5, 0.5]] +) +normals = np.c_[orientation_xyz, truth_gradient(orientation_xyz)] + +axis = np.linspace(0.1, 0.9, 9) +points = np.array(np.meshgrid(axis, axis, axis, indexing="ij")).reshape(3, -1).T +bounding_box = BoundingBox(origin=[0.0, 0.0, 0.0], maximum=[1.0, 1.0, 1.0]) + +# %% +# A helper that builds, solves and measures +# ----------------------------------------- +# ``InterpolatorBuilder`` chooses the solver with ``use_solver``. Whether that happens before or +# after ``setup_interpolator`` matters for the native solver, so the helper takes it as a +# parameter. ``tracemalloc`` measures the *Python* memory used while setting up and solving; the +# native core allocates its own memory outside of Python, which this does not see. + + +def run(solver, *, choose_solver="before", **solver_kwargs): + builder = InterpolatorBuilder( + "FDI", bounding_box=bounding_box, nelements=50000, buffer=0.1 + ) + builder.add_value_constraints(values).add_normal_constraints(normals) + if choose_solver == "before": + builder.use_solver(solver, **solver_kwargs) + + tracemalloc.start() + start = time.perf_counter() + builder.setup_interpolator() + if choose_solver == "after": + builder.use_solver(solver, **solver_kwargs) + builder.solve() + elapsed = time.perf_counter() - start + peak = tracemalloc.get_traced_memory()[1] + tracemalloc.stop() + + interpolator = builder.build() + timing = getattr(interpolator, "latest_solve_timing", None) or {} + return { + "builder": builder, + "interpolator": interpolator, + "field": interpolator.evaluate_value(points), + "seconds": elapsed, + "python peak MB": peak / 1e6, + "converged": bool(interpolator.up_to_date), + "iterations": timing.get("iterations"), + "unknowns": interpolator.dof, + } + + +# %% +# The Python solver as the reference +# ---------------------------------- + +python = run("cg") +print(f"{python['unknowns']} unknowns; Python CG: {python['seconds']:.1f} s") + +# %% +# The native solver with its defaults +# ----------------------------------- +# Same problem, native solver. It is faster, but look at ``converged`` and at the error. + +native_default = run("native") + +error = lambda result: np.sqrt(np.mean((result["field"] - truth(points)) ** 2)) +print(f"converged: {native_default['converged']} after {native_default['iterations']} iterations") +print(f"rms error against the truth: Python {error(python):.4f}, native {error(native_default):.4f}") +print(f"largest difference from the Python field: {np.max(np.abs(native_default['field'] - python['field'])):.2f}") + +# %% +# The native solver stopped at its default limit of 5000 iterations without converging, logged a +# warning, and left the interpolator marked as *not* up to date. The field it returned is an +# unfinished solution and is far from the Python one, even though the residual in the warning +# looks close to the tolerance (about 3e-8 against 1e-8): on this problem a small residual does not +# mean a small error in the field. Nothing raised, so the way to catch this is to read the warning +# or check ``interpolator.up_to_date``, not to judge by the residual. +# +# The limit is a keyword, ``maxiter``, passed through ``use_solver``. This problem needs a little +# more than 5000 iterations, so raising the limit is enough: + +native = run("native", maxiter=50000) +print(f"converged: {native['converged']} after {native['iterations']} iterations") +print(f"largest difference from the Python field: {np.max(np.abs(native['field'] - python['field'])):.1e}") +print(f"time: Python {python['seconds']:.1f} s, native {native['seconds']:.1f} s") + +# %% +# Now the native field agrees with the Python field to about 1e-8, which is the solver tolerance, +# and it is faster. How much faster depends on the machine and the problem size. In a separate run +# on a problem twice this size the gap was larger (about 3 times, against about 1.7 times here). + +# %% +# Memory: choose the solver before the setup +# ------------------------------------------ +# The Python memory used while setting up and solving, for the Python solver, for the native +# solver chosen before ``setup_interpolator``, and for the native solver chosen after it. + +after = run("native", choose_solver="after", maxiter=50000) +memory = pd.DataFrame( + { + "Python CG": python, + "native, chosen before setup": native, + "native, chosen after setup": after, + } +).T[["python peak MB", "seconds"]] +print(memory.astype(float).round(1).to_string()) +print( + "same field whichever way: largest difference between the two native runs " + f"{np.max(np.abs(after['field'] - native['field'])):.1e}" +) + +# %% +# Chosen before setup, the Python side never builds the sparse matrix and its peak memory is a +# small fraction of the Python solver's. Chosen after setup, the matrix has already been built, so +# the peak is many times higher than that (the solver logs a warning to say so), although still +# below the Python solver's here. The answer is the same either way. + +# %% +# Starting from the previous solution +# ----------------------------------- +# After a small change, for example a new regularisation weight, the previous field is a good +# starting point. ``warm_start=True`` starts the solver from it instead of from zero. + + +def resolve(warm_start): + builder = InterpolatorBuilder("FDI", bounding_box=bounding_box, nelements=50000, buffer=0.1) + builder.use_solver("native", maxiter=50000) + builder.add_value_constraints(values).add_normal_constraints(normals) + builder.setup_interpolator().solve() + first = builder.interpolator.latest_solve_timing["iterations"] + + builder.setup_interpolator(regularisation=1.1) # a slightly stronger smoothing weight + start = time.perf_counter() + builder.solve(warm_start=warm_start) + seconds = time.perf_counter() - start + return first, builder.interpolator.latest_solve_timing["iterations"], seconds + + +first, cold_iterations, cold_seconds = resolve(False) +_, warm_iterations, warm_seconds = resolve(True) +print(f"first solve: {first} iterations") +print(f"new weight, starting from zero: {cold_iterations} iterations, {cold_seconds:.2f} s") +print(f"new weight, warm start: {warm_iterations} iterations, {warm_seconds:.2f} s") + +# %% +# What to take from this +# ---------------------- +# * The native solver gives the same answer as the Python solver, to the solver tolerance, and +# is faster. The saving in Python memory only appears if the solver is chosen **before** +# ``setup_interpolator``. +# * **Check that it converged.** The default limit of 5000 iterations is not enough for a problem +# of this size, and what it returns when it stops is not the answer. Look for the warning, check +# ``interpolator.up_to_date``, and raise ``maxiter`` if needed. +# * ``warm_start=True`` saves iterations when you solve again after a small change. +# +# The timings and memory figures are from one machine and one problem; they show the shape of the +# behaviour, not a benchmark. For the full list of options and unsupported settings see the native +# solver page of the user guide. diff --git a/examples/README.rst b/examples/README.rst index 69348d2f5..3e32da5e6 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -14,7 +14,8 @@ learn them: logging, controlling data weighting, and comparing interpolators. 5. **Interpolation** - how the interpolators behave, measured against synthetic models with a known answer: choosing an interpolator, the regularisation weight, and - fold constraints with sparse data, refolded folds, and the kinds of orientation data. + fold constraints with sparse data, refolded folds, the kinds of orientation data, and + the native solver. Each example is a standalone, runnable Python script. Most examples load one of the sample datasets bundled in :code:`LoopStructural.datasets`, so From 0fa5424b6f22897d8222437126fbca0ab1cf48dc Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 17:28:07 +1030 Subject: [PATCH 06/10] ci: install the repo itself with the docs-check packages so pip does not resolve LoopStructural from PyPI Co-Authored-By: Claude Sonnet 5.5 --- .github/workflows/docs-check.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml index 4f4a23e67..4bb159acf 100644 --- a/.github/workflows/docs-check.yml +++ b/.github/workflows/docs-check.yml @@ -42,7 +42,7 @@ jobs: - name: Install the packages and the docs extra run: | - pip install -e packages/loop_common -e packages/loop_interpolation \ + pip install -e . -e packages/loop_common -e packages/loop_interpolation \ -e "packages/loopstructural_visualisation[all]" pip install -e ".[docs]" From 639ae96d5614a94e434164b0f706106ac4c2d532 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 17:58:41 +1030 Subject: [PATCH 07/10] ci: install the docs extra first, then the local packages over it The root pyproject pins loop-common and loop-interpolation below the versions in packages/, so installing everything in one resolve fails. Co-Authored-By: Claude Sonnet 5.5 --- .github/workflows/docs-check.yml | 6 ++++-- .github/workflows/documentation.yml | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml index 4bb159acf..19a23e9d1 100644 --- a/.github/workflows/docs-check.yml +++ b/.github/workflows/docs-check.yml @@ -42,9 +42,11 @@ jobs: - name: Install the packages and the docs extra run: | - pip install -e . -e packages/loop_common -e packages/loop_interpolation \ - -e "packages/loopstructural_visualisation[all]" + # The root pyproject pins loop-common/loop-interpolation below the versions in + # packages/, so install the docs extra first and then put the local packages over it. pip install -e ".[docs]" + pip install -e packages/loop_common -e packages/loop_interpolation \ + -e "packages/loopstructural_visualisation[all]" - name: Use the OSMesa build of VTK run: | diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 5763251b6..44243019f 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -42,6 +42,9 @@ jobs: - name: Install PyVista and dependencies run: | pip install -e .[docs] + # the root pyproject pins the packages below the versions in packages/; install those over it + pip install -e packages/loop_common -e packages/loop_interpolation \ + -e "packages/loopstructural_visualisation[all]" - name: Install custom OSMesa VTK variant run: | From f95929b6371d2f364fcc1df23e210dc342cdb8fe Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 18:00:17 +1030 Subject: [PATCH 08/10] docs: correct visualisation guide (naming, bounding_box, min/max, integer value) and remove stale LavaVu mentions Co-Authored-By: Claude Sonnet 5.5 --- docs/source/_static/wls_install.rst | 2 +- docs/source/user_guide/visualisation.rst | 23 +++++++++++++------ examples/1_basic/plot_1_data_preparation.py | 2 +- .../1_basic/plot_3_model_visualisation.py | 6 ++--- 4 files changed, 20 insertions(+), 13 deletions(-) diff --git a/docs/source/_static/wls_install.rst b/docs/source/_static/wls_install.rst index 63554ed72..2f0d4ef24 100644 --- a/docs/source/_static/wls_install.rst +++ b/docs/source/_static/wls_install.rst @@ -218,7 +218,7 @@ Windows using linux subsystem **Running LoopStructural** - To use the model viewing capabilities of LavaVu you need to use MobaXterm or another ssh/terminal client with x forwarding capabilities. + To use the interactive 3D model viewer you need to use MobaXterm or another ssh/terminal client with x forwarding capabilities. To run the included examples in LoopStructural you can simply run .. code-block:: diff --git a/docs/source/user_guide/visualisation.rst b/docs/source/user_guide/visualisation.rst index fbc92c480..428268863 100644 --- a/docs/source/user_guide/visualisation.rst +++ b/docs/source/user_guide/visualisation.rst @@ -22,8 +22,9 @@ Because ``Loop3DView`` *is* a PyVista plotter, everything PyVista can do (``add_ view.display() # opens the viewer (inline in a notebook) Some methods need a bounding box to know where to evaluate the feature. By default this is -the bounding box of the model passed to ``Loop3DView``; every ``plot_*`` method also accepts a -``bounding_box=`` argument to override it. +the bounding box of the model passed to ``Loop3DView``; ``plot_surface``, +``plot_scalar_field``, ``plot_vector_field``, ``plot_fault`` and ``plot_fault_displacements`` +also accept a ``bounding_box=`` argument to override it. Visualising geological features ------------------------------- @@ -37,7 +38,9 @@ There are four main ways to look at a geological feature. view.plot_scalar_field(feature, cmap="viridis") 2. Surfaces: isosurfaces of the implicit function. ``value`` is one isovalue, a list of - isovalues, or an integer number of regularly spaced surfaces. + isovalues, or a positive integer number of regularly spaced surfaces between the minimum and + maximum of the feature (``0`` is read as the isovalue 0.0, not a count). If ``value`` is + omitted the surface at the middle of the feature's range is drawn. .. code-block:: python @@ -67,7 +70,8 @@ The methods share a few arguments: * - ``opacity`` - transparency of the object, between 0 and 1 * - ``name`` - - name of the object in the viewer. Adding two objects with the same name replaces the first + - name of the object in the viewer. If the name is already used a numeric suffix is added + (``name_1``, ``name_2``, ...) so nothing is replaced * - ``pyvista_kwargs`` - a dictionary passed on to :meth:`pyvista.Plotter.add_mesh` * - ``bounding_box`` @@ -108,7 +112,9 @@ constraints as arrows. Arbitrary points and arrows can be added with the PyVista .. warning:: - Pass a distinct ``name`` for each object. Reusing a name replaces the earlier object. + Pass a distinct ``name`` for each object. ``Loop3DView`` does not replace an existing + object when a name is reused: the new object is given a numeric suffix (``name_1``, + ``name_2``, ...) and both stay in the scene. Visualising the geological model -------------------------------- @@ -153,8 +159,11 @@ location, so arbitrary quantities can be drawn with the same methods. view.plot_scalar_field(custom_feature, bounding_box=bb) The function is evaluated at the locations the visualisation needs, normally between the -origin and maximum of the bounding box. Set the ``min`` and ``max`` attributes of the feature -if the colour range should differ from the default. +origin and maximum of the bounding box. ``min()`` and ``max()`` on a feature are methods, not +settable attributes, and a feature without a model (such as a lambda feature) returns 0 for +both. To control the colour range pass ``vmin`` and ``vmax`` to the plotting method instead. +If you use the lambda feature as ``paint_with``, also pass ``vmin`` and ``vmax`` since the +default range comes from ``min()`` and ``max()``. Using the grid directly ----------------------- diff --git a/examples/1_basic/plot_1_data_preparation.py b/examples/1_basic/plot_1_data_preparation.py index 17fbf0dcc..a97b81712 100644 --- a/examples/1_basic/plot_1_data_preparation.py +++ b/examples/1_basic/plot_1_data_preparation.py @@ -134,7 +134,7 @@ # xyz coordinates. # # In the following example we will use matplotlib to visualise these results however, the -# next tutorial will show how to use the lavavu visualisation model. +# next tutorial will show how to use the 3D (PyVista based) visualisation module. import matplotlib.pyplot as plt diff --git a/examples/1_basic/plot_3_model_visualisation.py b/examples/1_basic/plot_3_model_visualisation.py index 9eb5318a1..ae0f17978 100644 --- a/examples/1_basic/plot_3_model_visualisation.py +++ b/examples/1_basic/plot_3_model_visualisation.py @@ -3,10 +3,8 @@ 1c. Visualising models =============================== The following tutorial will demonstrate how to use the Loop structural visualisation module. -This module provides a wrapper for the lavavu model that is written by -Owen Kaluza. - -Lavavu allows for interactive visualisation of 3D models within a jupyter +This module provides a wrapper around the PyVista plotter that allows +interactive visualisation of 3D models, including within a jupyter notebook environment. """ From 3b2970d951eff774c6e7cf1de31f0ce34a9807df Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 18:00:36 +1030 Subject: [PATCH 09/10] docs: note the expected fold_norm sign warning in the fold examples; number the interpolation gallery title Co-Authored-By: Claude Sonnet 5.5 --- examples/5_interpolation/README.rst | 4 ++-- examples/5_interpolation/plot_3_fold_constraints.py | 2 +- examples/5_interpolation/plot_4_refolded_folds.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/5_interpolation/README.rst b/examples/5_interpolation/README.rst index 62abebd2d..540174fc9 100644 --- a/examples/5_interpolation/README.rst +++ b/examples/5_interpolation/README.rst @@ -1,5 +1,5 @@ -Interpolation -------------- +5. Interpolation +---------------- These examples use synthetic geological models whose scalar field is known exactly, so every result can be measured against the truth. They show how the interpolators behave: which to diff --git a/examples/5_interpolation/plot_3_fold_constraints.py b/examples/5_interpolation/plot_3_fold_constraints.py index cff178fa5..0d125efad 100644 --- a/examples/5_interpolation/plot_3_fold_constraints.py +++ b/examples/5_interpolation/plot_3_fold_constraints.py @@ -84,7 +84,7 @@ def contact(x, y, level=0.5): # the axis. For this fold they are just ``x`` and ``y``. A ``FoldEvent`` combines the frame with # the fold axis and a function giving the limb rotation angle (the dip of the limb) in terms of # ``gx``. For a cylindrical fold the axis is constant. The signs follow the library's -# convention for the direction of the frame's third axis. +# convention for the direction of the frame's third axis. Note: the library may warn that ``fold_norm`` may need sign inversion; that warning is expected here because of this deliberate sign convention (the ``+ 180``/negated angles) and can be ignored. class Coordinate: diff --git a/examples/5_interpolation/plot_4_refolded_folds.py b/examples/5_interpolation/plot_4_refolded_folds.py index 4cec24b41..42b2bbf63 100644 --- a/examples/5_interpolation/plot_4_refolded_folds.py +++ b/examples/5_interpolation/plot_4_refolded_folds.py @@ -80,7 +80,7 @@ def contact(x, y, level=0.5): # # The axis rotation is the plunge of the early fold's axis, which the later fold varies with # ``gy``. The ``+ 180`` is a sign convention of the library: the frame's third axis must point -# against the orientation vectors. +# against the orientation vectors. Note: the library may warn that ``fold_norm`` may need sign inversion; that warning is expected here because of this deliberate sign convention (the ``+ 180``/negated angles) and can be ignored. class Coordinate: From 9528e8a8f30265852e181712dc70364abc83172b Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 29 Oct 2026 18:13:14 +1030 Subject: [PATCH 10/10] ci: install the docs dependencies with the workspace packages --no-deps The root pyproject pins loop-common and loop-interpolation below the versions in packages/, and pip also checks the pin when only the packages are installed afterwards. Install the third-party requirements first and the workspace packages without dependencies. Co-Authored-By: Claude Sonnet 5.5 --- .github/scripts/install_docs_deps.py | 40 ++++++++++++++++++++++++++++ .github/workflows/docs-check.yml | 6 +---- .github/workflows/documentation.yml | 5 +--- 3 files changed, 42 insertions(+), 9 deletions(-) create mode 100644 .github/scripts/install_docs_deps.py diff --git a/.github/scripts/install_docs_deps.py b/.github/scripts/install_docs_deps.py new file mode 100644 index 000000000..a08a17e2f --- /dev/null +++ b/.github/scripts/install_docs_deps.py @@ -0,0 +1,40 @@ +"""Install the documentation dependencies against the packages in this checkout. + +The root pyproject pins loop-common and loop-interpolation below the versions in packages/, and the +visualisation package depends on LoopStructural, so a single ``pip install -e .[docs]`` either pulls +the old releases from PyPI or fails to resolve. Install the third-party requirements from the +pyprojects, then install the workspace packages with --no-deps. +""" + +import re +import subprocess +import sys +import tomllib +from pathlib import Path + +ROOT = Path(__file__).resolve().parents[2] +WORKSPACE = ("loop-common", "loop-interpolation", "loopstructural", "loopstructuralvisualisation") +PACKAGES = ("packages/loop_common", "packages/loop_interpolation", "packages/loopstructural_visualisation") + + +def project(path): + with open(ROOT / path / "pyproject.toml", "rb") as f: + return tomllib.load(f)["project"] + + +def third_party(requirements): + """Drop the requirements that are provided by this checkout.""" + return [r for r in requirements if re.split(r"[\s\[<>=!~;]", r, maxsplit=1)[0].lower().replace("_", "-") not in WORKSPACE] + + +root = project(".") +requirements = third_party(root["dependencies"]) + third_party(root["optional-dependencies"]["docs"]) +for path in PACKAGES: + meta = project(path) + requirements += third_party(meta.get("dependencies", [])) + for extra in ("native", "all"): + requirements += third_party(meta.get("optional-dependencies", {}).get(extra, [])) + +pip = [sys.executable, "-m", "pip", "install"] +subprocess.check_call(pip + sorted(set(requirements))) +subprocess.check_call(pip + ["--no-deps", "-e", str(ROOT)] + [a for p in PACKAGES for a in ("-e", str(ROOT / p))]) diff --git a/.github/workflows/docs-check.yml b/.github/workflows/docs-check.yml index 19a23e9d1..6f83778ed 100644 --- a/.github/workflows/docs-check.yml +++ b/.github/workflows/docs-check.yml @@ -42,11 +42,7 @@ jobs: - name: Install the packages and the docs extra run: | - # The root pyproject pins loop-common/loop-interpolation below the versions in - # packages/, so install the docs extra first and then put the local packages over it. - pip install -e ".[docs]" - pip install -e packages/loop_common -e packages/loop_interpolation \ - -e "packages/loopstructural_visualisation[all]" + python .github/scripts/install_docs_deps.py - name: Use the OSMesa build of VTK run: | diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 44243019f..178c355fd 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -41,10 +41,7 @@ jobs: - name: Install PyVista and dependencies run: | - pip install -e .[docs] - # the root pyproject pins the packages below the versions in packages/; install those over it - pip install -e packages/loop_common -e packages/loop_interpolation \ - -e "packages/loopstructural_visualisation[all]" + python .github/scripts/install_docs_deps.py - name: Install custom OSMesa VTK variant run: |