From 4ab7406a01a85ba7d5b4eb87fff4b4353242f4b8 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Wed, 4 Nov 2026 13:01:30 +1030 Subject: [PATCH 01/14] chore(lint): fix baseline Python checks --- .../modelling/core/_model_feature_factory.py | 192 +++---- .../modelling/core/_model_relationships.py | 48 +- .../modelling/core/_model_serializer.py | 48 +- .../modelling/features/_analytical_feature.py | 4 +- .../_cross_product_geological_feature.py | 3 +- .../modelling/features/_feature_converters.py | 25 +- .../modelling/features/_geological_feature.py | 5 +- .../features/_lambda_geological_feature.py | 25 +- .../features/_projected_vector_feature.py | 3 +- .../builders/_folded_feature_builder.py | 14 +- .../builders/_geological_feature_builder.py | 38 +- .../builders/_structural_frame_builder.py | 13 +- LoopStructural/utils/dtm_creator.py | 5 +- .../_interpolator_builder.py | 8 +- .../tests/golden/conftest.py | 54 +- .../tests/golden/geomodels.py | 68 ++- .../loop_interpolation/tests/golden/render.py | 129 +++-- .../loop_interpolation/tests/golden/run.py | 33 +- .../tests/golden/test_golden.py | 41 +- .../tests/test_discrete_fold_interpolator.py | 7 +- .../tests/test_native_inequality.py | 8 +- .../tests/test_native_no_rows.py | 2 +- .../tests/test_native_parity.py | 10 +- .../tests/test_p2_fold_interpolator.py | 37 +- .../tests/test_quadrature.py | 11 +- .../tests/test_solver_strategy.py | 4 +- .../fixtures/make_fixture.py | 52 +- .../scripts/check_wasm_size.py | 2 + .../loopstructuralvisualisation/_2d_viewer.py | 6 +- .../loopstructuralvisualisation/_3d_viewer.py | 162 +++--- .../loopstructuralvisualisation/__init__.py | 4 +- .../_register_loop_ui.py | 1 - .../_stratigraphic_column.py | 7 +- .../trame/__init__.py | 2 +- .../trame/ui/vuetify3.py | 93 ++-- .../tests/test_plotting.py | 2 +- packages/map2loop/pyproject.toml | 7 + packages/map2loop/src/map2loop/__init__.py | 23 +- .../geodata_files/load_map2loop_data.py | 17 +- .../map2loop/src/map2loop/aus_state_urls.py | 15 +- packages/map2loop/src/map2loop/config.py | 35 +- .../src/map2loop/contact_extractor.py | 27 +- packages/map2loop/src/map2loop/data_checks.py | 495 +++++++++++------- .../src/map2loop/deformation_history.py | 48 +- .../src/map2loop/fault_orientation.py | 6 +- .../map2loop/src/map2loop/interpolators.py | 23 +- packages/map2loop/src/map2loop/logging.py | 4 +- packages/map2loop/src/map2loop/mapdata.py | 209 ++++---- packages/map2loop/src/map2loop/project.py | 231 ++++---- packages/map2loop/src/map2loop/sampler.py | 48 +- packages/map2loop/src/map2loop/sorter.py | 57 +- .../src/map2loop/stratigraphic_column.py | 8 +- .../src/map2loop/thickness_calculator.py | 187 ++++--- .../map2loop/src/map2loop/throw_calculator.py | 5 +- packages/map2loop/src/map2loop/topology.py | 38 +- .../map2loop/src/map2loop/utils/__init__.py | 28 +- .../src/map2loop/utils/load_map2loop_data.py | 16 +- .../src/map2loop/utils/utility_functions.py | 39 +- .../test_contact_extractor.py | 7 +- .../map2loop/tests/data_checks/test_config.py | 26 +- .../data_checks/test_input_data_faults.py | 7 +- .../tests/data_checks/test_input_data_fold.py | 79 ++- .../data_checks/test_input_data_geology.py | 17 +- .../data_checks/test_input_data_structure.py | 7 +- ..._mapdata_assign_random_colours_to_units.py | 1 + .../tests/mapdata/test_mapdata_dipdir.py | 18 +- .../mapdata/test_minimum_fault_length.py | 6 +- .../test_set_get_recreate_bounding_box.py | 7 +- .../test_set_get_working_projection.py | 13 +- .../tests/project/test_config_arguments.py | 62 ++- .../test_ignore_codes_setters_getters.py | 9 +- .../tests/project/test_plot_hamersley.py | 26 +- .../project/test_thickness_calculations.py | 18 +- .../tests/sampler/test_SamplerSpacing.py | 19 +- .../sampler/test_SamplerSpacing_featureId.py | 5 +- packages/map2loop/tests/test_import.py | 2 +- .../test_interpolated_structure.py | 12 +- .../test_ThicknessStructuralPoint.py | 4 +- .../test_ThicknessCalculatorAlpha.py | 6 +- .../map2loop/tests/topology/test_topology.py | 5 +- .../tests/utils/test_rgb_and_hex_functions.py | 5 +- 81 files changed, 1854 insertions(+), 1239 deletions(-) diff --git a/LoopStructural/modelling/core/_model_feature_factory.py b/LoopStructural/modelling/core/_model_feature_factory.py index 11bafc939..ac29bb074 100644 --- a/LoopStructural/modelling/core/_model_feature_factory.py +++ b/LoopStructural/modelling/core/_model_feature_factory.py @@ -13,22 +13,30 @@ from __future__ import annotations +from typing import TYPE_CHECKING, Callable + import numpy as np import pandas as pd from LoopStructural import LoopStructuralConfig -from ...modelling.features import FeatureType, UnconformityFeature +from ...modelling.features import FeatureType, GeologicalFeature, UnconformityFeature from ...modelling.features.builders import ( FaultBuilder, FoldedFeatureBuilder, GeologicalFeatureBuilder, StructuralFrameBuilder, ) +from ...modelling.features.fault import FaultSegment from ...modelling.features.fold import FoldEvent, FoldFrame -from ...modelling.intrusions import IntrusionBuilder, IntrusionFrameBuilder +from ...modelling.intrusions import IntrusionBuilder, IntrusionFeature, IntrusionFrameBuilder from ...utils import getLogger +if TYPE_CHECKING: + from numpy.typing import ArrayLike + + from .geological_model import GeologicalModel + logger = getLogger(__name__) # interpolator types that implement fold constraints @@ -38,18 +46,19 @@ class ModelFeatureFactory: @staticmethod def build_foliation( - model, + model: GeologicalModel, series_surface_name: str, *, index: int | None = None, data: pd.DataFrame | None = None, interpolatortype: str = "FDI", nelements: int = LoopStructuralConfig.nelements, - tol=None, - faults=None, - **kwargs, - ): - """ + tol: float | None = None, + faults: list[GeologicalFeature | str] | None = None, + **kwargs: object, + ) -> GeologicalFeature | None: + """Build a foliation feature from model data. + Parameters ---------- series_surface_name : string @@ -72,7 +81,7 @@ def build_foliation( the created geological feature Notes - ------ + ----- This function creates an instance of a :class:`LoopStructural.modelling.features.builders.GeologicalFeatureBuilder` and will return a :class:`LoopStructural.modelling.features.builders.GeologicalFeature` @@ -83,7 +92,6 @@ def build_foliation( An interpolator will be chosen by calling :meth:`LoopStructural.GeologicalModel.get_interpolator` """ - # if tol is not specified use the model default if tol is None: tol = model.tol @@ -103,7 +111,9 @@ def build_foliation( if data.shape[0] == 0: logger.warning("No data for {series_surface_data}, skipping") return - series_builder.add_data_from_data_frame(model.prepare_data(data, include_feature_name=False)) + series_builder.add_data_from_data_frame( + model.prepare_data(data, include_feature_name=False) + ) model._add_faults(series_builder, features=faults) # build feature @@ -118,18 +128,19 @@ def build_foliation( @staticmethod def build_fold_frame( - model, + model: GeologicalModel, fold_frame_name: str, *, index: int | None = None, - data=None, - interpolatortype="FDI", - nelements=LoopStructuralConfig.nelements, - tol=None, - buffer=0.1, - **kwargs, - ): - """ + data: pd.DataFrame | None = None, + interpolatortype: str = "FDI", + nelements: int = LoopStructuralConfig.nelements, + tol: float | None = None, + buffer: float = 0.1, + **kwargs: object, + ) -> FoldFrame | None: + """Build a fold frame from model data. + Parameters ---------- fold_frame_name : string @@ -156,7 +167,6 @@ def build_fold_frame( fold_frame : FoldFrame the created fold frame """ - if tol is None: tol = model.tol @@ -195,22 +205,22 @@ def build_fold_frame( @staticmethod def build_folded_foliation( - model, - foliation_name, + model: GeologicalModel, + foliation_name: str, *, index: int | None = None, - data=None, - interpolatortype="DFI", - nelements=LoopStructuralConfig.nelements, - buffer=0.1, - fold_frame=None, - svario=True, - tol=None, - invert_fold_norm=False, - **kwargs, - ): + data: pd.DataFrame | None = None, + interpolatortype: str = "DFI", + nelements: int = LoopStructuralConfig.nelements, + buffer: float = 0.1, + fold_frame: FoldFrame | None = None, + svario: bool = True, + tol: float | None = None, + invert_fold_norm: bool = False, + **kwargs: object, + ) -> GeologicalFeature | None: """ - Create a folded foliation field from data and a fold frame + Create a folded foliation field from data and a fold frame. Parameters ---------- @@ -236,7 +246,6 @@ def build_folded_foliation( Notes ----- - - The P2 mesh is larger than the model box because of ``buffer``. - Building a folded foliation uses the fold interpolation code from Laurent et al., 2016 and fold profile fitting from Grose et al., 2017. For more information about the fold modelling @@ -244,7 +253,6 @@ def build_folded_foliation( :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder` """ - if tol is None: tol = model.tol @@ -277,7 +285,9 @@ def build_folded_foliation( if data.shape[0] == 0: logger.warning(f"No data for {foliation_name}, skipping") return - series_builder.add_data_from_data_frame(model.prepare_data(data, include_feature_name=False)) + series_builder.add_data_from_data_frame( + model.prepare_data(data, include_feature_name=False) + ) model._add_faults(series_builder) # build feature @@ -294,18 +304,18 @@ def build_folded_foliation( @staticmethod def build_folded_fold_frame( - model, + model: GeologicalModel, fold_frame_name: str, *, index: int | None = None, data: pd.DataFrame | None = None, - interpolatortype="FDI", - nelements=LoopStructuralConfig.nelements, - fold_frame=None, - tol=None, - **kwargs, - ): - """ + interpolatortype: str = "FDI", + nelements: int = LoopStructuralConfig.nelements, + fold_frame: FoldFrame | None = None, + tol: float | None = None, + **kwargs: object, + ) -> FoldFrame: + """Build a folded fold frame from model data. Parameters ---------- @@ -348,7 +358,6 @@ def build_folded_fold_frame( see :class:`LoopStructural.modelling.features.fold.FoldEvent`, :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder` """ - if tol is None: tol = model.tol @@ -398,14 +407,15 @@ def build_folded_fold_frame( @staticmethod def validate_intrusion_inputs( - intrusion_name, - intrusion_frame_name, - intrusion_data, - intrusion_frame_data, - intrusion_frame_parameters, - ): - """Fail fast, at the `create_and_add_intrusion` boundary, with a clear - message naming the missing piece -- instead of a bare `KeyError` + intrusion_name: str, + intrusion_frame_name: str, + intrusion_data: pd.DataFrame, + intrusion_frame_data: pd.DataFrame, + intrusion_frame_parameters: dict[str, object], + ) -> None: + """Fail fast at the intrusion creation boundary with a clear error. + + This names the missing piece instead of raising a bare `KeyError` several calls deep inside `IntrusionFrameBuilder`/`IntrusionBuilder` once building has already started. See ``INTRUSIONS.md`` finding 4. """ @@ -442,17 +452,17 @@ def validate_intrusion_inputs( @staticmethod def build_intrusion( - model, - intrusion_name, - intrusion_frame_name, + model: GeologicalModel, + intrusion_name: str, + intrusion_frame_name: str, *, - intrusion_frame_parameters=None, - intrusion_lateral_extent_model=None, - intrusion_vertical_extent_model=None, - geometric_scaling_parameters=None, - **kwargs, - ): - """ + intrusion_frame_parameters: dict[str, object] | None = None, + intrusion_lateral_extent_model: Callable[..., object] | None = None, + intrusion_vertical_extent_model: Callable[..., object] | None = None, + geometric_scaling_parameters: dict[str, object] | None = None, + **kwargs: object, + ) -> IntrusionFeature: + """Build an intrusion feature from model data. Note ----- @@ -570,15 +580,16 @@ def build_intrusion( @staticmethod def build_domain_fault( - model, - fault_surface_data, + model: GeologicalModel, + fault_surface_data: pd.DataFrame | str, *, - nelements=LoopStructuralConfig.nelements, - interpolatortype="FDI", + nelements: int = LoopStructuralConfig.nelements, + interpolatortype: str = "FDI", index: int | None = None, - **kwargs, - ): - """ + **kwargs: object, + ) -> UnconformityFeature: + """Build a domain fault from model data. + Parameters ---------- fault_surface_data : string @@ -626,32 +637,33 @@ def build_domain_fault( @staticmethod def build_fault( - model, + model: GeologicalModel, fault_name: str, displacement: float, *, index: int | None = None, data: pd.DataFrame | None = None, - interpolatortype="FDI", - tol=None, - fault_slip_vector=None, - fault_normal_vector=None, - fault_center=None, - major_axis=None, - minor_axis=None, - intermediate_axis=None, - faultfunction="BaseFault", - faults=None, + interpolatortype: str = "FDI", + tol: float | None = None, + fault_slip_vector: ArrayLike | None = None, + fault_normal_vector: ArrayLike | None = None, + fault_center: ArrayLike | None = None, + major_axis: float | None = None, + minor_axis: float | None = None, + intermediate_axis: float | None = None, + faultfunction: str = "BaseFault", + faults: list[GeologicalFeature | str] | None = None, force_mesh_geometry: bool = False, points: bool = False, - fault_buffer=0.2, - fault_trace_anisotropy=0.0, - fault_dip=90, - fault_dip_anisotropy=0.0, - fault_pitch=None, - **kwargs, - ): - """ + fault_buffer: float = 0.2, + fault_trace_anisotropy: float = 0.0, + fault_dip: float = 90, + fault_dip_anisotropy: float = 0.0, + fault_pitch: float | None = None, + **kwargs: object, + ) -> FaultSegment: + """Build a fault from model data. + Parameters ---------- fault_name : string diff --git a/LoopStructural/modelling/core/_model_relationships.py b/LoopStructural/modelling/core/_model_relationships.py index 5076864d6..03ba1e6b3 100644 --- a/LoopStructural/modelling/core/_model_relationships.py +++ b/LoopStructural/modelling/core/_model_relationships.py @@ -10,7 +10,15 @@ feature-construction methods. """ -from ...modelling.features import FeatureType, UnconformityFeature +from __future__ import annotations + +from typing import TYPE_CHECKING + +from ...modelling.features import FeatureType, GeologicalFeature, UnconformityFeature +from ...modelling.features.builders import GeologicalFeatureBuilder, StructuralFrameBuilder + +if TYPE_CHECKING: + from .geological_model import GeologicalModel from ...utils import getLogger logger = getLogger(__name__) @@ -18,8 +26,12 @@ class FeatureRelationshipManager: @staticmethod - def add_faults(model, feature_builder, features=None): - """Adds all existing faults to a geological feature builder + def add_faults( + model: GeologicalModel, + feature_builder: GeologicalFeatureBuilder | StructuralFrameBuilder, + features: list[GeologicalFeature | str] | None = None, + ) -> None: + """Add all existing faults to a geological feature builder. Parameters ---------- @@ -38,9 +50,10 @@ def add_faults(model, feature_builder, features=None): feature_builder.add_fault(f) @staticmethod - def add_domain_fault_above(model, feature): - """ - Looks through the feature list and adds any domain faults to the feature. The domain fault masks everything + def add_domain_fault_above(model: GeologicalModel, feature: GeologicalFeatureBuilder) -> None: + """Add the active domain fault above the feature. + + The domain fault masks everything where the fault scalar field is < 0 as being active when added to feature. Parameters @@ -63,9 +76,10 @@ def add_domain_fault_above(model, feature): break @staticmethod - def add_domain_fault_below(model, domain_fault): - """ - Looks through the feature list and adds any the domain_fault to the features + def add_domain_fault_below(model: GeologicalModel, domain_fault: GeologicalFeature) -> None: + """Add the domain fault to features below it. + + The domain fault is added to features that already exist in the stack until an unconformity or another domain fault is reached (that boundary is not cropped). domain faults to the feature. The domain fault masks everything where the fault scalar field @@ -92,9 +106,10 @@ def add_domain_fault_below(model, domain_fault): f.add_region(region) @staticmethod - def add_unconformity_above(model, feature): - """ - Adds a region to the feature to prevent the value from being + def add_unconformity_above(model: GeologicalModel, feature: GeologicalFeature) -> None: + """Add a region that excludes values above the unconformity. + + This prevents the value from being interpolated where the unconformities exists above e.g. if there is another feature above and the unconformity is at 0 then the features added below (after) will only be visible where the @@ -105,7 +120,6 @@ def add_unconformity_above(model, feature): model : GeologicalModel feature - GeologicalFeature """ - if feature.type == FeatureType.FAULT: return for f in reversed(model.features): @@ -117,7 +131,9 @@ def add_unconformity_above(model, feature): break @staticmethod - def add_unconformity(model, feature, value, index=None): + def add_unconformity( + model: GeologicalModel, feature: GeologicalFeature, value: float, index: int | None = None + ) -> UnconformityFeature: """ Use an existing feature to add an unconformity to the model. @@ -163,7 +179,9 @@ def add_unconformity(model, feature, value, index=None): return uc_feature @staticmethod - def add_onlap_unconformity(model, feature, value, index=None): + def add_onlap_unconformity( + model: GeologicalModel, feature: GeologicalFeature, value: float, index: int | None = None + ) -> UnconformityFeature | None: """ Use an existing feature to add an unconformity to the model. diff --git a/LoopStructural/modelling/core/_model_serializer.py b/LoopStructural/modelling/core/_model_serializer.py index cbf4d1120..5989d5abb 100644 --- a/LoopStructural/modelling/core/_model_serializer.py +++ b/LoopStructural/modelling/core/_model_serializer.py @@ -10,8 +10,11 @@ delegate to the staticmethods here. """ +from __future__ import annotations + import json import pathlib +from typing import TYPE_CHECKING import pandas as pd @@ -21,12 +24,17 @@ from ..features.fault import FaultSegment from .stratigraphic_column import StratigraphicColumn +if TYPE_CHECKING: + from .geological_model import GeologicalModel + logger = getLogger(__name__) class ModelSerializer: @staticmethod - def feature_recipe_kind(feature): + def feature_recipe_kind( + feature: GeologicalFeature | StructuralFrame | UnconformityFeature | FaultSegment, + ) -> str: if isinstance(feature, GeologicalFeature): return "foliation" if isinstance(feature, StructuralFrame): @@ -38,7 +46,7 @@ def feature_recipe_kind(feature): return feature.__class__.__name__.lower() @staticmethod - def to_dict(model): + def to_dict(model: GeologicalModel) -> dict[str, object]: result = {} result["model"] = {} result["model"]["features"] = [f.name for f in model.features] @@ -47,7 +55,9 @@ def to_dict(model): return result @staticmethod - def to_recipe_dict(model, data_reference=None): + def to_recipe_dict( + model: GeologicalModel, data_reference: str | pathlib.Path | None = None + ) -> dict[str, object]: recipe = { "schema": "LoopStructural.GeologicalModelRecipe", "version": 1, @@ -83,7 +93,7 @@ def to_recipe_dict(model, data_reference=None): return recipe @staticmethod - def from_recipe_dict(cls, recipe): + def from_recipe_dict(model_class: type[GeologicalModel], recipe: object) -> GeologicalModel: if not isinstance(recipe, dict): raise TypeError("recipe must be a dictionary") @@ -94,7 +104,7 @@ def from_recipe_dict(cls, recipe): if not isinstance(bounding_box, BoundingBox): raise TypeError("recipe must include a bounding_box dictionary") - model = cls(bounding_box) + model = model_class(bounding_box) data_source = model_data.get("data_source") if isinstance(data_source, dict): @@ -152,22 +162,28 @@ def from_recipe_dict(cls, recipe): return model @staticmethod - def to_recipe_json(model, data_reference=None, indent=2): + def to_recipe_json( + model: GeologicalModel, data_reference: str | pathlib.Path | None = None, indent: int = 2 + ) -> str: recipe = ModelSerializer.to_recipe_dict(model, data_reference=data_reference) return json.dumps(recipe, indent=indent) @staticmethod - def from_recipe_json(cls, json_str): + def from_recipe_json(model_class: type[GeologicalModel], json_str: str) -> GeologicalModel: if not isinstance(json_str, str): raise TypeError("json_str must be a string") try: recipe = json.loads(json_str) except json.JSONDecodeError as e: raise TypeError(f"json_str is not valid JSON: {e}") - return ModelSerializer.from_recipe_dict(cls, recipe) + return ModelSerializer.from_recipe_dict(model_class, recipe) @staticmethod - def save_recipe(model, filename, data_reference=None): + def save_recipe( + model: GeologicalModel, + filename: str | pathlib.Path, + data_reference: str | pathlib.Path | None = None, + ) -> None: filename = pathlib.Path(filename) recipe = ModelSerializer.to_recipe_dict(model, data_reference=data_reference) with open(filename, "w") as f: @@ -175,17 +191,19 @@ def save_recipe(model, filename, data_reference=None): logger.info(f"Recipe saved to {filename}") @staticmethod - def load_recipe(cls, filename): + def load_recipe( + model_class: type[GeologicalModel], filename: str | pathlib.Path + ) -> GeologicalModel: filename = pathlib.Path(filename) if not filename.exists(): raise FileNotFoundError(f"Recipe file not found: {filename}") with open(filename, "r") as f: recipe = json.load(f) logger.info(f"Recipe loaded from {filename}") - return ModelSerializer.from_recipe_dict(cls, recipe) + return ModelSerializer.from_recipe_dict(model_class, recipe) @staticmethod - def to_file(model, file): + def to_file(model: GeologicalModel, file: str | pathlib.Path) -> None: try: import dill as pickle except ImportError: @@ -199,7 +217,9 @@ def to_file(model, file): logger.error("Error saving file") @staticmethod - def from_file(cls, file, allow_pickle: bool = True): + def from_file( + model_class: type[GeologicalModel], file: str | pathlib.Path, allow_pickle: bool = True + ) -> GeologicalModel | None: if not allow_pickle: raise LoopValueError( "Pickle-based loading is disabled (allow_pickle=False). " @@ -232,7 +252,7 @@ def from_file(cls, file, allow_pickle: bool = True): except Exception as e: logger.error(f"Failed to load model from {file}: {e}") raise LoopValueError(f"Failed to load model from {file}: {e}") from e - if isinstance(model, cls): + if isinstance(model, model_class): logger.info("GeologicalModel initialised from file") return model else: diff --git a/LoopStructural/modelling/features/_analytical_feature.py b/LoopStructural/modelling/features/_analytical_feature.py index 561d4f70a..ebeca9c25 100644 --- a/LoopStructural/modelling/features/_analytical_feature.py +++ b/LoopStructural/modelling/features/_analytical_feature.py @@ -2,9 +2,9 @@ import numpy as np -from ._base_geological_feature import BaseFeature -from . import FeatureType from ...utils import getLogger +from . import FeatureType +from ._base_geological_feature import BaseFeature logger = getLogger(__name__) diff --git a/LoopStructural/modelling/features/_cross_product_geological_feature.py b/LoopStructural/modelling/features/_cross_product_geological_feature.py index ce6622c29..e3ee6ca6d 100644 --- a/LoopStructural/modelling/features/_cross_product_geological_feature.py +++ b/LoopStructural/modelling/features/_cross_product_geological_feature.py @@ -1,10 +1,11 @@ """ """ + from __future__ import annotations import numpy as np -from ._base_geological_feature import BaseFeature from ...utils import getLogger +from ._base_geological_feature import BaseFeature logger = getLogger(__name__) diff --git a/LoopStructural/modelling/features/_feature_converters.py b/LoopStructural/modelling/features/_feature_converters.py index e60148746..49fbf9641 100644 --- a/LoopStructural/modelling/features/_feature_converters.py +++ b/LoopStructural/modelling/features/_feature_converters.py @@ -2,21 +2,24 @@ from LoopStructural.modelling.features.fold import FoldEvent, FoldFrame +class _InvalidFoldFrameError(TypeError, ValueError): + """Report an invalid fold frame and keep ValueError compatibility.""" + + def add_fold_to_feature(feature, fold_frame, **kwargs): if not isinstance(fold_frame, FoldFrame): - raise ValueError("fold_frame must be a FoldFrame instance") + raise _InvalidFoldFrameError("fold_frame must be a FoldFrame instance") - fold = FoldEvent(fold_frame, name=f"Fold_{feature.name}", invert_norm=kwargs.get('invert_fold_norm', False)) - - builder = FoldedFeatureBuilder.from_feature_builder( - feature.builder, - fold, - **kwargs + fold = FoldEvent( + fold_frame, name=f"Fold_{feature.name}", invert_norm=kwargs.get('invert_fold_norm', False) ) + + builder = FoldedFeatureBuilder.from_feature_builder(feature.builder, fold, **kwargs) feature = builder.feature feature.fold = fold return feature + def convert_feature_to_structural_frame(feature, **kwargs): """ Convert a geological feature to a structural frame by adding the feature to the frame @@ -25,7 +28,7 @@ def convert_feature_to_structural_frame(feature, **kwargs): ---------- feature : GeologicalFeature the geological feature to convert - + Returns ------- StructuralFrame @@ -33,9 +36,5 @@ def convert_feature_to_structural_frame(feature, **kwargs): """ builder = feature.builder - new_builder = StructuralFrameBuilder.from_feature_builder( - builder, - **kwargs - ) + new_builder = StructuralFrameBuilder.from_feature_builder(builder, **kwargs) return new_builder.frame - \ No newline at end of file diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py index 4660591b5..41c540bd9 100644 --- a/LoopStructural/modelling/features/_geological_feature.py +++ b/LoopStructural/modelling/features/_geological_feature.py @@ -3,6 +3,7 @@ This module contains classes for representing geometrical elements in geological models such as foliations, fault planes, and fold rotation angles. """ + from __future__ import annotations import numpy as np @@ -10,9 +11,9 @@ from LoopStructural.utils.maths import gradient_from_tetrahedron, regular_tetraherdron_for_points from ...geometry import ValuePoints, VectorPoints -from ._base_geological_feature import BaseFeature -from . import FeatureType from ...utils import LoopValueError, getLogger +from . import FeatureType +from ._base_geological_feature import BaseFeature logger = getLogger(__name__) diff --git a/LoopStructural/modelling/features/_lambda_geological_feature.py b/LoopStructural/modelling/features/_lambda_geological_feature.py index 13695fb2a..cc76fd789 100644 --- a/LoopStructural/modelling/features/_lambda_geological_feature.py +++ b/LoopStructural/modelling/features/_lambda_geological_feature.py @@ -1,6 +1,7 @@ """ Geological features """ + from __future__ import annotations from typing import Callable @@ -9,9 +10,9 @@ from LoopStructural.utils.maths import gradient_from_tetrahedron, regular_tetraherdron_for_points -from ._base_geological_feature import BaseFeature -from . import FeatureType from ...utils import LoopValueError, getLogger +from . import FeatureType +from ._base_geological_feature import BaseFeature logger = getLogger(__name__) @@ -50,7 +51,14 @@ def __init__( builder : optional the builder used to create this feature, by default None """ - BaseFeature.__init__(self, name, model, faults if faults is not None else [], regions if regions is not None else [], builder) + BaseFeature.__init__( + self, + name, + model, + faults if faults is not None else [], + regions if regions is not None else [], + builder, + ) self.type = FeatureType.LAMBDA self.function = function self.gradient_function = gradient_function @@ -95,10 +103,12 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray: pos = f.apply_to_points(pos, precomputed_gx=precomputed_gx.get(id(f))) if self.function is not None: - v[mask] = self.function(pos[mask,:]) + v[mask] = self.function(pos[mask, :]) return v - def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False,element_scale_parameter=None) -> np.ndarray: + def evaluate_gradient( + self, pos: np.ndarray, ignore_regions=False, element_scale_parameter=None + ) -> np.ndarray: """Evaluate the gradient of the underlying function at locations, applying any faults associated with this feature @@ -165,9 +175,7 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False,element_scale_ tetrahedron.shape ) - values = self.function(tetrahedron_faulted.reshape(-1, 3)).reshape( - (-1, 4) - ) + values = self.function(tetrahedron_faulted.reshape(-1, 3)).reshape((-1, 4)) v[mask, :] = gradient_from_tetrahedron(tetrahedron[mask, :, :], values[mask]) return v @@ -190,5 +198,6 @@ def copy(self, name: str | None = None): self.faults, self.builder, ) + def is_valid(self): return not (self.function is None and self.gradient_function is None) diff --git a/LoopStructural/modelling/features/_projected_vector_feature.py b/LoopStructural/modelling/features/_projected_vector_feature.py index 1b4cc4653..37a1cb20c 100644 --- a/LoopStructural/modelling/features/_projected_vector_feature.py +++ b/LoopStructural/modelling/features/_projected_vector_feature.py @@ -1,10 +1,11 @@ """ """ + from __future__ import annotations import numpy as np -from ._base_geological_feature import BaseFeature from ...utils import getLogger +from ._base_geological_feature import BaseFeature logger = getLogger(__name__) diff --git a/LoopStructural/modelling/features/builders/_folded_feature_builder.py b/LoopStructural/modelling/features/builders/_folded_feature_builder.py index d2c0a8f9e..6552eeb74 100644 --- a/LoopStructural/modelling/features/builders/_folded_feature_builder.py +++ b/LoopStructural/modelling/features/builders/_folded_feature_builder.py @@ -3,10 +3,10 @@ from ....geometry import BoundingBox from ....interpolators import InterpolatorType from ....modelling.features import FeatureType -from ._geological_feature_builder import GeologicalFeatureBuilder from ....modelling.features.fold.fold_function import FoldRotationType, get_fold_rotation_profile from ....utils import InterpolatorError, getLogger from ....utils._api_registry import public_api +from ._geological_feature_builder import GeologicalFeatureBuilder logger = getLogger(__name__) @@ -66,6 +66,7 @@ def __init__( self.svario = svario self.axis_profile_type = axis_profile_type self.limb_profile_type = limb_profile_type + @classmethod def from_feature_builder(cls, feature_builder, fold, **kwargs): """Create a FoldedFeatureBuilder from an existing feature builder""" @@ -83,10 +84,11 @@ def from_feature_builder(cls, feature_builder, fold, **kwargs): fold=fold, nelements=feature_builder.interpolator.n_elements, name=feature_builder.name, - **kwargs + **kwargs, ) builder.data = feature_builder.data return builder + @property def fold_axis_rotation(self): if self.fold.fold_axis_rotation is None: @@ -119,7 +121,9 @@ def set_fold_axis(self): fold_axis_rotation = get_fold_rotation_profile(self.axis_profile_type, far, fad) if "axis_function" in kwargs: # allow predefined function to be used - logger.error("axis_function is deprecated, use a specific fold rotation angle profile type") + logger.error( + "axis_function is deprecated, use a specific fold rotation angle profile type" + ) else: fold_axis_rotation.fit(params={'wavelength': kwargs.get("axis_wl", None)}) self.fold.fold_axis_rotation = fold_axis_rotation @@ -137,7 +141,9 @@ def set_fold_limb_rotation(self): fold_limb_rotation = get_fold_rotation_profile(self.limb_profile_type, flr, fld) if "limb_function" in kwargs: # allow for predefined functions to be used - logger.error("limb_function is deprecated, use a specific fold rotation angle profile type") + logger.error( + "limb_function is deprecated, use a specific fold rotation angle profile type" + ) else: fold_limb_rotation.fit(params={'wavelength': kwargs.get("limb_wl", None)}) diff --git a/LoopStructural/modelling/features/builders/_geological_feature_builder.py b/LoopStructural/modelling/features/builders/_geological_feature_builder.py index 11f1ea0cc..2ed03f535 100644 --- a/LoopStructural/modelling/features/builders/_geological_feature_builder.py +++ b/LoopStructural/modelling/features/builders/_geological_feature_builder.py @@ -7,7 +7,6 @@ from ....interpolators import DiscreteInterpolator, GeologicalInterpolator, InterpolatorFactory from ....modelling.features import GeologicalFeature -from ._base_builder import BaseBuilder from ....utils import getLogger from ....utils._api_registry import public_api from ....utils.helper import ( @@ -24,6 +23,7 @@ weight_name, xyz_names, ) +from ._base_builder import BaseBuilder logger = getLogger(__name__) @@ -57,7 +57,7 @@ def __init__( nelements=nelements, buffer=kwargs.get("buffer", 0.2), ) - + if not issubclass(type(interpolator), GeologicalInterpolator): raise TypeError( f"interpolator is {type(interpolator)} and must be a GeologicalInterpolator" @@ -84,7 +84,15 @@ def __init__( self._orthogonal_features = {} self._equality_constraints = {} # add default parameters - self.update_build_arguments({'cpw':1.0,'npw':1.0,'regularisation':1.,'nelements':self.interpolator.n_elements}) + self.update_build_arguments( + { + 'cpw': 1.0, + 'npw': 1.0, + 'regularisation': 1.0, + 'nelements': self.interpolator.n_elements, + } + ) + def set_not_up_to_date(self, caller): logger.info( f"Setting {self.name} to not up to date from an instance of {caller.__class__.__name__}" @@ -118,9 +126,7 @@ def add_data_from_data_frame(self, data_frame, overwrite=False): """ self.data = data_frame.copy() - def add_orthogonal_feature( - self, feature, w=1.0, region=None, step=1, B=0, sample_rule=None - ): + def add_orthogonal_feature(self, feature, w=1.0, region=None, step=1, B=0, sample_rule=None): """ Add a constraint to the interpolator so that the gradient of an exisitng feature is orthogonal to the feature being built. E.g. dot product @@ -429,6 +435,7 @@ def get_interface_constraints(self): ) else: return np.zeros((0, 5)) + # def get_inequality_pair_constraints(self): # mask = np.all(~np.isnan(self.data.loc[:,pairs_name()].to_numpy(float())),axis=1) # if mask.shape[0] > 0: @@ -504,7 +511,8 @@ def check_interpolation_geometry(self, data, buffer=0.3): ] self.interpolator.support.origin = origin self.interpolator.support.maximum = maximum - self.update_build_arguments({'nelements':self.interpolator.n_elements}) + self.update_build_arguments({'nelements': self.interpolator.n_elements}) + def build(self, data_region=None, **kwargs): """ Runs the interpolation and builds the geological feature @@ -533,7 +541,9 @@ def build(self, data_region=None, **kwargs): logger.info(f'Kwargs has {kwargs["nelements"]} elements') if self.interpolator.n_elements != kwargs['nelements']: logger.info('Setting nelements to {} for {}'.format(kwargs['nelements'], self.name)) - self.build_arguments['nelements'] = self.interpolator.set_nelements(kwargs['nelements']) + self.build_arguments['nelements'] = self.interpolator.set_nelements( + kwargs['nelements'] + ) logger.info(f'Interpolator nelements {self.interpolator.n_elements}') self._up_to_date = False if domain: @@ -545,24 +555,20 @@ def build(self, data_region=None, **kwargs): # if self.model.reuse_supports == False: if np.any(np.min(bb[0, :]) < self.interpolator.support.origin): neworigin = np.min([self.interpolator.support.origin, bb[0, :]], axis=0) - logger.info( - f"Changing origin of support for {self.name} from \ + logger.info(f"Changing origin of support for {self.name} from \ {self.interpolator.support.origin[0]} \ {self.interpolator.support.origin[1]} \ {self.interpolator.support.origin[2]} \ - to {neworigin[0]} {neworigin[1]} {neworigin[2]}" - ) + to {neworigin[0]} {neworigin[1]} {neworigin[2]}") self.interpolator.support.origin = neworigin if np.any(np.max(bb[0, :]) < self.interpolator.support.maximum): newmax = np.max([self.interpolator.support.maximum, bb[0, :]], axis=0) - logger.info( - f"Changing origin of support for {self.name} from \ + logger.info(f"Changing origin of support for {self.name} from \ from {self.interpolator.support.maximum[0]} \ {self.interpolator.support.maximum[1]} \ {self.interpolator.support.maximum[2]} \ - to {newmax[0]} {newmax[1]} {newmax[2]}" - ) + to {newmax[0]} {newmax[1]} {newmax[2]}") self.interpolator.support.maximum = newmax self.interpolator.set_region(region=region) diff --git a/LoopStructural/modelling/features/builders/_structural_frame_builder.py b/LoopStructural/modelling/features/builders/_structural_frame_builder.py index 4cf81d4b1..e925ba210 100644 --- a/LoopStructural/modelling/features/builders/_structural_frame_builder.py +++ b/LoopStructural/modelling/features/builders/_structural_frame_builder.py @@ -1,6 +1,7 @@ """ structural frame builder """ + from __future__ import annotations import copy @@ -18,9 +19,9 @@ from ....modelling.features import StructuralFrame +from ._base_builder import BaseBuilder from ._folded_feature_builder import FoldedFeatureBuilder from ._geological_feature_builder import GeologicalFeatureBuilder -from ._base_builder import BaseBuilder class StructuralFrameBuilder(BaseBuilder): @@ -144,13 +145,13 @@ def from_feature_builder(cls, feature_builder, **kwargs): if hasattr(feature_builder, 'fold'): logger.warning("feature builder has a fold - using this to create a folded frame") kwargs['fold'] = copy.deepcopy(feature_builder.fold) - builder = cls( - interpolatortype=[feature_builder.interpolator.type]*3, + builder = cls( + interpolatortype=[feature_builder.interpolator.type] * 3, bounding_box=feature_builder.model.bounding_box, - nelements=[feature_builder.interpolator.n_elements]*3, + nelements=[feature_builder.interpolator.n_elements] * 3, name=feature_builder.name, model=feature_builder.model, - **kwargs + **kwargs, ) builder.add_data_from_data_frame(feature_builder.data) return builder @@ -158,9 +159,11 @@ def from_feature_builder(cls, feature_builder, **kwargs): @property def build_arguments(self): return self.builders[0].build_arguments + def update_build_arguments(self, kwargs): for i in range(3): self.builders[i].update_build_arguments(kwargs) + @property def frame(self): return self._frame diff --git a/LoopStructural/utils/dtm_creator.py b/LoopStructural/utils/dtm_creator.py index 801980d7b..4edb2861c 100644 --- a/LoopStructural/utils/dtm_creator.py +++ b/LoopStructural/utils/dtm_creator.py @@ -1,12 +1,13 @@ +from __future__ import annotations + from pathlib import Path -from typing import Union from .logging import getLogger logger = getLogger(__name__) -def create_dtm_with_rasterio(dtm_path: Union[str, Path]): +def create_dtm_with_rasterio(dtm_path: str | Path): try: import rasterio except ImportError: diff --git a/packages/loop_interpolation/src/loop_interpolation/_interpolator_builder.py b/packages/loop_interpolation/src/loop_interpolation/_interpolator_builder.py index 556681972..df65e7130 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_interpolator_builder.py +++ b/packages/loop_interpolation/src/loop_interpolation/_interpolator_builder.py @@ -4,13 +4,15 @@ InterpolatorFactory and provides a chainable API for adding constraints, configuring setup options, and solving. """ + from __future__ import annotations import numpy as np from loop_common.geometry import BoundingBox -from ._interpolatortype import InterpolatorType +from ._geological_interpolator import GeologicalInterpolator from ._interpolator_factory import InterpolatorFactory +from ._interpolatortype import InterpolatorType class InterpolatorBuilder: @@ -228,9 +230,7 @@ def add_tangent_constraints(self, tangent_constraints: np.ndarray) -> Interpolat """ return self._set_constraint("set_tangent_constraints", tangent_constraints) - def add_inequality_constraints( - self, inequality_constraints: np.ndarray - ) -> InterpolatorBuilder: + def add_inequality_constraints(self, inequality_constraints: np.ndarray) -> InterpolatorBuilder: return self._set_constraint("set_value_inequality_constraints", inequality_constraints) def add_inequality_pair_constraints( diff --git a/packages/loop_interpolation/tests/golden/conftest.py b/packages/loop_interpolation/tests/golden/conftest.py index 5d94c5f24..9633b690c 100644 --- a/packages/loop_interpolation/tests/golden/conftest.py +++ b/packages/loop_interpolation/tests/golden/conftest.py @@ -1,5 +1,3 @@ - - def pytest_addoption(parser): parser.addoption( "--update-golden", @@ -16,16 +14,18 @@ def pytest_addoption(parser): # --------------------------------------------------------------------------- HTML test report -import base64 # noqa: E402 -import html as _html # noqa: E402 -import time # noqa: E402 +import base64 +import html as _html +import time _RECORDS: list[dict] = [] _FOLD_ENGINES = [("fold_sparse", e, False) for e in ("FD", "P1", "DFI", "FDFold", "P2Fold")] _BY_NAME = { "test_fold_constraints_help": _FOLD_ENGINES, - "test_refolded_fold_constraints_help": [("refolded_sparse", e, False) for e in ("FD", "P1", "DFI", "FDFold", "P2Fold")], + "test_refolded_fold_constraints_help": [ + ("refolded_sparse", e, False) for e in ("FD", "P1", "DFI", "FDFold", "P2Fold") + ], "test_p2_not_worse_than_p1_on_smooth_dome": [("dome", "P1", False), ("dome", "P2", False)], "test_known_p2_normal_weight_has_effect": [("planar_noisy", "P2", False)], "test_dfi_is_repeatable_with_reseed": [("fold_sparse", "DFI", False)], @@ -42,7 +42,10 @@ def _pairs_for(item): return [(params["case"], params["engine"], False)] if name == "test_known_default_solve_is_converged": return [("planar_dip", params["engine"], True), ("planar_dip", params["engine"], False)] - if name in ("test_known_update_keeps_fold_weights", "test_known_points_on_lower_boundary_are_finite"): + if name in ( + "test_known_update_keeps_fold_weights", + "test_known_points_on_lower_boundary_are_finite", + ): return [("fold_sparse", params["engine"], False)] return [] @@ -58,14 +61,27 @@ def pytest_runtest_makereport(item, call): xfail = item.get_closest_marker("xfail") if xfail is not None and call.excinfo.errisinstance(Exception) and call.when == "call": outcome, detail = "xfailed (known bug)", xfail.kwargs.get("reason", "") - elif call.excinfo.errisinstance(__import__("_pytest.outcomes", fromlist=["Skipped"]).Skipped): + elif call.excinfo.errisinstance( + __import__("_pytest.outcomes", fromlist=["Skipped"]).Skipped + ): outcome, detail = "skipped", str(call.excinfo.value) else: outcome, detail = "FAILED", str(call.excinfo.value)[:1500] - elif item.get_closest_marker("xfail") is not None and item.get_closest_marker("xfail").kwargs.get("strict"): - outcome, detail = "XPASS (bug fixed? remove the marker)", item.get_closest_marker("xfail").kwargs.get("reason", "") - _RECORDS.append({"item": item, "outcome": outcome, "detail": detail, "t": call.duration, - "doc": (item.function.__doc__ or "").strip()}) + elif item.get_closest_marker("xfail") is not None and item.get_closest_marker( + "xfail" + ).kwargs.get("strict"): + outcome, detail = "XPASS (bug fixed? remove the marker)", item.get_closest_marker( + "xfail" + ).kwargs.get("reason", "") + _RECORDS.append( + { + "item": item, + "outcome": outcome, + "detail": detail, + "t": call.duration, + "doc": (item.function.__doc__ or "").strip(), + } + ) def _b64(data: bytes) -> str: @@ -76,7 +92,6 @@ def pytest_sessionfinish(session, exitstatus): if not session.config.getoption("--report") or not _RECORDS: return import tempfile - from pathlib import Path from . import harness, render, test_golden @@ -105,7 +120,9 @@ def images(case, engine, default): figs = "" for case, engine, default in _pairs_for(item): surf, sec, m = images(case, engine, default) - label = f"{case} / {engine}" + (" / package-default solver" if default else " / converged solver") + label = f"{case} / {engine}" + ( + " / package-default solver" if default else " / converged solver" + ) figs += ( f'
{_html.escape(label)} — value rms {m["value_rms"]:.3g}, ' f'gradient error median {m["angle_median"]:.1f} deg
' @@ -116,7 +133,8 @@ def images(case, engine, default): f'{_html.escape(item.name)} {rec["t"]:.2f}s' + (f"

{_html.escape(rec['doc'])}

" if rec["doc"] else "") + (f'
{_html.escape(rec["detail"])}
' if rec["detail"] else "") - + figs + "" + + figs + + "" ) summary = ", ".join(f"{v} {k}" for k, v in sorted(counts.items())) legend = ( @@ -134,8 +152,10 @@ def images(case, engine, default): "img.s{width:520px;vertical-align:top}img.x{width:min(1100px,calc(100% - 540px));vertical-align:top}" "@media(max-width:1200px){img.s,img.x{width:100%}}" f"

Interpolation golden test report

{time.strftime('%Y-%m-%d %H:%M')} — {summary}

" - + legend + "".join(blocks) - + "" ) path = harness.OUT_DIR / "test_report.html" diff --git a/packages/loop_interpolation/tests/golden/geomodels.py b/packages/loop_interpolation/tests/golden/geomodels.py index eeb2105d8..8b91666ac 100644 --- a/packages/loop_interpolation/tests/golden/geomodels.py +++ b/packages/loop_interpolation/tests/golden/geomodels.py @@ -5,6 +5,7 @@ analytic field, so the error of a solution can be measured against the truth. All random numbers use fixed seeds: the data are identical on every run. """ + from __future__ import annotations from dataclasses import dataclass, field @@ -52,7 +53,7 @@ def _horizon_points(truth, xy: np.ndarray, levels=HORIZONS) -> np.ndarray: out = [] for x, y in xy: - def f(z, level=None): + def f(z, level=None, x=x, y=y): return truth(np.array([[x, y, z]]))[0] - level for level in levels: @@ -83,7 +84,10 @@ def planar_dip() -> Case: "Planar dipping contacts (dip 30, 3 horizons, 6 boreholes)", "Linear truth. Every interpolator must reproduce it. A failure here is a basic bug " "(sign, scaling, weights, data placement).", - s, g, values, np.c_[npts, g(npts)], + s, + g, + values, + np.c_[npts, g(npts)], ) @@ -98,7 +102,10 @@ def planar_extrapolation() -> Case: "Planar contacts, data only in x < 0.35", "Tests regularisation as an extrapolator into a data-free region (x > 0.65). " "The regularisation should continue the plane, not flatten or curl it.", - s, g, values, np.c_[npts, g(npts)], + s, + g, + values, + np.c_[npts, g(npts)], ) @@ -134,16 +141,28 @@ def g(p): rng = np.random.default_rng(4) 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))] + [ + [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 = _horizon_points(s, xy) - npts = 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]]) + npts = 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]] + ) return Case( "dome", "Gaussian dome (relief 0.2, sigma 0.2, 3 horizons, 9 boreholes)", "Smooth curved truth. Tests whether regularisation flattens the crest and whether " "gradient data bend the horizons correctly.", - s, g, values, np.c_[npts, g(npts)], + s, + g, + values, + np.c_[npts, g(npts)], ) @@ -188,7 +207,11 @@ def fold_section() -> Case: "Plain interpolation of a fold from one cross-section. Along the axis (y) the " "regularisation must carry the section outward. Error far from y = 0.5 is expected; " "the golden file detects a change in the pattern.", - s, g, values, np.c_[npts, g(npts)], fold=fold, + s, + g, + values, + np.c_[npts, g(npts)], + fold=fold, ) @@ -203,7 +226,12 @@ def fold_sparse() -> Case: "The Bermagui situation: little data, strong fold. The fold constraints must supply the " "gradient direction everywhere. The FD and P1 columns have no fold constraint and show " "what is gained. A drop in DFI/FDFold accuracy is a fold-constraint regression.", - s, g, values, np.c_[p, g(p)], fold=fold, engines=("FD", "P1", "P2", "DFI", "FDFold", "P2Fold"), + s, + g, + values, + np.c_[p, g(p)], + fold=fold, + engines=("FD", "P1", "P2", "DFI", "FDFold", "P2Fold"), ) @@ -226,7 +254,9 @@ def s(p): return p[:, 2] - a2 * np.sin(k2 * p[:, 0]) - b3 * np.sin(k3 * p[:, 1]) def g(p): - return np.c_[-a2 * k2 * np.cos(k2 * p[:, 0]), -b3 * k3 * np.cos(k3 * p[:, 1]), np.ones(len(p))] + return np.c_[ + -a2 * k2 * np.cos(k2 * p[:, 0]), -b3 * k3 * np.cos(k3 * p[:, 1]), np.ones(len(p)) + ] def limb(gx): return -np.rad2deg(np.arctan(a2 * k2 * np.cos(k2 * gx))) @@ -245,11 +275,27 @@ def axis_rotation(gy): "The Bermagui geometry in miniature: the hinge plunge changes along the axis, so a constant " "axis is wrong. Tests the fold_axis_rotation path of the fold interpolators (DFI, FDFold). " "FD, P1 and P2 have no fold constraint and show what the fold event adds.", - s, g, values, np.c_[p, g(p)], fold=fold, engines=("FD", "P1", "P2", "DFI", "FDFold", "P2Fold"), + s, + g, + values, + np.c_[p, g(p)], + fold=fold, + engines=("FD", "P1", "P2", "DFI", "FDFold", "P2Fold"), ) -ALL = {f.__name__: f for f in (planar_dip, planar_extrapolation, planar_noisy, dome, fold_section, fold_sparse, refolded_sparse)} +ALL = { + f.__name__: f + for f in ( + planar_dip, + planar_extrapolation, + planar_noisy, + dome, + fold_section, + fold_sparse, + refolded_sparse, + ) +} def build(name: str) -> Case: diff --git a/packages/loop_interpolation/tests/golden/render.py b/packages/loop_interpolation/tests/golden/render.py index 846c205e2..7bbda9694 100644 --- a/packages/loop_interpolation/tests/golden/render.py +++ b/packages/loop_interpolation/tests/golden/render.py @@ -7,19 +7,19 @@ 3. Section: unit gradient arrows, solved (blue) over truth (grey). 4. Plan (XY at z = 0.5): signed value error. Blue = solved too low, red = too high. """ + from __future__ import annotations -import textwrap from pathlib import Path import matplotlib matplotlib.use("Agg") -import matplotlib.pyplot as plt # noqa: E402 -import numpy as np # noqa: E402 +import matplotlib.pyplot as plt +import numpy as np -from . import geomodels # noqa: E402 -from .harness import OUT_DIR, angle_deg # noqa: E402 +from . import geomodels +from .harness import OUT_DIR, angle_deg N = 61 # section resolution SECTION_Y = 0.5 @@ -40,8 +40,15 @@ def _plane_grid(kind: str): def _horizons(ax, a, b, field, colors, styles, widths=1.2): if np.isfinite(field).sum() < 4: return - ax.contour(a, b, field.reshape(a.shape), levels=geomodels.HORIZONS, colors=colors, - linestyles=styles, linewidths=widths) + ax.contour( + a, + b, + field.reshape(a.shape), + levels=geomodels.HORIZONS, + colors=colors, + linestyles=styles, + linewidths=widths, + ) def _data(ax, case, kind): @@ -53,7 +60,9 @@ def _data(ax, case, kind): if mn.any(): nn = n[mn] u = nn[:, 3:6] / np.linalg.norm(nn[:, 3:6], axis=1, keepdims=True) - ax.quiver(nn[:, 0], nn[:, 2], u[:, 0], u[:, 2], color="k", scale=22, width=0.006, zorder=6) + ax.quiver( + nn[:, 0], nn[:, 2], u[:, 0], u[:, 2], color="k", scale=22, width=0.006, zorder=6 + ) ax.plot(n[~mn, 0], n[~mn, 2], "^", ms=3, mfc="none", mec="0.5", mew=0.6, zorder=5) else: mv, mn = np.abs(v[:, 2] - PLAN_Z) < SLAB, np.abs(n[:, 2] - PLAN_Z) < SLAB @@ -61,15 +70,17 @@ def _data(ax, case, kind): ax.plot(n[mn, 0], n[mn, 1], "^", ms=4, mfc="white", mec="k", mew=0.8, zorder=5) -def render_case(case_name: str, results: dict[str, dict], out_dir: Path | None = None, - suffix: str = "") -> Path: +def render_case( + case_name: str, results: dict[str, dict], out_dir: Path | None = None, suffix: str = "" +) -> Path: """`results` maps engine name -> harness.compute() result.""" case = geomodels.build(case_name) out_dir = Path(out_dir or OUT_DIR / "figures") out_dir.mkdir(parents=True, exist_ok=True) engines = list(results) - fig, axes = plt.subplots(len(engines), 4, figsize=(15, 3.55 * len(engines)), squeeze=False, - constrained_layout=True) + fig, axes = plt.subplots( + len(engines), 4, figsize=(15, 3.55 * len(engines)), squeeze=False, constrained_layout=True + ) xa, za, xz = _plane_grid("xz") xb, yb, xy = _plane_grid("xy") truth_xz, tgrad_xz = case.truth(xz), case.truth_grad(xz) @@ -84,8 +95,14 @@ def render_case(case_name: str, results: dict[str, dict], out_dir: Path | None = v_xy = np.asarray(interp.evaluate_value(xy), float) ax1, ax2, ax3, ax4 = axes[row] - im = ax1.contourf(xa, za, v_xz.reshape(xa.shape), levels=np.linspace(fmin, fmax, 21), - cmap="cividis", extend="both") + im = ax1.contourf( + xa, + za, + v_xz.reshape(xa.shape), + levels=np.linspace(fmin, fmax, 21), + cmap="cividis", + extend="both", + ) _horizons(ax1, xa, za, truth_xz, "k", "--", 1.4) _horizons(ax1, xa, za, v_xz, "#ff7f0e", "-", 1.6) _data(ax1, case, "xz") @@ -94,17 +111,34 @@ def render_case(case_name: str, results: dict[str, dict], out_dir: Path | None = ang = angle_deg(g_xz, tgrad_xz).reshape(xa.shape) im2 = ax2.pcolormesh(xa, za, ang, cmap="magma_r", vmin=0, vmax=45, shading="auto") _horizons(ax2, xa, za, truth_xz, "k", "--", 1.0) - ax2.set_title(f"gradient error median {m['angle_median']:.1f} p90 {m['angle_p90']:.1f} deg", fontsize=9) + ax2.set_title( + f"gradient error median {m['angle_median']:.1f} p90 {m['angle_p90']:.1f} deg", + fontsize=9, + ) gu = g_xz / np.linalg.norm(g_xz, axis=1, keepdims=True) tu = tgrad_xz / np.linalg.norm(tgrad_xz, axis=1, keepdims=True) _horizons(ax3, xa, za, truth_xz, "0.6", "-", 0.8) sx, sz = xa[arrow, arrow], za[arrow, arrow] shape = xa.shape - ax3.quiver(sx, sz, tu[:, 0].reshape(shape)[arrow, arrow], tu[:, 2].reshape(shape)[arrow, arrow], - color="0.55", scale=18, width=0.007) - ax3.quiver(sx, sz, gu[:, 0].reshape(shape)[arrow, arrow], gu[:, 2].reshape(shape)[arrow, arrow], - color="#1f77b4", scale=18, width=0.004) + ax3.quiver( + sx, + sz, + tu[:, 0].reshape(shape)[arrow, arrow], + tu[:, 2].reshape(shape)[arrow, arrow], + color="0.55", + scale=18, + width=0.007, + ) + ax3.quiver( + sx, + sz, + gu[:, 0].reshape(shape)[arrow, arrow], + gu[:, 2].reshape(shape)[arrow, arrow], + color="#1f77b4", + scale=18, + width=0.004, + ) ax3.set_title("unit gradient: solved (blue) over truth (grey)", fontsize=9) e = (v_xy - truth_xy).reshape(xb.shape) @@ -133,8 +167,14 @@ def render_case(case_name: str, results: dict[str, dict], out_dir: Path | None = return path -def render_diff(case_name: str, engine: str, points: np.ndarray, golden: dict, current: dict, - out_dir: Path | None = None) -> Path: +def render_diff( + case_name: str, + engine: str, + points: np.ndarray, + golden: dict, + current: dict, + out_dir: Path | None = None, +) -> Path: """Where did the field move? Maps the |current - golden| value and angle change on two sections.""" out_dir = Path(out_dir or OUT_DIR / "failures") out_dir.mkdir(parents=True, exist_ok=True) @@ -144,15 +184,26 @@ def render_diff(case_name: str, engine: str, points: np.ndarray, golden: dict, c da = angle_deg(current["gradient"], golden["gradient"]).reshape(shape) fig, axes = plt.subplots(2, 3, figsize=(12, 7), constrained_layout=True) i = n // 2 - views = {"XZ (y mid)": (dv[:, i, :], da[:, i, :]), "XY (z mid)": (dv[i, :, :], da[i, :, :]), - "YZ (x mid)": (dv[:, :, i], da[:, :, i])} + views = { + "XZ (y mid)": (dv[:, i, :], da[:, i, :]), + "XY (z mid)": (dv[i, :, :], da[i, :, :]), + "YZ (x mid)": (dv[:, :, i], da[:, :, i]), + } for col, (name, (v2, a2)) in enumerate(views.items()): lim = max(1e-12, float(np.nanmax(np.abs(v2)))) - p = axes[0, col].imshow(v2, origin="lower", extent=(0, 1, 0, 1), cmap="RdBu_r", vmin=-lim, vmax=lim) + p = axes[0, col].imshow( + v2, origin="lower", extent=(0, 1, 0, 1), cmap="RdBu_r", vmin=-lim, vmax=lim + ) axes[0, col].set_title(f"{name}: value change (current - golden)") fig.colorbar(p, ax=axes[0, col], shrink=0.8) - q = axes[1, col].imshow(a2, origin="lower", extent=(0, 1, 0, 1), cmap="magma_r", vmin=0, - vmax=max(1e-6, float(np.nanmax(a2)))) + q = axes[1, col].imshow( + a2, + origin="lower", + extent=(0, 1, 0, 1), + cmap="magma_r", + vmin=0, + vmax=max(1e-6, float(np.nanmax(a2))), + ) axes[1, col].set_title(f"{name}: gradient direction change (deg)") fig.colorbar(q, ax=axes[1, col], shrink=0.8) fig.suptitle(f"{case_name} / {engine}: current vs golden", fontsize=11) @@ -174,7 +225,9 @@ def render_surfaces(case_name: str, engine: str, interp, n: int = 41) -> bytes: u = np.linspace(eps, 1 - eps, n) zz, yy, xx = np.meshgrid(u, u, u, indexing="ij") pts = np.c_[xx.ravel(), yy.ravel(), zz.ravel()] - grid = pv.ImageData(dimensions=(n, n, n), spacing=((1 - 2 * eps) / (n - 1),) * 3, origin=(eps,) * 3) + grid = pv.ImageData( + dimensions=(n, n, n), spacing=((1 - 2 * eps) / (n - 1),) * 3, origin=(eps,) * 3 + ) grid["solved"] = np.asarray(interp.evaluate_value(pts), float) grid["truth"] = case.truth(pts) @@ -189,12 +242,22 @@ def render_surfaces(case_name: str, engine: str, interp, n: int = 41) -> bytes: truth = grid.contour([level], scalars="truth") if truth.n_points: pl.add_mesh(truth, color="black", style="wireframe", opacity=0.3, line_width=1) - pl.add_mesh(pv.PolyData(case.values[:, :3]), color="white", point_size=9, - render_points_as_spheres=True) - pl.add_mesh(pv.PolyData(case.values[:, :3]), color="black", point_size=11, - render_points_as_spheres=True, opacity=0.35) - pl.add_arrows(case.normals[:, :3], case.normals[:, 3:6] / np.linalg.norm(case.normals[:, 3:6], axis=1, keepdims=True), - mag=0.12, color="red") + pl.add_mesh( + pv.PolyData(case.values[:, :3]), color="white", point_size=9, render_points_as_spheres=True + ) + pl.add_mesh( + pv.PolyData(case.values[:, :3]), + color="black", + point_size=11, + render_points_as_spheres=True, + opacity=0.35, + ) + pl.add_arrows( + case.normals[:, :3], + case.normals[:, 3:6] / np.linalg.norm(case.normals[:, 3:6], axis=1, keepdims=True), + mag=0.12, + color="red", + ) pl.camera_position = [(2.1, -1.7, 1.5), (0.5, 0.5, 0.5), (0, 0, 1)] img = pl.screenshot(return_img=True) pl.close() diff --git a/packages/loop_interpolation/tests/golden/run.py b/packages/loop_interpolation/tests/golden/run.py index 05d298d5f..4f00f7744 100644 --- a/packages/loop_interpolation/tests/golden/run.py +++ b/packages/loop_interpolation/tests/golden/run.py @@ -1,22 +1,29 @@ """Command line for the golden suite. - python -m tests.golden.run run everything, write figures + report.html to out/ - python -m tests.golden.run --update also (re)write the golden files from the current code - python -m tests.golden.run --case dome one case - python -m tests.golden.run --engine P2 one engine - python -m tests.golden.run --vti also write one .vti per run (open in ParaView / pyvista) +python -m tests.golden.run run everything, write figures + report.html to out/ +python -m tests.golden.run --update also (re)write the golden files from the current code +python -m tests.golden.run --case dome one case +python -m tests.golden.run --engine P2 one engine +python -m tests.golden.run --vti also write one .vti per run (open in ParaView / pyvista) """ + from __future__ import annotations import argparse import html -import json import numpy as np from . import geomodels, harness, render -METRIC_COLS = ["value_rms", "angle_median", "angle_p90", "flipped_fraction", "data_value_rms", "data_normal_angle_max"] +METRIC_COLS = [ + "value_rms", + "angle_median", + "angle_p90", + "flipped_fraction", + "data_value_rms", + "data_normal_angle_max", +] def write_vti(case_name: str, engine: str, result: dict, n: int = 33): @@ -52,7 +59,9 @@ def report(all_results: dict, golden_status: dict): st = golden_status.get((case_name, eng), "") cells = "".join(f"{r['metrics'][c]:.4g}" for c in METRIC_COLS) body += f"{eng}{r['info']['dof']}{cells}{st}" - rows.append(f"{head}{body}
enginedofvs golden
") + rows.append( + f"{head}{body}
enginedofvs golden
" + ) page = ( "Interpolation golden report" "