From 4cc492fed436f922fbb279fc4e4e1645354babb8 Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Fri, 28 Aug 2026 11:40:41 +0200 Subject: [PATCH 1/4] [ENH] Add tests and example model for finite-fault handling with diagnostics Introduces a test case and example model (`last_model.gempy`) for finite-fault scalar field handling and visualization diagnostics. Includes validation of finite-fault descriptions, boundary conditions, and ellipsoid distance functions. Adds diagnostic plots for fault layer meshes and finite-fault scalar distribution. --- tests/test_common/test_core/last_model.gempy | Bin 0 -> 8037 bytes .../test_core/test_last_model_finite_fault.py | 208 ++++++++++++++++++ 2 files changed, 208 insertions(+) create mode 100644 tests/test_common/test_core/last_model.gempy create mode 100644 tests/test_common/test_core/test_last_model_finite_fault.py diff --git a/tests/test_common/test_core/last_model.gempy b/tests/test_common/test_core/last_model.gempy new file mode 100644 index 0000000000000000000000000000000000000000..a3055223b38bf4aeba54cbfe5c41ba890ea87d06 GIT binary patch literal 8037 zcmeHMTd3qn8Sa@g>dx%C3p(ya9XXaGC~P-Lx|4JUnF;#lx^Pg@hhZtIlS+4a(n)V6 zIWuQ;U>+0{d=SJ3VPFLXUqunYmz@Ri!KY=JXV=Ab!ChF^3wU|azbcomB$b}aBKR;F z&ZJZI)n9-8_1Ax^gJ zngG;rS({!!rT?|=}Ol|^9G`1{fI)~Y;)&|$q>vyF%^=J2u zgf`L`=bR*2l=3xDZj`PI7T-JQp_DUF5UJNfOLj}`frxeDy_1NufW(C3h|xI2B*&o5 zRBCp$Z^Jr-K%D_p>!(&}uJ>A>>l0LM5T~)+ZXm2tUsT%88YFaL5S77FHbpKBY1?EGASa6uoObBBa4~p{v z9wBz-mCcO#I{Y_n9fXZ?mC7ycLo{IXi-f` zXcIf6n3Iqd!@i}4o6f|tcJt*-yxDzm60*%*voW1`cMvCV>fE11ae>$l-hi zyF43DON)y`x&%!zCdROMP_~I#Lcwz;IZxxQHykwieT$%R04Pi2LWDt;PPe4r!S^+} z2@XILV6jV1HgHBn#PvT6L`vVgLsIt_hp-NAdRxY#mO1WXvt;Fxm_s2X^?YDaj!?(08SN8(vN4R8qcq9Yxoglm0ggt z(5pGbOt8HE3|jALYFo99_Xyp|W-lje>ovl{K;S!wOd~N}Q=5_IOVvM!FA|uq+PL2|W#t%Mpyiy=#h1?{n77w+OT zC=2Q-d?&&}+vWA4T3MDQ4!J79nkhlep`qz!b)2!FN{Lz`LuHnvrt8Y-+N2q);;^=E zkL?LmNQ}AbI;Q0rmbaCZI1&*0d$YMSoq4unSnky6=m-p-xl?cM8OGF{o0jP|%Tf*( zszXuWxIKOEtmqoFB!!z1>0OVWW0+G|jpI#?nP=E;dqF_|o&nnVv9PGL%wtudMZ{?W97eT1ZqR@@T^p zC@2=aJf_Bh(T7xC!(foEp(Y^nl~5`+E@ua5+r3Hg2)V7x*1hSJ*9wL22)4!{cvMKgW<0X%-22~46gQ*T>9$w{eUA!kvJ>GK7P*+L{kGNsZi zRfTi9ETAUDuGM(Z0;9Xt+3^}&KZm@kMP0-Lt5np5P?9E82_+;~gQLh*47Kr+LutXf zCUIzFp2jdN353Xi&n}xuv0PH9YOObAvAJFnil7i)rmL{GN3dE+qq94u_(DdONt!|0 zEI89`dx--!4uNFM0^ya=6~dG5`V}A-y9qeW-jygZUZ|c5Y)p#kNU11a$%PG{^h(d` zXsK~xU47Kza8>LEZoBv@`U(nuMBo;I;(P@rV3o#VX=Af@*eg8cJcY{-1v>~`vjj_a zfndjrz+FqW!VrCN;s1tIy90s>e0C`0y0waTAPYA(iQH1nxlv>vtP0rwj9n9|XdB|` zs8ely%2WqX8iH*McyuUZmBFP>@xaP0(o9S!d$sBj(nTQ(Ii2xHIv824hfc-dDjhLy z?sAm!;n3aHB2EJkKvlZQaDKcIoGC{VWo)@l;T9x|Q`uLFsf_lopSga}x@&v+(eR0b z6Qhy%I}LYjGE$j9O5uO;OOL(r<0t$Zzp_t`Mo%{HU%&A79^h{L@}I9?9KCb;k+1yz z!@mFh18CcC-YUlb#$&&G_FwP!M|k38!8mw)<;Pddo=U--i9$A8#he)Nyu z{_Nr{^nurZ^1Y{g|H(59pWpXyee35pfAi*tYB~}FbHa1hm#@ydxcSzlvR(e~zxebY zUaZ#Qzjp1Lmqvej=r6bFIqg_wT=nyvOVt|X*tXbJ^h59=0(2quHtMC<{{8t+!N!Qc z#{m^Ha55pLgy$=n6smWknW!iR(44Jgd~mr6^3;S|bHBKaWTf+usUKQ6-eM6s(Rk)A zyzz`0!I(@9GM<}(8!aNJ0?fUx*=^mtLIpHxMO1>gVo$ZRQrj~Z#9fT-$Ti2ag>8@L zUSN-1XmlgTa><*fS?ht)9g r7JFVlH%|QVqpAk9x$Nq|cGT%}LmF}L%zNG|aHrtkoA9O|6JP%Y*6kng literal 0 HcmV?d00001 diff --git a/tests/test_common/test_core/test_last_model_finite_fault.py b/tests/test_common/test_core/test_last_model_finite_fault.py new file mode 100644 index 00000000..5622a6f8 --- /dev/null +++ b/tests/test_common/test_core/test_last_model_finite_fault.py @@ -0,0 +1,208 @@ +import copy +import os +from pathlib import Path + +import numpy as np +import pytest + +from gempy_engine.core.data import TaperType +from gempy_engine.modules.faults.finite_faults import get_ellipsoid_distance, get_local_frame + + +MODEL_PATH = Path(__file__).with_name("last_model.gempy") + + +def _load_model(): + gempy = pytest.importorskip("gempy", reason="Loading .gempy files currently requires the GemPy package") + with pytest.warns(UserWarning, match="still in development"): + return gempy, gempy.load_model(str(MODEL_PATH)) + + +def _directional_radius(radius, positive: bool) -> float: + if isinstance(radius, tuple): + return radius[0 if positive else 1] + return radius + + +def _plot_finite_fault_diagnostics( + model, + no_fault_model, + finite_fault, + normal, + finite_fault_scalar, + output_dir: Path, +) -> tuple[Path, ...]: + try: + import matplotlib + except ImportError: + return () + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + output_dir.mkdir(parents=True, exist_ok=True) + u, v, w = get_local_frame(normal, angle_deg=finite_fault.rotation_deg) + center = np.asarray(finite_fault.center) + + max_radius = max( + _directional_radius(finite_fault.strike_radius, True), + _directional_radius(finite_fault.strike_radius, False), + _directional_radius(finite_fault.dip_radius, True), + _directional_radius(finite_fault.dip_radius, False), + ) + local_range = np.linspace(-1.15 * max_radius, 1.15 * max_radius, 301) + strike, dip = np.meshgrid(local_range, local_range) + plane_points = center + strike[..., None] * u + dip[..., None] * v + distance = get_ellipsoid_distance( + points=plane_points.reshape(-1, 3), + center=center, + u=u, + v=v, + a=finite_fault.strike_radius, + b=finite_fault.dip_radius, + ).reshape(strike.shape) + grid_relative = model.grid.values - center + grid_strike = grid_relative @ u + grid_dip = grid_relative @ v + grid_normal = np.abs(grid_relative @ w) + section = np.argsort(grid_normal)[:max(1024, len(grid_normal) // 32)] + fault_element = model.structural_frame.get_group_by_name("fault_series").get_element_by_name("fault") + surface_relative = fault_element.surface_points.xyz - center + + footprint_path = output_dir / "last_model_computed_finite_fault_taper.png" + fig, ax = plt.subplots(figsize=(9, 7), constrained_layout=True) + image = ax.scatter( + grid_strike[section], + grid_dip[section], + c=finite_fault_scalar[section], + s=20, + cmap="viridis", + vmin=0.0, + vmax=1.0, + ) + ax.contour(strike, dip, distance, levels=[1.0], colors="red", linewidths=2.0) + ax.scatter(surface_relative @ u, surface_relative @ v, marker="x", s=80, color="cyan", label="Fault points") + ax.scatter(0.0, 0.0, marker="+", s=120, color="white", label="Finite-fault center") + ax.set( + title=f"{model.meta.name}: computed finite-fault taper near the fault plane", + xlabel="Local strike coordinate u", + ylabel="Local dip coordinate v", + aspect="equal", + ) + ax.legend(loc="upper right") + fig.colorbar(image, ax=ax, label="Slip multiplier") + fig.savefig(footprint_path, dpi=160) + plt.close(fig) + + finite_layer = model.structural_frame.get_group_by_name("stratigraphic_series").get_element_by_name("layer") + no_fault_layer = no_fault_model.structural_frame.get_group_by_name("stratigraphic_series").get_element_by_name("layer") + finite_relative = finite_layer.vertices - center + no_fault_relative = no_fault_layer.vertices - center + + mesh_path = output_dir / "last_model_layer_mesh_comparison.png" + fig, ax = plt.subplots(figsize=(9, 7), constrained_layout=True) + ax.contour(strike, dip, distance, levels=[1.0], colors="red", linewidths=2.0) + ax.scatter( + no_fault_relative @ u, + no_fault_relative @ v, + s=8, + color="black", + alpha=0.35, + label="No-fault layer mesh", + ) + ax.scatter( + finite_relative @ u, + finite_relative @ v, + s=8, + color="tab:orange", + alpha=0.5, + label="Finite-fault layer mesh", + ) + ax.set( + title="Layer meshes projected onto the finite-fault plane", + xlabel="Local strike coordinate u", + ylabel="Local dip coordinate v", + aspect="equal", + ) + ax.legend() + fig.savefig(mesh_path, dpi=160) + plt.close(fig) + return footprint_path, mesh_path + + +def test_last_model_finite_fault_description_and_boundary(tmp_path, monkeypatch): + monkeypatch.setenv("SET_RAW_SCALAR_FIELDS_IN_SOLUTION", "True") + gempy, model = _load_model() + no_fault_model = copy.deepcopy(model) + fault_group = model.structural_frame.get_group_by_name("fault_series") + + assert fault_group.fault_type is gempy.data.FaultType.FINITE + assert fault_group.finite_fault_draft is None + finite_fault = fault_group.faults_input_data.finite_fault + assert finite_fault.center == pytest.approx((8.73, -4.2, 3.5515034198760986)) + assert finite_fault.strike_radius == pytest.approx((16.22, 6.94)) + assert finite_fault.dip_radius == pytest.approx((13.93, 16.66)) + assert finite_fault.taper is TaperType.QUADRATIC + assert finite_fault.rotation_deg == 0.0 + + from gempy.modules.data_manipulation import input_data_descriptor_from_geo_model + + engine_descriptor = input_data_descriptor_from_geo_model(model) + engine_finite_fault = engine_descriptor.stack_structure.faults_input_data[0].finite_fault + assert engine_finite_fault.center == pytest.approx((0.19135310, -0.55017848, -0.11561399)) + assert engine_finite_fault.strike_radius == pytest.approx((0.97543458, 0.41735610)) + assert engine_finite_fault.dip_radius == pytest.approx((0.83771909, 1.00189520)) + assert fault_group.faults_input_data.finite_fault is finite_fault + + fault_element = fault_group.get_element_by_name("fault") + normal = fault_element.orientations.grads[0] + u, v, w = get_local_frame(normal, angle_deg=finite_fault.rotation_deg) + assert np.allclose(np.stack((u, v, w)) @ np.stack((u, v, w)).T, np.eye(3)) + + center = np.asarray(finite_fault.center) + directions = ( + (u, _directional_radius(finite_fault.strike_radius, True)), + (-u, _directional_radius(finite_fault.strike_radius, False)), + (v, _directional_radius(finite_fault.dip_radius, True)), + (-v, _directional_radius(finite_fault.dip_radius, False)), + ) + for direction, radius in directions: + points = center + np.array([0.5, 1.0, 1.01])[:, None] * radius * direction + slip = finite_fault.calculate_slip(points, normal) + assert slip[0] == pytest.approx((1.0 - 0.5 ** 2) ** 2) + assert slip[1] == pytest.approx(0.0, abs=1e-28) + assert slip[2] == 0.0 + + grid_distance = get_ellipsoid_distance( + points=model.grid.values, + center=center, + u=u, + v=v, + a=finite_fault.strike_radius, + b=finite_fault.dip_radius, + ) + grid_slip = finite_fault.calculate_slip(model.grid.values, normal) + assert np.any(grid_distance < 1.0) + assert np.any(grid_distance >= 1.0) + assert np.all(grid_slip[grid_distance >= 1.0] == 0.0) + + no_fault_model.structural_frame.get_group_by_name("fault_series").fault_relations = ( + gempy.data.FaultsRelationSpecialCase.OFFSET_NONE + ) + gempy.compute_model(no_fault_model) + gempy.compute_model(model) + computed_finite_fault_scalar = model.solutions.raw_arrays.finite_fault_scalar_field_matrix[0] + assert computed_finite_fault_scalar.shape == (len(model.grid.values),) + assert np.any(computed_finite_fault_scalar > 0.0) + assert np.any(computed_finite_fault_scalar == 0.0) + assert np.all(computed_finite_fault_scalar >= 0.0) + + output_dir = Path(os.environ.get("FINITE_FAULT_PLOT_DIR", tmp_path)) + plot_paths = _plot_finite_fault_diagnostics( + model=model, + no_fault_model=no_fault_model, + finite_fault=finite_fault, + normal=normal, + finite_fault_scalar=computed_finite_fault_scalar, + output_dir=output_dir, + ) + assert all(path.is_file() and path.stat().st_size > 0 for path in plot_paths) From 9e4921f10fe7f51a72d8b91398b8527404d36ae5 Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Fri, 28 Aug 2026 11:40:51 +0200 Subject: [PATCH 2/4] Add tests and documentation for finite fault transformations Introduce tests ensuring finite fault transformations are correctly applied without mutating the model and validate isotropic transform enforcement. Update documentation to clarify the handling of finite fault coordinates and transformation constraints in GemPy Engine. --- docs/implementation/finite_faults.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/implementation/finite_faults.md b/docs/implementation/finite_faults.md index eed809f5..2ef3a29f 100644 --- a/docs/implementation/finite_faults.md +++ b/docs/implementation/finite_faults.md @@ -73,6 +73,12 @@ default profile. Supplying spline points for another taper is invalid. projects points onto the fault surface, so normal distance is not part of its two-dimensional footprint. A volumetric ellipsoid would be a separate model. +Direct GemPy Engine callers provide the center and radii in engine coordinates. +GemPy models persist these values in world coordinates and transform a runtime +copy together with the model inputs before calling the engine. The persisted +finite-fault definition is not modified. The current strike/dip representation +requires an isotropic transform that does not tilt the vertical axis. + ## Serialization Pydantic dataclasses use `TypeAdapter` for serialization and deserialization: From 075aad5aa0e476bad725125397b36c521fa730a1 Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Fri, 28 Aug 2026 15:18:39 +0200 Subject: [PATCH 3/4] [ENH] Add micro-point front-end and serialization design documentation Documents the proposed micro-point front-end data model, serialization schema, and local RBF correction methods. Includes the architectural pipeline, validation requirements, serialization updates, and tests for integration, ensuring micro-point observations enhance local compliance without affecting macro cokriging systems. --- MICRO_POINT_FRONTEND_DESIGN.md | 695 +++++++++++++++++++++++++++++++++ 1 file changed, 695 insertions(+) create mode 100644 MICRO_POINT_FRONTEND_DESIGN.md diff --git a/MICRO_POINT_FRONTEND_DESIGN.md b/MICRO_POINT_FRONTEND_DESIGN.md new file mode 100644 index 00000000..ca7aac23 --- /dev/null +++ b/MICRO_POINT_FRONTEND_DESIGN.md @@ -0,0 +1,695 @@ +# Micro Point Front-End and Serialization Design + +Status: proposed + +This document defines how high-density micro contacts should enter GemPy, cross +the GemPy-to-engine boundary, and be serialized. It also defines the local RBF +correction that consumes those contacts. + +Where this document conflicts with `MICRO_ANISOTROPIC_FIELD_DEFORMATION.md` or +`PLAN.md`, this document is authoritative. Those documents describe the +prototype and contain assumptions that are not valid for spatially varying +anisotropy. + +## 1. Purpose + +Micro points are dense observations used to make an existing macro geological +model locally comply with contacts without adding every contact to the global +cokriging system. + +The intended pipeline is: + +```text +GemPy macro observations + -> macro cokriging solve + -> macro scalar field and gradients + -> per-stack micro residual solve + -> macro field + local micro correction + -> activation, octree refinement, and mesh extraction +``` + +The macro model remains the structural hypothesis. Micro points are a local +compliance layer and do not replace surface points or orientations. + +## 2. Goals + +- Represent each micro point as a position, local geological frame, and local + anisotropic support. +- Use the same representation in a Python front end and a 3D editor. +- Associate every micro point with exactly one `StructuralElement`. +- Keep authored observations separate from options and solved runtime state. +- Preserve micro points through `gempy.save_model()` and `gempy.load_model()`. +- Apply a correction only to the stack and interface to which a point belongs. +- Use one mathematically consistent operator for fitting and evaluation. +- Preserve coordinate, dtype, device, and gradient consistency. + +## 3. Non-Goals + +- Micro points do not participate in the macro cokriging matrix. +- Solved residuals and RBF weights are not durable model input. +- Version 1 does not support one micro point constraining multiple interfaces. +- Version 1 does not define micro corrections across fault blocks. +- Version 1 does not accept perspective transforms or arbitrary projective + matrices. + +## 4. Terminology + +- **Macro point:** A standard `SurfacePointsTable` or `OrientationsTable` + observation used by the main interpolation system. +- **Micro point:** A dense contact observation used by the additive local + correction. +- **Support transform:** A `4x4` affine transform mapping normalized local + support coordinates to the model's world/input coordinate system. +- **Support scale:** The correlation lengths encoded in the linear part of a + support transform. It is not merely a display-gizmo scale. +- **Local RBF:** A radial basis function centered at one micro point and + evaluated using that point's support transform. + +## 5. Front-End Data Model + +### 5.1 Ownership + +`MicroPointsTable` should be owned by `StructuralElement`, in the same vein as +`SurfacePointsTable` and `OrientationsTable`: + +```python +class StructuralElement: + surface_points: SurfacePointsTable + orientations: OrientationsTable + micro_points: MicroPointsTable +``` + +This establishes the interface association without a separate mutable +model-level relation. Moving an element between groups moves its micro points; +removing an element removes them. Basement elements must have an empty micro +table. + +`StructuralFrame` should provide derived aggregate views: + +```text +micro_points_copy +number_of_micro_points_per_element +number_of_micro_points_per_group +``` + +The aggregate table is needed for binary serialization and engine conversion, +but it is not the authoritative mutable owner. + +The element key used in flattened binary rows must be stable across renaming. +The current fallback `StructuralElement.id` is derived from the element name +when `_id == -1`; that is not a sufficient durable foreign key. Before micro +points are persisted, the implementation must either materialize and serialize +an explicit element ID or introduce a stable element UUID. Dense surface and +stack indices remain runtime values derived from current structural order. + +### 5.2 Canonical Transform + +For point `i`, define: + +```text +x_world_h = H_world_from_support[i] @ x_support_h +``` + +with: + +```text +H_world_from_support = [ B_i p_i ] + [ 0 1 ] +``` + +- `p_i` is the micro-point position. +- The columns of `B_i` are local support axes expressed in world coordinates. +- The lengths of those columns are the kernel correlation lengths. +- The third local axis is the interface-normal direction. + +For the initial axisymmetric model: + +```text +B_i = R_i @ diag(lateral_range_i, lateral_range_i, normal_range_i) +``` + +where `normal_range_i < lateral_range_i` in the common case. Rotation around +the normal has no effect when both lateral ranges are equal, but retaining a +complete frame is convenient for 3D front ends and allows future triaxial +support. + +The position must not also be serialized as independent `X`, `Y`, and `Z` +fields. The translation column is authoritative. A convenience `xyz` property +may return `support_transforms[:, :3, 3]`. + +### 5.3 Scale Semantics + +The three support scales are physical correlation lengths in model coordinate +units. Applying the inverse transform produces dimensionless local coordinates. + +This removes the ambiguous double scaling in the prototype, where +`anisotropy_matrices` contain inverse ranges and the result is divided by a +second `kernel_range`. The durable representation has one source of geometric +range: the support transform. + +An optional global support multiplier may exist as an algorithm option, but it +must multiply all support lengths explicitly and must be applied identically +during fitting and evaluation. + +### 5.4 Proposed Table + +The conceptual public object is: + +```python +@dataclass +class MicroPointsTable: + data: np.ndarray + name_id_map: dict[str, int] | None = None + + @classmethod + def from_transforms( + cls, + support_transforms: np.ndarray, + names: Sequence[str] | str, + nugget: np.ndarray | None = None, + name_id_map: dict[str, int] | None = None, + ) -> "MicroPointsTable": ... + + @classmethod + def initialize_empty(cls) -> "MicroPointsTable": ... + + @property + def support_transforms(self) -> np.ndarray: ... + + @property + def xyz(self) -> np.ndarray: ... +``` + +The proposed version 1 structured dtype is: + +```python +np.dtype([ + ("support_transform", "", + "byte_order": "little" + } +} +``` + +A missing manifest means legacy version 1. + +Version 2 adds a dedicated member: + +```text +model.gempy +|-- header.json +|-- input.bin +|-- micro_points.bin +|-- grid.bin +`-- liquid_earth_meta.json +``` + +A dedicated member is preferable to appending data to `input.bin` because: + +- Existing readers currently ignore trailing `input.bin` bytes. +- A distinct member has an independently validated length and schema. +- Legacy surface-point and orientation layout remains unchanged. +- Future micro-table versions can evolve without changing macro table offsets. + +The structural-frame metadata should include: + +```json +{ + "micro_points": { + "dtype_version": 1, + "row_count": 42, + "byte_length": 6048 + } +} +``` + +The reader must use the fixed dtype selected by `dtype_version`; it must not +execute or blindly trust an arbitrary dtype supplied by the file. + +### 10.3 Save Flow + +1. Each `StructuralElement.micro_points.data` remains excluded from JSON. +2. `StructuralFrame.micro_points_copy` concatenates rows in structural order. +3. `model_to_bytes()` writes those bytes to `micro_points.bin`. +4. The ZIP member order and timestamps remain deterministic. +5. Serialization validation compares the original and loaded micro tables + directly, not through process-local `hash(bytes)` values. + +The ZIP writer should set `ZipInfo.compress_type` explicitly if compression is +expected. The current `make_info()` path creates stored members despite the +`ZipFile` compression setting. + +### 10.4 Load Flow + +1. Read and validate the serialization manifest. +2. Read `micro_points.bin` for format version 2. +3. Verify its exact byte length and row-size divisibility before `np.frombuffer`. +4. Inject it through the binary loading context with `input.bin` and `grid.bin`. +5. Construct `StructuralElement` objects with empty micro tables by default. +6. Decode the global micro table using its fixed little-endian dtype. +7. Reject rows whose element IDs are unknown or duplicated ambiguously. +8. Redistribute rows to elements by `element_id`. +9. Run normal table and affine-transform validation. + +Loading a version 1 model, or a transitional archive with no +`micro_points.bin`, produces empty micro tables. Old model files therefore +remain loadable. + +### 10.5 What Is Not Serialized + +Do not serialize: + +- Engine-coordinate transforms. +- Inverse `3x3` anisotropy operators. +- Macro scalar samples or gradients. +- Target scalar values. +- Residuals. +- RBF weights. +- Solver factorizations or condition estimates. + +These values depend on the current macro model and are recomputed. If solved +state is cached in the future, it must be a disposable cache keyed by a strong +fingerprint over all inputs and options, not authoritative model data. + +## 11. Required Tests + +The following tests are required when this design is implemented. + +### 11.1 `MicroPointsTable` + +- Empty initialization. +- Construction from one and multiple support transforms. +- Exact dtype names, byte order, offsets, and 144-byte row size. +- `xyz` extraction from the translation column. +- Selection by element name and ID. +- Copy and writable-view behavior. +- Rejection of incorrect array shapes. +- Rejection of non-affine last rows. +- Rejection of NaN and infinity. +- Rejection of singular or ill-conditioned support blocks. +- Rejection of zero or negative support scales. +- Rejection of shear in the version 1 authored format. +- Rejection of negative nuggets and mismatched array lengths. + +Suggested location: + +```text +gempy/test/test_core/test_micro_points.py +``` + +### 11.2 Structural Ownership + +- Each element owns an independent micro table. +- Flattening preserves structural order. +- Per-element and per-group counts are correct. +- Redistributing a flattened table restores exact ownership. +- Moving an element between groups moves its micro points. +- Removing an element cannot leave dangling micro rows. +- Basement contributes zero rows and rejects authored micro points. + +### 11.3 Serialization + +- Empty micro tables round-trip through `save_model` and `load_model`. +- Multiple elements with different row counts round-trip exactly. +- All 16 transform values, IDs, and nuggets compare exactly after loading. +- Binary rows shuffled before loading are redistributed by element ID. +- Existing version 1 fixtures load with empty micro tables. +- Missing required version 2 members raise a clear error. +- Unsupported major versions raise a clear error. +- Truncated rows and incorrect byte lengths are rejected. +- Unknown element IDs are rejected rather than silently discarded. +- Two saves of the same model produce identical archive bytes. +- Large micro arrays remain in binary and never enter `header.json`. +- Saving after compute does not persist residuals or weights. + +Suggested location: + +```text +gempy/test/test_modules/test_serialize_model.py +``` + +### 11.4 Coordinate Conversion + +- Identity conversion. +- Translation, rotation, and isotropic model scaling. +- Nonuniform model scaling. +- Grid rotation around a nonzero pivot. +- Combined grid and input transforms. +- Micro center conversion matches the normal point-conversion pipeline. +- Normalized support distance is invariant between world and engine frames. +- Source matrices are not mutated during conversion. + +### 11.5 Local RBF + +- One-point correction. +- Solve/evaluate round trip with different support transforms at every point. +- All supported kernels use the same type during solve and evaluation. +- Zero nugget reproduces residuals within solver tolerance. +- Nonzero nugget has documented smoothing behavior. +- `strength=0` returns the exact macro field. +- Analytic micro gradients match finite differences. +- Scalar output is identical whether gradients are requested or not. +- Duplicate and nearly duplicate points fail or regularize deterministically. +- NumPy and PyTorch preserve dtype, device, and numerical parity. + +### 11.6 Integration + +- Micro points only modify their associated stack. +- Points on two interfaces receive the correct target scalar values. +- Fault and unsupported stack types reject micro correction clearly. +- Macro-preservation policy has measurable, documented behavior. +- Corrected gradients reach dual contouring. +- Contact-driven octree refinement resolves isolated contacts. +- Extracted interfaces approach contacts within the configured tolerance. +- Plotting is opt-in and disabled in automated tests. + +## 12. Implementation Sequence + +The recommended implementation order is: + +1. Add and validate `MicroPointsTable` in the user-facing `gempy` package. +2. Attach an empty table to every `StructuralElement` and aggregate it through + `StructuralFrame`. +3. Introduce serialization version 2 and `micro_points.bin` with compatibility + tests for existing files. +4. Add GemPy APIs for adding, modifying, and deleting micro points. +5. Add the numerical micro data object to `InterpolationInput` and partition + metadata to `InputDataDescriptor`. +6. Compose support transforms into engine coordinates in the GemPy engine + factory. +7. Slice micro data per stack and derive residuals from immutable macro + interface values. +8. Replace the prototype solve with the consistent local RBF system. +9. Add micro gradients, backend parity, and conditioning diagnostics. +10. Add contact-driven octree refinement and mesh-level compliance tests. + +Each stage should leave models with no micro points behaviorally identical to +current models. + +## 13. Open Decisions + +- Whether macro preservation uses all macro surface points, one reference point + per interface, or nearby weighted anchors. +- Which nonsymmetric solver is used after the dense reference implementation. +- Whether per-point nugget is a direct diagonal regularizer or derived from a + separately named uncertainty measurement. +- Whether triaxial supports with unequal lateral scales are exposed in the + first public API or only accepted through complete support transforms. +- How micro correction is masked across finite-fault domains. +- Whether support transforms are editable through Euler/TRS convenience APIs + while retaining the raw matrix as the canonical representation. + +These decisions do not change the core contract: a micro observation is owned +by one structural element and persisted as a local-support-to-world `4x4` +transform whose scale defines anisotropic kernel support. From 7f44b564c671b106ed9e1cebf1516448b25e4c09 Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Fri, 4 Sep 2026 13:55:14 +0200 Subject: [PATCH 4/4] [DEL] --- MICRO_POINT_FRONTEND_DESIGN.md | 695 --------------------------------- 1 file changed, 695 deletions(-) delete mode 100644 MICRO_POINT_FRONTEND_DESIGN.md diff --git a/MICRO_POINT_FRONTEND_DESIGN.md b/MICRO_POINT_FRONTEND_DESIGN.md deleted file mode 100644 index ca7aac23..00000000 --- a/MICRO_POINT_FRONTEND_DESIGN.md +++ /dev/null @@ -1,695 +0,0 @@ -# Micro Point Front-End and Serialization Design - -Status: proposed - -This document defines how high-density micro contacts should enter GemPy, cross -the GemPy-to-engine boundary, and be serialized. It also defines the local RBF -correction that consumes those contacts. - -Where this document conflicts with `MICRO_ANISOTROPIC_FIELD_DEFORMATION.md` or -`PLAN.md`, this document is authoritative. Those documents describe the -prototype and contain assumptions that are not valid for spatially varying -anisotropy. - -## 1. Purpose - -Micro points are dense observations used to make an existing macro geological -model locally comply with contacts without adding every contact to the global -cokriging system. - -The intended pipeline is: - -```text -GemPy macro observations - -> macro cokriging solve - -> macro scalar field and gradients - -> per-stack micro residual solve - -> macro field + local micro correction - -> activation, octree refinement, and mesh extraction -``` - -The macro model remains the structural hypothesis. Micro points are a local -compliance layer and do not replace surface points or orientations. - -## 2. Goals - -- Represent each micro point as a position, local geological frame, and local - anisotropic support. -- Use the same representation in a Python front end and a 3D editor. -- Associate every micro point with exactly one `StructuralElement`. -- Keep authored observations separate from options and solved runtime state. -- Preserve micro points through `gempy.save_model()` and `gempy.load_model()`. -- Apply a correction only to the stack and interface to which a point belongs. -- Use one mathematically consistent operator for fitting and evaluation. -- Preserve coordinate, dtype, device, and gradient consistency. - -## 3. Non-Goals - -- Micro points do not participate in the macro cokriging matrix. -- Solved residuals and RBF weights are not durable model input. -- Version 1 does not support one micro point constraining multiple interfaces. -- Version 1 does not define micro corrections across fault blocks. -- Version 1 does not accept perspective transforms or arbitrary projective - matrices. - -## 4. Terminology - -- **Macro point:** A standard `SurfacePointsTable` or `OrientationsTable` - observation used by the main interpolation system. -- **Micro point:** A dense contact observation used by the additive local - correction. -- **Support transform:** A `4x4` affine transform mapping normalized local - support coordinates to the model's world/input coordinate system. -- **Support scale:** The correlation lengths encoded in the linear part of a - support transform. It is not merely a display-gizmo scale. -- **Local RBF:** A radial basis function centered at one micro point and - evaluated using that point's support transform. - -## 5. Front-End Data Model - -### 5.1 Ownership - -`MicroPointsTable` should be owned by `StructuralElement`, in the same vein as -`SurfacePointsTable` and `OrientationsTable`: - -```python -class StructuralElement: - surface_points: SurfacePointsTable - orientations: OrientationsTable - micro_points: MicroPointsTable -``` - -This establishes the interface association without a separate mutable -model-level relation. Moving an element between groups moves its micro points; -removing an element removes them. Basement elements must have an empty micro -table. - -`StructuralFrame` should provide derived aggregate views: - -```text -micro_points_copy -number_of_micro_points_per_element -number_of_micro_points_per_group -``` - -The aggregate table is needed for binary serialization and engine conversion, -but it is not the authoritative mutable owner. - -The element key used in flattened binary rows must be stable across renaming. -The current fallback `StructuralElement.id` is derived from the element name -when `_id == -1`; that is not a sufficient durable foreign key. Before micro -points are persisted, the implementation must either materialize and serialize -an explicit element ID or introduce a stable element UUID. Dense surface and -stack indices remain runtime values derived from current structural order. - -### 5.2 Canonical Transform - -For point `i`, define: - -```text -x_world_h = H_world_from_support[i] @ x_support_h -``` - -with: - -```text -H_world_from_support = [ B_i p_i ] - [ 0 1 ] -``` - -- `p_i` is the micro-point position. -- The columns of `B_i` are local support axes expressed in world coordinates. -- The lengths of those columns are the kernel correlation lengths. -- The third local axis is the interface-normal direction. - -For the initial axisymmetric model: - -```text -B_i = R_i @ diag(lateral_range_i, lateral_range_i, normal_range_i) -``` - -where `normal_range_i < lateral_range_i` in the common case. Rotation around -the normal has no effect when both lateral ranges are equal, but retaining a -complete frame is convenient for 3D front ends and allows future triaxial -support. - -The position must not also be serialized as independent `X`, `Y`, and `Z` -fields. The translation column is authoritative. A convenience `xyz` property -may return `support_transforms[:, :3, 3]`. - -### 5.3 Scale Semantics - -The three support scales are physical correlation lengths in model coordinate -units. Applying the inverse transform produces dimensionless local coordinates. - -This removes the ambiguous double scaling in the prototype, where -`anisotropy_matrices` contain inverse ranges and the result is divided by a -second `kernel_range`. The durable representation has one source of geometric -range: the support transform. - -An optional global support multiplier may exist as an algorithm option, but it -must multiply all support lengths explicitly and must be applied identically -during fitting and evaluation. - -### 5.4 Proposed Table - -The conceptual public object is: - -```python -@dataclass -class MicroPointsTable: - data: np.ndarray - name_id_map: dict[str, int] | None = None - - @classmethod - def from_transforms( - cls, - support_transforms: np.ndarray, - names: Sequence[str] | str, - nugget: np.ndarray | None = None, - name_id_map: dict[str, int] | None = None, - ) -> "MicroPointsTable": ... - - @classmethod - def initialize_empty(cls) -> "MicroPointsTable": ... - - @property - def support_transforms(self) -> np.ndarray: ... - - @property - def xyz(self) -> np.ndarray: ... -``` - -The proposed version 1 structured dtype is: - -```python -np.dtype([ - ("support_transform", "", - "byte_order": "little" - } -} -``` - -A missing manifest means legacy version 1. - -Version 2 adds a dedicated member: - -```text -model.gempy -|-- header.json -|-- input.bin -|-- micro_points.bin -|-- grid.bin -`-- liquid_earth_meta.json -``` - -A dedicated member is preferable to appending data to `input.bin` because: - -- Existing readers currently ignore trailing `input.bin` bytes. -- A distinct member has an independently validated length and schema. -- Legacy surface-point and orientation layout remains unchanged. -- Future micro-table versions can evolve without changing macro table offsets. - -The structural-frame metadata should include: - -```json -{ - "micro_points": { - "dtype_version": 1, - "row_count": 42, - "byte_length": 6048 - } -} -``` - -The reader must use the fixed dtype selected by `dtype_version`; it must not -execute or blindly trust an arbitrary dtype supplied by the file. - -### 10.3 Save Flow - -1. Each `StructuralElement.micro_points.data` remains excluded from JSON. -2. `StructuralFrame.micro_points_copy` concatenates rows in structural order. -3. `model_to_bytes()` writes those bytes to `micro_points.bin`. -4. The ZIP member order and timestamps remain deterministic. -5. Serialization validation compares the original and loaded micro tables - directly, not through process-local `hash(bytes)` values. - -The ZIP writer should set `ZipInfo.compress_type` explicitly if compression is -expected. The current `make_info()` path creates stored members despite the -`ZipFile` compression setting. - -### 10.4 Load Flow - -1. Read and validate the serialization manifest. -2. Read `micro_points.bin` for format version 2. -3. Verify its exact byte length and row-size divisibility before `np.frombuffer`. -4. Inject it through the binary loading context with `input.bin` and `grid.bin`. -5. Construct `StructuralElement` objects with empty micro tables by default. -6. Decode the global micro table using its fixed little-endian dtype. -7. Reject rows whose element IDs are unknown or duplicated ambiguously. -8. Redistribute rows to elements by `element_id`. -9. Run normal table and affine-transform validation. - -Loading a version 1 model, or a transitional archive with no -`micro_points.bin`, produces empty micro tables. Old model files therefore -remain loadable. - -### 10.5 What Is Not Serialized - -Do not serialize: - -- Engine-coordinate transforms. -- Inverse `3x3` anisotropy operators. -- Macro scalar samples or gradients. -- Target scalar values. -- Residuals. -- RBF weights. -- Solver factorizations or condition estimates. - -These values depend on the current macro model and are recomputed. If solved -state is cached in the future, it must be a disposable cache keyed by a strong -fingerprint over all inputs and options, not authoritative model data. - -## 11. Required Tests - -The following tests are required when this design is implemented. - -### 11.1 `MicroPointsTable` - -- Empty initialization. -- Construction from one and multiple support transforms. -- Exact dtype names, byte order, offsets, and 144-byte row size. -- `xyz` extraction from the translation column. -- Selection by element name and ID. -- Copy and writable-view behavior. -- Rejection of incorrect array shapes. -- Rejection of non-affine last rows. -- Rejection of NaN and infinity. -- Rejection of singular or ill-conditioned support blocks. -- Rejection of zero or negative support scales. -- Rejection of shear in the version 1 authored format. -- Rejection of negative nuggets and mismatched array lengths. - -Suggested location: - -```text -gempy/test/test_core/test_micro_points.py -``` - -### 11.2 Structural Ownership - -- Each element owns an independent micro table. -- Flattening preserves structural order. -- Per-element and per-group counts are correct. -- Redistributing a flattened table restores exact ownership. -- Moving an element between groups moves its micro points. -- Removing an element cannot leave dangling micro rows. -- Basement contributes zero rows and rejects authored micro points. - -### 11.3 Serialization - -- Empty micro tables round-trip through `save_model` and `load_model`. -- Multiple elements with different row counts round-trip exactly. -- All 16 transform values, IDs, and nuggets compare exactly after loading. -- Binary rows shuffled before loading are redistributed by element ID. -- Existing version 1 fixtures load with empty micro tables. -- Missing required version 2 members raise a clear error. -- Unsupported major versions raise a clear error. -- Truncated rows and incorrect byte lengths are rejected. -- Unknown element IDs are rejected rather than silently discarded. -- Two saves of the same model produce identical archive bytes. -- Large micro arrays remain in binary and never enter `header.json`. -- Saving after compute does not persist residuals or weights. - -Suggested location: - -```text -gempy/test/test_modules/test_serialize_model.py -``` - -### 11.4 Coordinate Conversion - -- Identity conversion. -- Translation, rotation, and isotropic model scaling. -- Nonuniform model scaling. -- Grid rotation around a nonzero pivot. -- Combined grid and input transforms. -- Micro center conversion matches the normal point-conversion pipeline. -- Normalized support distance is invariant between world and engine frames. -- Source matrices are not mutated during conversion. - -### 11.5 Local RBF - -- One-point correction. -- Solve/evaluate round trip with different support transforms at every point. -- All supported kernels use the same type during solve and evaluation. -- Zero nugget reproduces residuals within solver tolerance. -- Nonzero nugget has documented smoothing behavior. -- `strength=0` returns the exact macro field. -- Analytic micro gradients match finite differences. -- Scalar output is identical whether gradients are requested or not. -- Duplicate and nearly duplicate points fail or regularize deterministically. -- NumPy and PyTorch preserve dtype, device, and numerical parity. - -### 11.6 Integration - -- Micro points only modify their associated stack. -- Points on two interfaces receive the correct target scalar values. -- Fault and unsupported stack types reject micro correction clearly. -- Macro-preservation policy has measurable, documented behavior. -- Corrected gradients reach dual contouring. -- Contact-driven octree refinement resolves isolated contacts. -- Extracted interfaces approach contacts within the configured tolerance. -- Plotting is opt-in and disabled in automated tests. - -## 12. Implementation Sequence - -The recommended implementation order is: - -1. Add and validate `MicroPointsTable` in the user-facing `gempy` package. -2. Attach an empty table to every `StructuralElement` and aggregate it through - `StructuralFrame`. -3. Introduce serialization version 2 and `micro_points.bin` with compatibility - tests for existing files. -4. Add GemPy APIs for adding, modifying, and deleting micro points. -5. Add the numerical micro data object to `InterpolationInput` and partition - metadata to `InputDataDescriptor`. -6. Compose support transforms into engine coordinates in the GemPy engine - factory. -7. Slice micro data per stack and derive residuals from immutable macro - interface values. -8. Replace the prototype solve with the consistent local RBF system. -9. Add micro gradients, backend parity, and conditioning diagnostics. -10. Add contact-driven octree refinement and mesh-level compliance tests. - -Each stage should leave models with no micro points behaviorally identical to -current models. - -## 13. Open Decisions - -- Whether macro preservation uses all macro surface points, one reference point - per interface, or nearby weighted anchors. -- Which nonsymmetric solver is used after the dense reference implementation. -- Whether per-point nugget is a direct diagonal regularizer or derived from a - separately named uncertainty measurement. -- Whether triaxial supports with unequal lateral scales are exposed in the - first public API or only accepted through complete support transforms. -- How micro correction is masked across finite-fault domains. -- Whether support transforms are editable through Euler/TRS convenience APIs - while retaining the raw matrix as the canonical representation. - -These decisions do not change the core contract: a micro observation is owned -by one structural element and persisted as a local-support-to-world `4x4` -transform whose scale defines anisotropic kernel support.