diff --git a/MICRO_ANISOTROPIC_FIELD_DEFORMATION.md b/MICRO_ANISOTROPIC_FIELD_DEFORMATION.md new file mode 100644 index 00000000..539c97c1 --- /dev/null +++ b/MICRO_ANISOTROPIC_FIELD_DEFORMATION.md @@ -0,0 +1,531 @@ +# Micro Anisotropic Field Deformation Notes + +This document summarizes the current prototype and the recommended next steps for moving the micro-correction idea into GemPy Engine with minimal disruption to the existing macro interpolation path. + +## Goal + +Preserve the existing GemPy macro model as the structural/geological hypothesis, then apply an optional local scalar-field correction that improves high-density contact compliance without inserting all contacts into the global cokriging system. + +The target architecture is: + +```text +macro cokriging solve -> macro scalar field -> optional micro correction -> final scalar field +``` + +The macro solve should remain unchanged as much as possible: + +- Do not add dense borehole contacts to the global covariance system. +- Do not modify `SolverInput` for macro interpolation. +- Do not modify covariance assembly, universal drift, fault drift, or cokriging weights. +- Keep the micro layer as an optional evaluator-side additive correction. + +## Current Prototype + +The prototype currently exists in: + +- `gempy_engine/core/data/options/micro_anisotropic_options.py` +- `gempy_engine/modules/evaluator/micro_anisotropic_evaluator.py` +- `tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_evaluator.py` +- `tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_macro_integration.py` + +The prototype supports: + +- Pure NumPy micro correction evaluation. +- Symmetric micro covariance solve. +- 2D and 3D anisotropy matrix construction from macro gradients. +- Optional additive correction in both evaluator paths: + - `symbolic_evaluator.py` + - `generic_evaluator.py` +- A pytest integration demo with optional plotting via `GEMPY_PLOT_MICRO=1`. + +Run the visual test with: + +```bash +MPLBACKEND=QtAgg GEMPY_PLOT_MICRO=1 DEFAULT_BACKEND=numpy \ + /home/leguark/.venv/2025/bin/pytest \ + tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_macro_integration.py -s +``` + +Run the headless tests with: + +```bash +MPLBACKEND=Agg DEFAULT_BACKEND=numpy \ + /home/leguark/.venv/2025/bin/pytest \ + tests/test_common/test_modules/test_evaluator/ -v +``` + +At the time of writing, the evaluator test group passes: + +```text +13 passed +``` + +## Micro Field Formula + +The final field is currently: + +```text +V_final(x) = V_macro(x) + V_micro(x) +``` + +where: + +```text +V_micro(x) = sum_i w_i * exp(-||A_i (x - p_i)|| / range) +``` + +Definitions: + +- `p_i`: micro constraint point. +- `w_i`: micro weight solved from the micro system. +- `A_i`: local anisotropy transform built from the macro gradient at `p_i`. +- `range`: micro kernel range, intended to be small relative to macro kernel range. + +The current prototype uses `micro_kernel_range = 0.5` in the 2D demo. + +## Micro Covariance Solve + +Weights are solved from: + +```text +K w = residuals +``` + +with a symmetric anisotropic distance: + +```text +Dist^2(i, j) = (p_i - p_j)^T * M_ij * (p_i - p_j) +M_ij = (A_i^T A_i + A_j^T A_j) / 2 +K_ij = exp(-Dist(i, j) / range) +``` + +The symmetric distance is important because the micro covariance matrix must stay symmetric and numerically stable. + +The passive evaluation uses the cheaper one-sided distance: + +```text +||A_i (x - p_i)|| +``` + +This means exact round-trip equality is guaranteed for isotropic/identical transforms, but not necessarily for strongly varying anisotropy. This is expected from the current formulation. + +## Target Scalar Logic + +The main correction made during prototyping was how targets are computed. + +Incorrect prototype logic: + +```text +target_scalar = macro scalar at one arbitrary micro contact +residual_i = target_scalar - V_macro(contact_i) +``` + +Corrected logic: + +1. Evaluate macro scalar values at original macro surface points. +2. Split those values by surface using `TensorsStructure.number_of_points_per_surface`. +3. Compute one target scalar per surface/interface. +4. Assign every micro contact to a surface/interface id. +5. Compute per-contact residuals using that contact's assigned surface target. + +Current test uses median targets: + +```python +target_per_surface = [ + median(V_macro(surface_0_points)), + median(V_macro(surface_1_points)), +] + +target_values_at_contacts = target_per_surface[micro_surface_ids] +residuals = target_values_at_contacts - V_macro(micro_contacts) +``` + +Median was chosen because it is robust and simple. It is also suitable while the system is still experimental. + +## Macro Preservation Constraint + +The most useful refinement so far is adding original macro surface points as zero-residual constraints in the micro solve. + +Instead of solving only with micro contacts: + +```text +points = contacts +residuals = target_surface - V_macro(contact) +``` + +the current prototype solves with: + +```text +points = [contacts, macro_surface_points] +residuals = [target_surface - V_macro(contact), 0] +``` + +This tells the micro field: + +- Correct the field at micro contacts. +- Preserve the already-authored macro field at original macro surface points. + +In the 2D prototype this produced: + +```text +RMS before: 2.459 +RMS after: 0.311 +Macro point drift max: 0.699 +Macro point drift mean: 0.185 +``` + +This is a good first result: contact compliance improves strongly while macro points move much less than the largest contact residuals. + +## Anisotropy Construction + +Anisotropy matrices are built from the macro gradient sampled at each micro constraint point. + +The local frame is: + +- Last local axis: normalized macro gradient, treated as stratigraphic up / vertical. +- Remaining local axes: lateral directions. + +The scale matrix is: + +```text +2D: S = diag(1 / r_lateral, 1 / r_vertical) +3D: S = diag(1 / r_lateral, 1 / r_lateral, 1 / r_vertical) +``` + +The transform is: + +```text +A_i = S * R_i^T +``` + +where `R_i` contains the local basis vectors as columns. + +Interpretation: + +- Smaller `r_vertical` means faster decay across stratigraphy. +- Larger `r_lateral` means wider influence along the layer. + +The visual test draws ellipses for contact anisotropy using: + +```text +||A_i d|| = micro_kernel_range +``` + +In 3D, the equivalent visualization would be ellipsoids or principal axes arrows. + +## Minimal Production Architecture + +Keep the production change centered around an optional evaluator overlay. + +### Data Object + +Continue with an option object similar to: + +```python +class MicroAnisotropicOptions(BaseModel): + enabled: bool = False + points: Optional[np.ndarray] = None + residuals: Optional[np.ndarray] = None + anisotropy_matrices: Optional[np.ndarray] = None + weights: Optional[np.ndarray] = None + kernel_range: float = 1.0 + nugget: float = 0.0 +``` + +Potential additions: + +```python +strength: float = 1.0 +preserve_macro_points: bool = True +r_vertical: float = 0.5 +r_lateral: float = 5.0 +``` + +If `strength` is added, evaluation becomes: + +```text +V_final = V_macro + strength * V_micro +``` + +This is useful as a diagnostic/tuning knob, but it is not a replacement for macro zero constraints. + +### Evaluator Hook + +Keep the hook after macro scalar evaluation: + +```python +scalar_field = scalar_field + evaluate_micro_correction(...) +``` + +Do this in both: + +- `symbolic_evaluator` +- `generic_evaluator` + +Reason: PyKeOps currently has known LazyTensor compatibility issues in parts of the evaluator stack. The generic path must remain capable of exercising the micro workflow. + +### Avoid Touching + +Avoid touching these until necessary: + +- `compute_weights()` +- `_solve_interpolation()` +- covariance matrix assembly +- kernel constructor internals +- `SolverInput` semantics +- stack loop / octree loop + +## Moving To 3D + +The 3D implementation should follow the same steps as the 2D test, but with 3-component coordinates and gradients. + +### 3D Pipeline + +1. Run macro solve normally. +2. Evaluate macro scalar and gradient at micro contacts. +3. Evaluate macro scalar and gradient at macro surface points. +4. Compute target scalar per surface: + +```python +target_per_surface[s] = median(V_macro(points_of_surface_s)) +``` + +5. Assign each micro contact a surface id: + +```python +micro_surface_ids: np.ndarray # shape (N_contacts,) +``` + +6. Compute contact residuals: + +```python +contact_residuals = target_per_surface[micro_surface_ids] - V_macro(micro_contacts) +``` + +7. Build augmented constraints: + +```python +constraint_points = np.vstack([micro_contacts, macro_surface_points]) +constraint_residuals = np.concatenate([contact_residuals, zeros_for_macro_points]) +constraint_gradients = np.vstack([grad_at_contacts, grad_at_macro_points]) +``` + +8. Build anisotropy matrices: + +```python +A = compute_anisotropy_matrices_from_gradients( + constraint_points, + constraint_gradients, + r_vertical=..., + r_lateral=..., +) +``` + +9. Solve micro weights: + +```python +micro_weights = solve_micro_weights( + constraint_points, + constraint_residuals, + A, + kernel_range=..., + nugget=..., +) +``` + +10. Store on `options.evaluation_options.micro_anisotropic` and evaluate final field. + +### 3D Frame Construction + +For each gradient `g`: + +```text +z_axis = normalize(g) +ref = [0, 1, 0] +if abs(dot(z_axis, ref)) > 0.99: + ref = [1, 0, 0] +x_axis = normalize(cross(z_axis, ref)) +y_axis = normalize(cross(z_axis, x_axis)) +R = [x_axis, y_axis, z_axis] +A = S @ R.T +``` + +This exists in the prototype and should be kept unless a more geologically meaningful strike direction is available. + +### 3D Tests To Add + +Add a test using `simple_model` or another lightweight 3D fixture: + +- Macro solve. +- Choose synthetic 3D micro contacts assigned to one surface. +- Compute per-surface median target scalar. +- Add macro surface points as zero constraints. +- Solve micro correction. +- Evaluate at contacts and macro points. + +Assertions: + +```python +contact_rms_after < contact_rms_before +max_macro_point_drift < tolerance +np.all(np.isfinite(final_field)) +``` + +Start with a loose macro drift tolerance and tighten it after visual inspection. + +## Octree And Mesh Extraction + +The current prototype only corrects scalar values at evaluation locations. For mesh extraction to capture micro contacts reliably, the octree must eventually refine around micro contact locations. + +Without this, a micro contact can lie inside a coarse cell that never gets evaluated finely enough for dual contouring to capture the corrected crossing. + +Recommended staged approach: + +### Stage 1: Scalar Evaluation Only + +Current state. Validate math and field behavior. + +### Stage 2: Corners / Dense Grid Evaluation + +Evaluate micro correction on the grid or octree corners used for mesh extraction. + +No octree logic changes yet. + +### Stage 3: Forced Refinement Around Contacts + +Add an optional octree refinement criterion: + +```text +refine cell if it intersects macro isosurface OR contains/near a micro contact +``` + +Try to implement this as a narrow optional hook in octree refinement, not in the macro interpolation loop. + +Potential option fields: + +```python +micro.force_octree_refinement: bool = False +micro.refinement_radius: float = ... +``` + +### Stage 4: Dual Contouring Validation + +Once micro correction is evaluated at final corner locations, run dual contouring and verify that extracted surfaces move toward micro contacts. + +## PyKeOps / GPU Path + +The prototype uses dense NumPy for the micro solve. This is correct for math validation. + +Next GPU path should be added behind the same public functions, not by changing evaluator call sites. + +Suggested evolution: + +```python +solve_micro_weights(..., backend="numpy" | "pykeops") +evaluate_micro_correction(..., backend="numpy" | "pykeops") +``` + +or use `BackendTensor` dispatch internally. + +Be careful: the project currently has known PyKeOps LazyTensor issues around NumPy ufuncs and some matrix operations. Keep the NumPy implementation as the reference path. + +Known environment/test pitfalls: + +- Some PyKeOps tests try to write to `/home/miguel`, causing permission errors. +- Some LazyTensor expressions fail with standard NumPy ufuncs like `sqrt`/`exp`. +- Run new tests with `MPLBACKEND=Agg` unless plotting intentionally. + +## Suggested Helper To Promote Later + +Once the 3D test looks good, promote the repeated integration-test logic into a helper, probably under `modules/evaluator` or a small `modules/micro_correction` package. + +Candidate function: + +```python +def build_micro_anisotropic_constraints( + macro_solver_input: SolverInput, + macro_weights: np.ndarray, + options: InterpolationOptions, + micro_contacts: np.ndarray, + micro_surface_ids: np.ndarray, + macro_surface_points: np.ndarray, + n_points_per_surface: np.ndarray, + r_vertical: float, + r_lateral: float, + kernel_range: float, + nugget: float, +) -> MicroAnisotropicOptions: + ... +``` + +Responsibilities: + +- Evaluate macro scalar/gradient at contacts. +- Evaluate macro scalar/gradient at macro surface points. +- Compute per-surface target scalars. +- Build augmented constraints with zero macro residuals. +- Build anisotropy matrices. +- Solve micro weights. +- Return populated `MicroAnisotropicOptions`. + +This helper should not call `compute_model()` and should not mutate macro input data. + +## Open Design Questions + +1. **Surface assignment source** + + The current prototype uses explicit `micro_surface_ids`. Production code needs this information from borehole/contact metadata. + +2. **Macro preservation strength** + + Zero macro constraints help, but there can still be drift. If stricter preservation is required, consider: + + - Smaller kernel range. + - Larger macro constraint weight via repeated constraints or lower nugget. + - A weighted solve. + - Adding nearby orientation constraints later. + +3. **Contact vs macro weighting** + + The current system treats contact residuals and zero macro residuals equally. Production may need weights: + + ```text + high confidence borehole contacts vs expert macro control points + ``` + +4. **Per-surface micro solve** + + Solving each surface independently may reduce cross-interface interference. This is attractive for production. + +5. **Faults** + + The prototype ignores fault-specific behavior. A later version should evaluate whether micro corrections should be isolated by fault block or respect fault masks. + +6. **Gradient quality** + + Anisotropy relies on macro gradients. Need safeguards for near-zero gradients, noisy gradients, and points outside the stable macro field region. + +## Recommended Next Steps + +1. Add a 3D pytest integration demo analogous to the current 2D test. +2. Add optional visual 3D diagnostics: + - Slices through macro/micro scalar fields. + - Micro contact residuals before/after. + - Principal axes of anisotropy ellipsoids. +3. Promote repeated constraint-building logic into a helper. +4. Add `strength` as a diagnostic/tuning parameter, but keep zero macro constraints as the preferred preservation mechanism. +5. Add per-surface or per-stack solve mode to reduce cross-talk. +6. Add optional contact-driven octree refinement. +7. Only after the NumPy path is stable, add PyKeOps versions of the micro solve/evaluation behind the same function signatures. + +## Current Best Mental Model + +Treat the macro model as the authored geological hypothesis and the micro field as a local, anisotropic, constrained residual corrector. + +The micro solve should answer: + +```text +What is the smallest local correction field that moves assigned contacts onto their intended surface scalar values while keeping macro control points unchanged? +``` + +That framing keeps the macro path intact and makes the micro layer optional, testable, and removable. diff --git a/PLAN.md b/PLAN.md new file mode 100644 index 00000000..ad272a58 --- /dev/null +++ b/PLAN.md @@ -0,0 +1,157 @@ +# Technical Implementation Report: Multiscale Anisotropic Field Deformation for Automated Data Compliance in Implicit Geological Modeling + +## Executive Summary + +Traditional implicit geological modeling frameworks face a fundamental trade-off between structural consistency and high-density data compliance. Incorporating thousands of borehole contact points directly into a global Universal Co-Kriging dual matrix introduces a cubic computational complexity bottleneck ($O(N^3)$) and risks severe structural distortions such as artificial dimples and broken fault mechanics. + +This report outlines a next-generation, hybrid **Multiscale Anisotropic Field Deformation** framework. By treating the geological model as an expert-driven structural hypothesis (the Macro model) and delegating high-density borehole snapping to a localized, GPU-accelerated geostatistical compliance layer (the Micro step), this architecture achieves exact data precision while preserving structural logic. The execution leverages **PyKeOps** for symbolic matrix operations and a **Conjugate Gradient (CG)** solver on the GPU, culminating in real-time execution ready for extraction via Dual Contouring and visualization in **Unity 6**. + +--- + +## 1. Architectural Philosophy: The Hypothesis-First Paradigm + +Rather than treating geological modeling as an automated black-box data-fitting problem, this framework decouples the modeling process into two distinct spatial frequencies: + +1. **The Macro-Scale Framework (Low Frequency):** The geologist uses a clean, sparse subset of structural data (regional dips, major fault geometries) to establish the primary architectural trend using GemPy. This acts as the regional geological hypothesis, ensuring tectonic and structural rules are strictly enforced. +2. **The Micro-Scale Compliance Layer (High Frequency):** An automated, local optimization pass adjusts the generated continuous scalar field to ensure the target isosurfaces intersect perfectly with thousands of borehole contacts without affecting the regional framework outside a localized "damage radius." + +--- + +## 2. Mathematical Formulation & Algorithmic Steps + +### Step 1: Macro Field Evaluation and Gradient Extraction + +The baseline GemPy model is evaluated to produce a continuous global structural trend field, $V_{macro}(\mathbf{x})$. At each of the $N$ borehole contact points $\mathbf{x}_i = (x_i, y_i, z_i)$, two operations are performed: + +1. **Trend Lookup:** The baseline scalar value is sampled: $V_{pred, i} = V_{macro}(\mathbf{x}_i)$. +2. **Gradient Sampling:** The normalized mathematical gradient vector is extracted: + +$$\mathbf{g}_i = \nabla V_{macro}(\mathbf{x}_i)$$ + + + +### Step 2: Local Horizon Anchoring and Residual Calculation + +Because GemPy operates on relative drift constraints (where the absolute scalar value of an interface emerges from gradient orientations rather than static inputs), residuals cannot be evaluated against an arbitrary global constant. + +Contacts are grouped by their respective geological horizons. For each horizon group, a localized anchor point $\mathbf{x}_{anchor}$ is selected. The target scalar value for that specific horizon is locked to the macro value at that anchor: $C_{target} = V_{macro}(\mathbf{x}_{anchor})$. Local scalar residuals are then computed for all points within the horizon group: + + +$$\Delta V_i = C_{target} - V_{macro}(\mathbf{x}_i)$$ + +### Step 3: Constructing the Geologically Aligned Anisotropy Tensors + +To prevent spherical "bullseye" artifacts and cross-layer data contamination, distances must be evaluated in a warped, local coordinate system. For each borehole point $i$, a localized anisotropic transformation matrix $\mathbf{A}_i$ is constructed using a reverse Translation-Rotation-Scale (TRS) process: + +$$\mathbf{A}_i = \mathbf{S} \cdot \mathbf{R}_i^T$$ + +* **Rotation ($\mathbf{R}_i^T$):** Aligns the world coordinate axes to the local geology using the sampled macro-gradient $\mathbf{g}_i$ as the local "Up" vector, while the strike and dip vectors establish the local "Right" and "Forward" planes. +* **Scale ($\mathbf{S}$):** Imposes a steep distance penalty perpendicular to the stratigraphy: + +$$\mathbf{S} = \begin{bmatrix} \frac{1}{r_{lateral}} & 0 & 0 \\ 0 & \frac{1}{r_{vertical}} & 0 \\ 0 & 0 & \frac{1}{r_{lateral}} \end{bmatrix}$$ + + + +Where $r_{vertical} \ll r_{lateral}$. The vertical range is strictly constrained to be smaller than the minimum stratigraphic thickness between adjacent horizons, mathematically isolating independent layers. + +--- + +## 3. High-Performance GPU Implementation via PyKeOps + +### Symmetric Distance Computation + +To guarantee that the global covariance matrix remains strictly positive semi-definite (preventing invalid mathematical spaces or solver divergence), a symmetric distance metric is applied between any two interacting data points $i$ and $j$: + +$$\text{Dist}^2(i, j) = (\mathbf{x}_i - \mathbf{x}_j)^T \left( \frac{\mathbf{A}_i^T\mathbf{A}_i + \mathbf{A}_j^T\mathbf{A}_j}{2} \right) (\mathbf{x}_i - \mathbf{x}_j)$$ + +### Symbolic Kernel and Optimization Loop + +The $N \times N$ covariance matrix $\mathbf{K}$ is initialized symbolically inside **PyKeOps** as a `LazyTensor`. This allows the GPU to compute covariance coefficients on the fly in registers, reducing memory consumption from $O(N^2)$ to $O(N)$ storage and entirely bypassing VRAM bottlenecks. + +```python +import torch +from pykeops.torch import LazyTensor + +X = torch.tensor(X_boreholes, dtype=torch.float32).cuda() +A = torch.tensor(A_matrices, dtype=torch.float32).cuda() + +# Map coordinates to localized transformation spaces +AX = torch.einsum('nij,nj->ni', A, X) + +x_i = LazyTensor(AX[:, None, :]) # (N, 1, 3) +x_j = LazyTensor(AX[None, :, :]) # (1, N, 3) + +dist_squared = ((x_i - x_j) ** 2).sum(-1) +K = (- (dist_squared.sqrt())).exp() # Exponential Covariance Kernel + +``` + +Because the macro model provides a highly accurate starting baseline, the residual vector $\mathbf{y} = [\Delta V_1, \dots, \Delta V_N]^T$ sits close to the final solution space. A PyKeOps-backed **Conjugate Gradient (CG)** solver solves the linear system $\mathbf{K}\mathbf{w} = \mathbf{y}$ for the weights vector $\mathbf{w}$ within a fraction of a second. + +--- + +## 4. Optimized Mesh Extraction via Dual-Criteria Octree Subdivision + +To minimize the computational overhead of the micro-evaluation pass, the final field deformation is restricted exclusively to the finest level of an adaptive octree structure, directly targeting the boundary cells where the Dual Contouring mesh is extracted. + +However, relying solely on the macro-geology field to guide octree subdivision introduces a critical geometric vulnerability: if a dense borehole contact point deviates significantly from the macro trend, it may fall outside the fine cells generated by the macro engine, landing in a massive, un-subdivided coarse block. Consequently, the micro-correction field at that location would be bypassed, and the extracted mesh would fail to capture the data point. + +To eliminate this data-omission risk, the framework utilizes a **Dual-Criteria Octree Subdivision** protocol. A spatial volume cell is forced to subdivide to its highest resolution tier if **either** of the following conditions is met: + +1. **The Macro Isosurface Criterion:** The cell intersects a target structural threshold of the baseline GemPy field ($V_{macro}(\mathbf{x}) = C_{target}$). +2. **The High-Density Data Criterion:** The cell boundaries enclose one or more raw borehole contact coordinates ($\mathbf{x}_i$). + +### Passive Metric Grid Evaluation + +When it is time to extract the mesh via Dual Contouring, the framework loops through the active cells at the finest resolution level. For each corner coordinate $\mathbf{x}_{corner}$ of a fine boundary cell, it looks up the baseline value and adds the fast, passive PyKeOps anisotropic distance lookup: + +$$V_{final}(\mathbf{x}_{corner}) = V_{macro}(\mathbf{x}_{corner}) + \sum_{i=1}^{N} w_i \cdot e^{-\|\mathbf{A}_i (\mathbf{x}_{corner} - \mathbf{x}_i)\|}$$ + +This eliminates the need to calculate new gradients at millions of grid nodes, keeping the final evaluation pass highly parallelized and computationally cheap. + +--- + +## 5. End-to-End Execution Pipeline + +By implementing this dual-gated subdivision rule, the downstream mesh extraction operates as a streamlined, highly parallelized graphics workflow: + +``` +[GemPy Macro Model] + [Borehole Coordinates] + │ + ▼ + [Dual-Criteria Octree Generation] + (Fine cells at macro horizons & well sites) + │ + ▼ + [Fine Cell Corner Evaluation] + (Passive Anisotropic PyKeOps Vector Lookup) + │ + ▼ + [Dual Contouring Extraction] + (Forced Isosurface Crossing at Well Caps) + │ + ▼ + [Unity 6 Mesh Buffers] + +``` + +1. **Sparse Global Architecture:** The octree remains broad and lightweight across the vast majority of the asset volume, preventing unnecessary scalar evaluation loops in homogeneous rock masses. +2. **Targeted Precision Hooks:** The grid is guaranteed to maintain ultra-high-resolution cell matrices immediately surrounding every well track, providing the necessary mathematical "hooks" for data snapping. +3. **Forced Crossings:** When the Dual Contouring engine processes the corners of these fine well-bounding cells, it samples the combined $V_{final}(\mathbf{x}_{corner})$ field. The micro-weights ($\mathbf{w}$) seamlessly shift the scalar values across the boundary threshold within that cell, forcing the extracted vertex to lock onto the physical borehole coordinate with millimeter precision before piping the clean topology directly to the **Unity 6** render buffers. + +--- + +## 6. Architectural Comparison Matrix + +| Property | Standard Global Engine | Pure Post-Process Mesh Snapping | Proposed Multiscale PyKeOps Framework | +| --- | --- | --- | --- | +| **Computational Complexity** | Cubic $O(N^3)$ — chokes on dense datasets. | Linear $O(N)$ — executed entirely on the client. | **Ultra-Fast $O(N)$** — GPU-accelerated symbolic linear iterations. | +| **Volumetric Field Consistency** | High. | Broken — visual mesh diverges from the underlying scalar volume. | **Perfect** — corrections are applied directly inside the 3D scalar volume. | +| **Structural Integrity (Faults/Dips)** | Maintained globally, but lacks local data compliance. | Poor — risks crossing surfaces, artifact generation, and fault smearing. | **Maintained** — locked to the macro-gradient; strong vertical anisotropy prevents layer cross-talk. | +| **Visual Artifacts** | None. | High risk of sharp conical dimples ("tents over poles"). | **None** — corrections stretch naturally into smooth geological ovals. | + +--- + +## Conclusion + +The **Multiscale Anisotropic Field Deformation** framework successfully bridges the gap between expert geological intuition and strict data compliance. By utilizing the macro-gradient of GemPy to dictate localized, anisotropic data transformations and leveraging PyKeOps for memory-efficient GPU parallelization, this implementation delivers an auditable, structurally sound, and lightning-fast modeling pipeline that satisfies both mathematical rigor and real-time interactive rendering demands. \ No newline at end of file diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py index e19e3f1c..67e3de64 100644 --- a/gempy_engine/core/data/options/evaluation_options.py +++ b/gempy_engine/core/data/options/evaluation_options.py @@ -1,10 +1,12 @@ import enum -from dataclasses import dataclass +from dataclasses import dataclass, field from typing import Annotated from typing_extensions import deprecated +from .micro_anisotropic_options import MicroAnisotropicOptions + class MeshExtractionMaskingOptions(enum.Enum): NOTHING = enum.auto() # * This is only for testing @@ -27,6 +29,8 @@ class EvaluationOptions: evaluation_chunk_size: int = 500_000 + micro_anisotropic: MicroAnisotropicOptions = field(default_factory=MicroAnisotropicOptions) + compute_scalar: bool = True compute_scalar_gradient: bool = False diff --git a/gempy_engine/core/data/options/micro_anisotropic_options.py b/gempy_engine/core/data/options/micro_anisotropic_options.py new file mode 100644 index 00000000..805e8bcf --- /dev/null +++ b/gempy_engine/core/data/options/micro_anisotropic_options.py @@ -0,0 +1,21 @@ +from typing import Literal, Optional + +import numpy as np +from pydantic import BaseModel, ConfigDict + +MicroKernelType = Literal["exponential", "matern_3_2", "matern_5_2"] + + +class MicroAnisotropicOptions(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + + enabled: bool = False + points: Optional[np.ndarray] = None # (N, 3) micro constraint points + residuals: Optional[np.ndarray] = None # (N,) target residual values + anisotropy_matrices: Optional[np.ndarray] = None # (N, 3, 3) per-point anisotropy transforms + weights: Optional[np.ndarray] = None # (N,) solved micro weights + kernel_range: float = 1.0 # range for the micro kernel + kernel_type: MicroKernelType = "matern_5_2" # kernel function for micro solve + eval + nugget: float = 0.0 # diagonal nugget for the micro solve + preserve_macro_points: bool = True # include macro SP as zero-residual constraints + strength: float = 1.0 # global strength multiplier (1.0 = full correction) diff --git a/gempy_engine/modules/evaluator/generic_evaluator.py b/gempy_engine/modules/evaluator/generic_evaluator.py index 49ba777c..c945a24f 100644 --- a/gempy_engine/modules/evaluator/generic_evaluator.py +++ b/gempy_engine/modules/evaluator/generic_evaluator.py @@ -2,6 +2,7 @@ import gc from typing import Optional +import gempy_engine.config from gempy_engine.core.backend_tensor import BackendTensor from gempy_engine.core.data import InterpolationOptions from gempy_engine.core.data.exported_fields import ExportedFields @@ -66,6 +67,21 @@ def generic_evaluator( if n_chunks > 5: print(f"Chunking done: {n_chunks} chunks") + micro = options.evaluation_options.micro_anisotropic + if micro.enabled and micro.weights is not None and micro.points is not None and micro.anisotropy_matrices is not None: + from gempy_engine.modules.evaluator.micro_anisotropic_evaluator import evaluate_micro_correction + if BackendTensor.engine_backend != gempy_engine.config.AvailableBackends.numpy: + scalar_field = BackendTensor.t.to_numpy(scalar_field) + correction = evaluate_micro_correction( + xyz_to_interpolate=solver_input.xyz_to_interpolate, + micro_points=micro.points, + micro_weights=micro.weights, + anisotropy_matrices=micro.anisotropy_matrices, + kernel_range=micro.kernel_range, + kernel_type=micro.kernel_type, + ) + scalar_field = scalar_field + correction + return ExportedFields(scalar_field, gx_field, gy_field, gz_field) diff --git a/gempy_engine/modules/evaluator/micro_anisotropic_evaluator.py b/gempy_engine/modules/evaluator/micro_anisotropic_evaluator.py new file mode 100644 index 00000000..80f63265 --- /dev/null +++ b/gempy_engine/modules/evaluator/micro_anisotropic_evaluator.py @@ -0,0 +1,189 @@ +import numpy as np +from typing import Literal, Optional + +MicroKernelType = Literal["exponential", "matern_3_2", "matern_5_2"] + + +def _kernel_value(r: np.ndarray, kernel_type: MicroKernelType) -> np.ndarray: + """Evaluate the micro radial kernel K(r) where r = anisotropic_distance / kernel_range. + + All kernels satisfy K(0) = 1 and are positive and finite for r >= 0. + + exponential — Matérn 1/2: K(r) = exp(-r) + matern_3_2 — Matérn 3/2: K(r) = (1 + sqrt(3) r) exp(-sqrt(3) r) + matern_5_2 — Matérn 5/2: K(r) = (1 + sqrt(5) r + 5r²/3) exp(-sqrt(5) r) + """ + if kernel_type == "exponential": + return np.exp(-r) + elif kernel_type == "matern_3_2": + a = np.sqrt(3.0) * r + return (1.0 + a) * np.exp(-a) + elif kernel_type == "matern_5_2": + a = np.sqrt(5.0) * r + return (1.0 + a + (5.0 / 3.0) * r * r) * np.exp(-a) + else: + raise ValueError(f"Unknown micro kernel type: {kernel_type}") + + +def evaluate_micro_correction( + xyz_to_interpolate: np.ndarray, # (M, 3) + micro_points: np.ndarray, # (N, 3) + micro_weights: np.ndarray, # (N,) + anisotropy_matrices: np.ndarray, # (N, 3, 3) + kernel_range: float = 1.0, + kernel_type: MicroKernelType = "exponential", +) -> np.ndarray: + """Evaluate the micro correction field at target points. + + V(x) = sum_i w_i * K(||A_i (x - p_i)|| / range) + """ + M = xyz_to_interpolate.shape[0] + N = micro_points.shape[0] + correction = np.zeros(M, dtype=np.float64) + + for j in range(N): + Aj = anisotropy_matrices[j] + wj = micro_weights[j] + pj = micro_points[j] + diffs = xyz_to_interpolate - pj[np.newaxis, :] + transformed = np.einsum('ij,mj->mi', Aj, diffs) + dists = np.linalg.norm(transformed, axis=1) + r = dists / kernel_range + correction += wj * _kernel_value(r, kernel_type) + + return correction + + +def build_micro_covariance( + micro_points: np.ndarray, # (N, 3) + anisotropy_matrices: np.ndarray, # (N, 3, 3) + kernel_range: float = 1.0, + kernel_type: MicroKernelType = "exponential", + nugget: float = 0.0, +) -> np.ndarray: + """Build the symmetric NxN covariance matrix for the micro solve. + + K[i,j] = K(||A_i (p_i - p_j)|| / range) + + where K is the selected micro kernel and distances use the symmetric + metric M_ij = (A_i^T A_i + A_j^T A_j) / 2. + """ + N = micro_points.shape[0] + K = np.zeros((N, N), dtype=np.float64) + + ATA = np.einsum('nki,nkj->nij', anisotropy_matrices, anisotropy_matrices) + + for i in range(N): + for j in range(i, N): + M_ij = 0.5 * (ATA[i] + ATA[j]) + diff = micro_points[i] - micro_points[j] + dist_sq = diff @ M_ij @ diff + dist = np.sqrt(max(dist_sq, 0.0)) + r = dist / kernel_range + val = float(_kernel_value(np.array(r), kernel_type)) + K[i, j] = val + K[j, i] = val + + if nugget > 0: + np.fill_diagonal(K, K.diagonal() + nugget) + + return K + + +def solve_micro_weights( + micro_points: np.ndarray, # (N, 3) + residuals: np.ndarray, # (N,) + anisotropy_matrices: np.ndarray, # (N, 3, 3) + kernel_range: float = 1.0, + kernel_type: MicroKernelType = "exponential", + nugget: float = 0.0, +) -> np.ndarray: + """Solve K @ w = residuals for the micro correction weights. + + Returns weights array of shape (N,). + """ + K = build_micro_covariance(micro_points, anisotropy_matrices, kernel_range, kernel_type, nugget) + weights = np.linalg.solve(K, residuals) + return weights + + +def compute_macro_values_at_micro_points( + xyz_to_interpolate: np.ndarray, + weights: np.ndarray, + solver_input: 'SolverInput', + options: 'InterpolationOptions', +) -> np.ndarray: + """Extract the macro scalar field at micro point locations. + + This evaluates the macro interpolation exactly at the micro contact points + to compute residuals = target_values - macro_values. + """ + from gempy_engine.modules.evaluator.symbolic_evaluator import symbolic_evaluator + from gempy_engine.core.data.internal_structs import SolverInput + + proxy_input = SolverInput( + sp_internal=solver_input.sp_internal, + ori_internal=solver_input.ori_internal, + xyz_to_interpolate=xyz_to_interpolate, + fault_internal=solver_input._fault_internal, + ) + + exported = symbolic_evaluator(proxy_input, weights, options) + return exported.scalar_field + + +def compute_anisotropy_matrices_from_gradients( + micro_points: np.ndarray, # (N, D) + gradients: np.ndarray, # (N, D) gradient vectors + r_vertical: float = 1.0, + r_lateral: float = 10.0, +) -> np.ndarray: + """Build per-point anisotropy matrices from macro gradient directions. + + A_i = S * R_i^T + + R_i^T projects world coordinates into a local frame aligned with the gradient + (last axis = gradient direction = stratigraphic up). + S = diag(lateral scale repeated, vertical scale) + + Works for 2D and 3D. + """ + N, D = micro_points.shape + assert gradients.shape == (N, D), f"gradients shape {gradients.shape} != (N, D) {(N, D)}" + + if D == 2: + lateral_scales = np.array([1.0 / r_lateral], dtype=np.float64) + scales = np.concatenate([lateral_scales, [1.0 / r_vertical]]) + S = np.diag(scales) + else: + S = np.diag(np.array([1.0 / r_lateral, 1.0 / r_lateral, 1.0 / r_vertical])) + + matrices = np.zeros((N, D, D), dtype=np.float64) + + for i in range(N): + grad = gradients[i].astype(np.float64) + grad_norm = np.linalg.norm(grad) + if grad_norm < 1e-10: + grad = np.zeros(D, dtype=np.float64) + grad[-1] = 1.0 + + z_axis = grad / np.linalg.norm(grad) + + if D == 2: + x_axis = np.array([z_axis[1], -z_axis[0]], dtype=np.float64) + R = np.column_stack([x_axis, z_axis]) + else: + ref = np.array([0.0, 1.0, 0.0], dtype=np.float64) + if abs(np.dot(z_axis, ref)) > 0.99: + ref = np.array([1.0, 0.0, 0.0], dtype=np.float64) + + x_axis = np.cross(z_axis, ref) + x_axis = x_axis / np.linalg.norm(x_axis) + y_axis = np.cross(z_axis, x_axis) + y_axis = y_axis / np.linalg.norm(y_axis) + + R = np.column_stack([x_axis, y_axis, z_axis]) + + matrices[i] = S @ R.T + + return matrices diff --git a/gempy_engine/modules/evaluator/symbolic_evaluator.py b/gempy_engine/modules/evaluator/symbolic_evaluator.py index 675f5718..a29a19e7 100644 --- a/gempy_engine/modules/evaluator/symbolic_evaluator.py +++ b/gempy_engine/modules/evaluator/symbolic_evaluator.py @@ -77,9 +77,51 @@ def symbolic_evaluator(solver_input: SolverInput, weights: np.ndarray, options: else: raise ValueError("Number of dimensions have to be 2 or 3") + scalar_field = _apply_micro_correction(scalar_field, solver_input, options) + return ExportedFields(scalar_field, gx_field, gy_field, gz_field) +def _apply_micro_correction(scalar_field: np.ndarray, solver_input: SolverInput, options: InterpolationOptions) -> np.ndarray: + micro = options.evaluation_options.micro_anisotropic + if not micro.enabled: + return scalar_field + if micro.weights is None or micro.points is None or micro.anisotropy_matrices is None: + return scalar_field + + from gempy_engine.modules.evaluator.micro_anisotropic_evaluator import evaluate_micro_correction + + correction = evaluate_micro_correction( + xyz_to_interpolate=solver_input.xyz_to_interpolate, + micro_points=micro.points, + micro_weights=micro.weights, + anisotropy_matrices=micro.anisotropy_matrices, + kernel_range=micro.kernel_range, + kernel_type=micro.kernel_type, + ) + return scalar_field + correction + + +def _apply_micro_correction_stacked(scalar_field: np.ndarray, eval_input: EvaluatorInput, options: InterpolationOptions) -> np.ndarray: + micro = options.evaluation_options.micro_anisotropic + if not micro.enabled: + return scalar_field + if micro.weights is None or micro.points is None or micro.anisotropy_matrices is None: + return scalar_field + + from gempy_engine.modules.evaluator.micro_anisotropic_evaluator import evaluate_micro_correction + + correction = evaluate_micro_correction( + xyz_to_interpolate=eval_input.xyz_to_interpolate, + micro_points=micro.points, + micro_weights=micro.weights, + anisotropy_matrices=micro.anisotropy_matrices, + kernel_range=micro.kernel_range, + kernel_type=micro.kernel_type, + ) + return scalar_field + correction + + def _build_block_sparse_ranges(M_sizes: list[int], N_sizes: list[int]): """Build PyKeOps block-sparse ranges tuple for block-diagonal evaluation.""" keep_i = np.cumsum([0] + N_sizes) @@ -271,13 +313,15 @@ def _run_prep(args): gx_field = gx_fields[idx] gy_field = gy_fields[idx] gz_field = gz_fields[idx] - + if BackendTensor.engine_backend == gempy_engine.config.AvailableBackends.numpy: s_field = BackendTensor.t.to_numpy(s_field) if gx_field is not None: gx_field = BackendTensor.t.to_numpy(gx_field) if gy_field is not None: gy_field = BackendTensor.t.to_numpy(gy_field) if gz_field is not None: gz_field = BackendTensor.t.to_numpy(gz_field) + s_field = _apply_micro_correction_stacked(s_field, eval_inputs[idx], options_list[idx]) + results.append(ExportedFields(s_field, gx_field, gy_field, gz_field)) return results diff --git a/tests/test_common/test_modules/test_evaluator/__init__.py b/tests/test_common/test_modules/test_evaluator/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_evaluator.py b/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_evaluator.py new file mode 100644 index 00000000..7b5e4e41 --- /dev/null +++ b/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_evaluator.py @@ -0,0 +1,212 @@ +import numpy as np + +from gempy_engine.modules.evaluator.micro_anisotropic_evaluator import ( + evaluate_micro_correction, + build_micro_covariance, + solve_micro_weights, + compute_anisotropy_matrices_from_gradients, +) + + +def _make_identity_anisotropy(N: int) -> np.ndarray: + return np.tile(np.eye(3, dtype=np.float64), (N, 1, 1)) + + +# ---------------------------------------------------------------- +# evaluate_micro_correction +# ---------------------------------------------------------------- + +def test_evaluate_single_point_identity(): + """One micro point at origin, identity anisotropy, weight=1.0, range=1.0.""" + xyz = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]], dtype=np.float64) + points = np.array([[0.0, 0.0, 0.0]], dtype=np.float64) + weights = np.array([1.0], dtype=np.float64) + A = _make_identity_anisotropy(1) + + correction = evaluate_micro_correction(xyz, points, weights, A, kernel_range=1.0) + + np.testing.assert_allclose(correction[0], 1.0, rtol=1e-10) + np.testing.assert_allclose(correction[1], np.exp(-1.0), rtol=1e-6) + np.testing.assert_allclose(correction[2], np.exp(-2.0), rtol=1e-6) + + +def test_evaluate_monotonic_decay(): + """Correction magnitude decreases monotonically with distance.""" + xyz = np.linspace(0, 5, 100)[:, np.newaxis] * np.array([1.0, 0.0, 0.0]) + points = np.array([[0.0, 0.0, 0.0]], dtype=np.float64) + weights = np.array([1.0], dtype=np.float64) + A = _make_identity_anisotropy(1) + + correction = evaluate_micro_correction(xyz, points, weights, A, kernel_range=1.0) + + diffs = np.diff(correction) + assert np.all(diffs <= 0), "Correction should decrease monotonically with distance" + + +def test_evaluate_vertical_anisotropy(): + """Vertical anisotropy: decay should be faster in Z than in XY.""" + range_val = 1.0 + r_vertical = 0.5 + r_lateral = 5.0 + + xyz = np.array([ + [1.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + ], dtype=np.float64) + points = np.array([[0.0, 0.0, 0.0]], dtype=np.float64) + weights = np.array([1.0], dtype=np.float64) + + grad = np.array([[0.0, 0.0, 1.0]], dtype=np.float64) + A = compute_anisotropy_matrices_from_gradients(points, grad, r_vertical, r_lateral) + + correction = evaluate_micro_correction(xyz, points, weights, A, kernel_range=range_val) + + xy_val = correction[0] + z_val = correction[1] + assert xy_val > z_val, ( + f"Lateral correction ({xy_val}) should be larger than vertical ({z_val}) " + f"since r_lateral > r_vertical" + ) + + +# ---------------------------------------------------------------- +# build_micro_covariance +# ---------------------------------------------------------------- + +def test_build_covariance_identity_symmetric_psd(): + """Covariance matrix with identity anisotropy is symmetric and PSD.""" + N = 5 + rng = np.random.default_rng(42) + points = rng.uniform(-5, 5, (N, 3)).astype(np.float64) + A = _make_identity_anisotropy(N) + + K = build_micro_covariance(points, A, kernel_range=2.0) + + np.testing.assert_allclose(K, K.T, atol=1e-14) + eigvals = np.linalg.eigvalsh(K) + assert np.all(eigvals >= -1e-10), f"K is not PSD: min eigenvalue = {eigvals.min()}" + + +def test_build_covariance_diagonal_max(): + """Diagonal entries are the largest in each row (kernel is maximum at zero distance).""" + N = 10 + rng = np.random.default_rng(123) + points = rng.uniform(-5, 5, (N, 3)).astype(np.float64) + A = _make_identity_anisotropy(N) + + K = build_micro_covariance(points, A, kernel_range=2.0) + + for i in range(N): + assert K[i, i] >= np.max(K[i, :]) - 1e-14, f"Row {i}: diag {K[i,i]:.6f} < max {np.max(K[i,:]):.6f}" + + +def test_build_covariance_nugget(): + """Nugget increases diagonal by exactly the nugget value.""" + N = 5 + rng = np.random.default_rng(99) + points = rng.uniform(-5, 5, (N, 3)).astype(np.float64) + A = _make_identity_anisotropy(N) + + K_no_nugget = build_micro_covariance(points, A, kernel_range=2.0, nugget=0.0) + K_with_nugget = build_micro_covariance(points, A, kernel_range=2.0, nugget=0.1) + + np.testing.assert_allclose( + np.diag(K_with_nugget) - np.diag(K_no_nugget), + 0.1, + atol=1e-14 + ) + + +# ---------------------------------------------------------------- +# solve_micro_weights +# ---------------------------------------------------------------- + +def test_solve_and_evaluate_roundtrip_identity(): + """Solve K@w = residuals, then evaluate back at micro points -> should match residuals.""" + N = 4 + rng = np.random.default_rng(42) + points = rng.uniform(-3, 3, (N, 3)).astype(np.float64) + A = _make_identity_anisotropy(N) + residuals = np.array([0.5, -0.3, 1.2, -0.8], dtype=np.float64) + + weights = solve_micro_weights(points, residuals, A, kernel_range=2.0, nugget=1e-6) + correction = evaluate_micro_correction(points, points, weights, A, kernel_range=2.0) + + np.testing.assert_allclose(correction, residuals, rtol=1e-5) + + +def test_solve_micro_weights_produces_finite_weights(): + """Solving with anisotropic matrices should produce finite, non-NaN weights.""" + N = 4 + rng = np.random.default_rng(99) + points = rng.uniform(-3, 3, (N, 3)).astype(np.float64) + residuals = np.array([0.5, -0.3, 1.2, -0.8], dtype=np.float64) + + grad = rng.normal(0, 1, (N, 3)).astype(np.float64) + grad = grad / np.linalg.norm(grad, axis=1, keepdims=True) + A = compute_anisotropy_matrices_from_gradients(points, grad, r_vertical=0.5, r_lateral=5.0) + + weights = solve_micro_weights(points, residuals, A, kernel_range=2.0, nugget=1e-4) + assert np.all(np.isfinite(weights)), "Weights should be finite" + assert np.all(weights != 0), "Weights should be non-zero" + + correction = evaluate_micro_correction(points, points, weights, A, kernel_range=2.0) + assert np.all(np.isfinite(correction)), "Correction evaluation should be finite" + + +def test_solve_micro_weights_far_apart_points(): + """When points are far apart, weights should approximate residuals (diagonal-dominant K).""" + N = 3 + points = np.array([ + [0.0, 0.0, 0.0], + [100.0, 0.0, 0.0], + [0.0, 100.0, 0.0], + ], dtype=np.float64) + A = _make_identity_anisotropy(N) + residuals = np.array([1.0, 2.0, 3.0], dtype=np.float64) + + weights = solve_micro_weights(points, residuals, A, kernel_range=1.0) + np.testing.assert_allclose(weights, residuals, rtol=1e-3) + + +# ---------------------------------------------------------------- +# compute_anisotropy_matrices_from_gradients +# ---------------------------------------------------------------- + +def test_anisotropy_matrices_shape(): + N = 3 + points = np.random.default_rng(1).uniform(0, 1, (N, 3)).astype(np.float64) + grad = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], dtype=np.float64) + + A = compute_anisotropy_matrices_from_gradients(points, grad) + + assert A.shape == (N, 3, 3) + + +def test_anisotropy_vertical_gradient_produces_expected_scaling(): + """With vertical gradient, X and Y axes get lateral scale, Z gets vertical scale.""" + r_v, r_l = 0.5, 5.0 + points = np.array([[0.0, 0.0, 0.0]], dtype=np.float64) + grad = np.array([[0.0, 0.0, 1.0]], dtype=np.float64) + + A = compute_anisotropy_matrices_from_gradients(points, grad, r_v, r_l) + + x_transformed = A[0] @ np.array([1.0, 0.0, 0.0]) + z_transformed = A[0] @ np.array([0.0, 0.0, 1.0]) + + np.testing.assert_allclose(np.linalg.norm(x_transformed), 1.0 / r_l, rtol=1e-6) + np.testing.assert_allclose(np.linalg.norm(z_transformed), 1.0 / r_v, rtol=1e-6) + + +def test_anisotropy_matrix_is_minimum_stretch(): + """Anisotropy matrices produce ovals not flat lines (determinant > 0).""" + N = 5 + rng = np.random.default_rng(42) + points = rng.uniform(0, 1, (N, 3)).astype(np.float64) + grad = rng.normal(0, 1, (N, 3)).astype(np.float64) + + A = compute_anisotropy_matrices_from_gradients(points, grad, r_vertical=0.3, r_lateral=5.0) + + for i in range(N): + det = np.linalg.det(A[i]) + assert det > 0, f"Matrix {i} has non-positive determinant: {det}" diff --git a/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_macro_integration.py b/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_macro_integration.py new file mode 100644 index 00000000..ceef8c70 --- /dev/null +++ b/tests/test_common/test_modules/test_evaluator/test_micro_anisotropic_macro_integration.py @@ -0,0 +1,730 @@ +import os + +import numpy as np + +from gempy_engine.core.data import Orientations +from gempy_engine.core.data.internal_structs import SolverInput +from gempy_engine.API.interp_single._interp_scalar_field import ( + _solve_interpolation, + _evaluate_sys_eq, +) +from gempy_engine.modules.data_preprocess._input_preparation import ( + surface_points_preprocess, + orientations_preprocess, +) +from gempy_engine.modules.evaluator.micro_anisotropic_evaluator import ( + compute_anisotropy_matrices_from_gradients, + solve_micro_weights, +) + +PLOT = os.getenv("GEMPY_PLOT_MICRO", "1") == "1" + +_MICRO_SURFACE_COLORS = {0: "#00bfff", 1: "#ff6b35"} +_MACRO_SURFACE_COLORS = {0: "#0099cc", 1: "#cc5500"} + + +def _build_grid_2d(x_range, y_range, nx, ny): + x = np.linspace(*x_range, nx) + y = np.linspace(*y_range, ny) + xv, yv = np.meshgrid(x, y) + return np.column_stack([xv.ravel(), yv.ravel()]) + + +def _eval_at_points(sp_internal, ori_internal, options, weights, xyz): + eval_in = SolverInput(sp_internal, ori_internal, xyz_to_interpolate=xyz, fault_internal=None) + options.evaluation_options.compute_scalar_gradient = True + return _evaluate_sys_eq(eval_in, weights, options) + + +def test_micro_correction_moves_contacts_closer_to_target(simple_model_2): + sp, orientations, options, data_descriptor = simple_model_2 + orientations.dip_positions = np.array([[ 0., 4.], [ 4., 1.]]) + orientations.dip_gradients = np.array([[ -.2, .8], [ 0, 1.]]) + options.kernel_options.range = 20 + + options.evaluation_options.compute_scalar_gradient = True + + sp_internal = surface_points_preprocess(sp, data_descriptor.tensors_structure) + ori_internal = orientations_preprocess(orientations) + + n_per_surface = data_descriptor.tensors_structure.number_of_points_per_surface + macro_sp_coords = sp.sp_coords + macro_ori_positions = orientations.dip_positions + macro_sp_surface_ids = np.concatenate([ + np.full(n_per_surface[0], 0, dtype=int), + np.full(n_per_surface[1], 1, dtype=int), + ]) + + solver_input = SolverInput(sp_internal, ori_internal, xyz_to_interpolate=None, fault_internal=None) + macro_weights = _solve_interpolation(solver_input, options.kernel_options) + + # --- target scalars: median macro scalar at original surface points --- + exported_macro_sp = _eval_at_points(sp_internal, ori_internal, options, macro_weights, macro_sp_coords) + macro_at_sp = exported_macro_sp.scalar_field + macro_sp_gx = exported_macro_sp.gx_field + macro_sp_gy = exported_macro_sp.gy_field + macro_sp_gradients = np.column_stack([macro_sp_gx, macro_sp_gy]) + target_per_surface = [ + float(np.median(macro_at_sp[:n_per_surface[0]])), + float(np.median(macro_at_sp[n_per_surface[0]:])), + ] + print(f"target S0 = {target_per_surface[0]:.3f} | target S1 = {target_per_surface[1]:.3f}") + + # --- micro contacts --- + contacts = np.array([ + [1.0, 2.3], [2.0, 2.5], [3.0, 1.5], + [0.5, 1.8], [1.4, 0.2], [2.5, 1.0], + ], dtype=np.float64) + contact_surface_ids = np.array([1, 1, 0, 1, 0, 0], dtype=int) + + exported_contacts = _eval_at_points(sp_internal, ori_internal, options, macro_weights, contacts) + macro_values_at_contacts = exported_contacts.scalar_field + contact_gx = exported_contacts.gx_field + contact_gy = exported_contacts.gy_field + contact_gradients = np.column_stack([contact_gx, contact_gy]) + + target_values_at_contacts = np.array([target_per_surface[sid] for sid in contact_surface_ids]) + contact_residuals = target_values_at_contacts - macro_values_at_contacts + + for i in range(len(contacts)): + print(f" contact {i} (S{contact_surface_ids[i]}): target={target_values_at_contacts[i]:.3f} " + f"macro={macro_values_at_contacts[i]:.3f} residual={contact_residuals[i]:.3f}") + + # --- build augmented micro system (Option 3): contacts + macro points as zero constraints --- + constraint_points = np.vstack([contacts, macro_sp_coords]) + constraint_gradients = np.vstack([contact_gradients, macro_sp_gradients]) + constraint_residuals = np.concatenate([ + contact_residuals, + np.zeros(len(macro_sp_coords)), + ]) + n_contacts = len(contacts) + n_macro = len(macro_sp_coords) + + micro_kernel_range = 0.5 + A = compute_anisotropy_matrices_from_gradients( + constraint_points, constraint_gradients, r_vertical=.5, r_lateral=5.0, + ) + all_weights = solve_micro_weights(constraint_points, constraint_residuals, A, + kernel_range=micro_kernel_range, nugget=1e-6, + kernel_type="matern_5_2") + + print(f" micro weights: contacts {np.array2string(all_weights[:n_contacts], precision=3)}, " + f"macro {np.array2string(all_weights[n_contacts:], precision=3)}") + + # --- grid evaluation --- + grid_xy = _build_grid_2d((-1, 5), (-1, 5), 40, 40) + options.evaluation_options.compute_scalar_gradient = False + macro_fields = _eval_at_points(sp_internal, ori_internal, options, macro_weights, grid_xy) + + micro = options.evaluation_options.micro_anisotropic + micro.enabled = True + micro.points = constraint_points + micro.weights = all_weights + micro.anisotropy_matrices = A + micro.kernel_range = micro_kernel_range + + micro_fields = _eval_at_points(sp_internal, ori_internal, options, macro_weights, grid_xy) + + macro_field_2d = macro_fields.scalar_field.reshape(40, 40) + micro_field_2d = micro_fields.scalar_field.reshape(40, 40) + diff_field = micro_field_2d - macro_field_2d + + assert np.all(np.isfinite(micro_field_2d)), "Micro field is not finite" + assert np.all(np.isfinite(diff_field)), "Diff field is not finite" + max_abs_diff = np.max(np.abs(diff_field)) + assert max_abs_diff > 1e-6, f"Micro correction should produce nonzero change, got max abs diff = {max_abs_diff}" + + # --- contact compliance --- + micro_exported = _eval_at_points(sp_internal, ori_internal, options, macro_weights, contacts) + options.evaluation_options.micro_anisotropic.enabled = False + corrected_contacts = micro_exported.scalar_field + rms_before = np.sqrt(np.mean(contact_residuals ** 2)) + rms_after = np.sqrt(np.mean((target_values_at_contacts - corrected_contacts) ** 2)) + assert rms_after < rms_before, ( + f"Micro correction should reduce contact RMS error. " + f"Before: {rms_before:.6f}, After: {rms_after:.6f}" + ) + + # --- macro point preservation --- + options.evaluation_options.compute_scalar_gradient = False + micro.enabled = True # re-enable for this eval + macro_after_exported = _eval_at_points(sp_internal, ori_internal, options, macro_weights, macro_sp_coords) + micro.enabled = False + macro_after_sp = macro_after_exported.scalar_field + macro_drift = np.abs(macro_after_sp - macro_at_sp) + max_macro_drift = np.max(macro_drift) + mean_macro_drift = np.mean(macro_drift) + assert max_macro_drift < 1.0, ( + f"Macro points shifted too much by micro correction. " + f"Max drift: {max_macro_drift:.4f}, Mean: {mean_macro_drift:.4f}" + ) + + print(f"RMS before: {rms_before:.6f}, RMS after: {rms_after:.6f}") + print(f"Macro point drift — max: {max_macro_drift:.4f}, mean: {mean_macro_drift:.4f}") + + if PLOT: + _plot_results( + grid_xy, macro_field_2d, micro_field_2d, diff_field, + contacts, contact_surface_ids, + macro_values_at_contacts, corrected_contacts, target_values_at_contacts, + macro_sp_coords, macro_sp_surface_ids, n_per_surface, macro_ori_positions, + micro_kernel_range, A, target_per_surface, + n_contacts, macro_before=macro_at_sp, macro_after=macro_after_sp, + ) + + +# ---------------------------------------------------------------- +# plotting +# ---------------------------------------------------------------- +def _plot_results(grid_xy, macro_field, micro_field, diff, + contacts, contact_surface_ids, + macro_vals, corrected_vals, target_vals, + macro_sp_coords, macro_sp_surface_ids, n_per_surface, + macro_ori_positions, micro_kernel_range, A_matrices, + target_per_surface, n_contacts, + macro_before=None, macro_after=None): + import matplotlib.pyplot as plt + + x = grid_xy[:, 0].reshape(macro_field.shape) + y = grid_xy[:, 1].reshape(macro_field.shape) + xlim = (-1, 6) + ylim = (-1, 5) + + vmin = min(macro_field.min(), micro_field.min()) + vmax = max(macro_field.max(), micro_field.max()) + levels = np.linspace(vmin, vmax, 25) + + fig, axes = plt.subplots(2, 2, figsize=(14, 11)) + + # --- Top-left: Macro field --- + ax = axes[0, 0] + ax.set_title("Macro scalar field") + ax.contourf(x, y, macro_field, levels=levels, cmap="viridis", extend="both") + _draw_macro_input(ax, macro_sp_coords, n_per_surface, macro_ori_positions) + + for sv, label in zip(target_per_surface, ["S0 target", "S1 target"]): + ax.contour(x, y, macro_field, levels=[sv], colors="white", linewidths=1.5, linestyles="-") + _draw_contacts(ax, contacts, contact_surface_ids, contact_vals=macro_vals) + + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.legend(loc="upper right", fontsize=6) + ax.set_aspect("equal") + + # --- Top-right: Micro field --- + ax = axes[0, 1] + ax.set_title(f"Micro-adjusted scalar field (range={micro_kernel_range})") + ax.contourf(x, y, micro_field, levels=levels, cmap="viridis", extend="both") + for sv in target_per_surface: + ax.contour(x, y, micro_field, levels=[sv], colors="white", linewidths=1.5, linestyles="-") + + _draw_anisotropy_ellipses(ax, contacts, A_matrices[:n_contacts], micro_kernel_range, + color="yellow", alpha=0.25) + _draw_macro_input(ax, macro_sp_coords, n_per_surface, macro_ori_positions) + _draw_contacts(ax, contacts, contact_surface_ids, contact_vals=corrected_vals) + + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.legend(loc="upper right", fontsize=6) + ax.set_aspect("equal") + + # --- Bottom-left: Difference field --- + ax = axes[1, 0] + c = ax.contourf(x, y, diff, cmap="RdBu_r", levels=20, extend="both") + plt.colorbar(c, ax=ax, shrink=0.9) + ax.set_title("Micro - Macro difference") + _draw_macro_input(ax, macro_sp_coords, n_per_surface, macro_ori_positions) + _draw_contacts(ax, contacts, contact_surface_ids) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.set_aspect("equal") + + # --- Bottom-right: Residual bar chart + macro drift --- + ax = axes[1, 1] + before_abs = np.abs(macro_vals - target_vals) + after_abs = np.abs(corrected_vals - target_vals) + n = len(contacts) + x_idx = np.arange(n) + width = 0.35 + colors_before = [_MICRO_SURFACE_COLORS[sid] for sid in contact_surface_ids] + colors_after = [_MACRO_SURFACE_COLORS[sid] for sid in contact_surface_ids] + ax.bar(x_idx - width/2, before_abs, width, color=colors_before, label="|macro - target|") + ax.bar(x_idx + width/2, after_abs, width, color=colors_after, label="|corrected - target|") + ax.set_xticks(x_idx) + ax.set_xticklabels([f"c{i}\n(S{contact_surface_ids[i]})" for i in range(n)]) + ax.set_title("Contact residual error (abs)") + ax.legend(loc="upper left", fontsize=7) + + if macro_before is not None and macro_after is not None: + drift_text = (f"macro pt drift:\n" + f" max: {np.max(np.abs(macro_after - macro_before)):.4f}\n" + f" mean: {np.mean(np.abs(macro_after - macro_before)):.4f}") + ax.text(0.95, 0.90, drift_text, transform=ax.transAxes, + fontsize=7, verticalalignment="top", horizontalalignment="right", + bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.7)) + + plt.tight_layout() + plt.show() + + +def _draw_contacts(ax, contacts, contact_surface_ids, contact_vals=None): + for sid in [0, 1]: + mask = contact_surface_ids == sid + if mask.any(): + ax.plot(contacts[mask, 0], contacts[mask, 1], "o", + color=_MICRO_SURFACE_COLORS[sid], markersize=8, + markeredgecolor="black", label=f"micro contacts S{sid}") + if contact_vals is not None: + for i in range(len(contacts)): + ax.annotate(f"{contact_vals[i]:.2f}", (contacts[i, 0], contacts[i, 1]), + textcoords="offset points", xytext=(5, 5), fontsize=7, + color=_MICRO_SURFACE_COLORS[contact_surface_ids[i]]) + + +def _draw_anisotropy_ellipses(ax, points, A_matrices, kernel_range, color="yellow", alpha=0.3): + from matplotlib.patches import Ellipse + + for i, p in enumerate(points): + A = A_matrices[i] + ATA = A.T @ A + eigvals, eigvecs = np.linalg.eigh(ATA) + eigvals = np.maximum(eigvals, 1e-12) + semi_axes = kernel_range / np.sqrt(eigvals) + angle = np.degrees(np.arctan2(eigvecs[1, 0], eigvecs[0, 0])) + ell = Ellipse( + xy=(p[0], p[1]), + width=2 * semi_axes[0], + height=2 * semi_axes[1], + angle=angle, + facecolor=color, + edgecolor="black", + alpha=alpha, + linewidth=0.5, + ) + ax.add_patch(ell) + + +def _draw_macro_input(ax, macro_sp_coords, n_per_surface, macro_ori_positions): + s0 = slice(0, n_per_surface[0]) + s1 = slice(n_per_surface[0], n_per_surface[0] + n_per_surface[1]) + ax.plot(macro_sp_coords[s0, 0], macro_sp_coords[s0, 1], "s", + color=_MACRO_SURFACE_COLORS[0], markersize=5, markeredgecolor="black", + label=f"macro S0 pts") + ax.plot(macro_sp_coords[s1, 0], macro_sp_coords[s1, 1], "s", + color=_MACRO_SURFACE_COLORS[1], markersize=5, markeredgecolor="black", + label=f"macro S1 pts") + ax.plot(macro_ori_positions[:, 0], macro_ori_positions[:, 1], "^", + color="magenta", markersize=6, markeredgecolor="black", label="macro orientations") + + +# ---------------------------------------------------------------- +# 3D integration test +# ---------------------------------------------------------------- +def _build_grid_3d(x_range, y_range, z_range, nx, ny, nz): + x = np.linspace(*x_range, nx) + y = np.linspace(*y_range, ny) + z = np.linspace(*z_range, nz) + xv, yv, zv = np.meshgrid(x, y, z, indexing="ij") + return np.column_stack([xv.ravel(), yv.ravel(), zv.ravel()]) + + +def test_micro_correction_moves_3d_contacts_closer_to_target(simple_model): + """3D analog of the 2D integration test. + + Uses the existing simple_model fixture (7 surface points, 2 orientations, + 3D cubic kernel) and exercises the same pipeline: macro solve -> macro scalar/gradient + at contacts and surface points -> per-surface median target -> augmented micro solve + with macro zero constraints -> grid evaluation -> assertions. + + The dense NumPy micro solve here is a reference implementation. Production + path for 3D is intended to be PyKeOps matvec + CG, not dense K assembly. + """ + sp, orientations, options, data_descriptor = simple_model + + options.evaluation_options.compute_scalar_gradient = True + + from gempy_engine.core.data.internal_structs import SolverInput + from gempy_engine.API.interp_single._interp_scalar_field import ( + _solve_interpolation, + _evaluate_sys_eq, + ) + from gempy_engine.modules.data_preprocess._input_preparation import ( + surface_points_preprocess, + orientations_preprocess, + ) + + sp_internal = surface_points_preprocess(sp, data_descriptor.tensors_structure) + ori_internal = orientations_preprocess(orientations) + + n_per_surface = data_descriptor.tensors_structure.number_of_points_per_surface + macro_sp_coords = sp.sp_coords + + # --- macro solve --- + solver_input = SolverInput(sp_internal, ori_internal, xyz_to_interpolate=None, fault_internal=None) + macro_weights = _solve_interpolation(solver_input, options.kernel_options) + + # --- target scalars: median macro scalar at original surface points --- + def _eval_3d(xyz): + proxy = SolverInput(sp_internal, ori_internal, xyz_to_interpolate=xyz, fault_internal=None) + options.evaluation_options.compute_scalar_gradient = True + return _evaluate_sys_eq(proxy, macro_weights, options) + + exported_macro_sp = _eval_3d(macro_sp_coords) + macro_at_sp = exported_macro_sp.scalar_field + macro_sp_gx = exported_macro_sp.gx_field + macro_sp_gy = exported_macro_sp.gy_field + macro_sp_gz = exported_macro_sp.gz_field + macro_sp_gradients = np.column_stack([macro_sp_gx, macro_sp_gy, macro_sp_gz]) + + target_per_surface = [float(np.median(macro_at_sp[:n_per_surface[0]]))] + print(f"target S0 = {target_per_surface[0]:.3f}") + + # --- micro contacts (3D, placed close to the macro interface) --- + # Contacts are near existing macro surface points but shifted in z + # (stratigraphic direction) to produce small, realistic residuals. + contacts = np.array([ + [0.48, 0.55, 0.38], # near SP1 [0.50, 0.50, 0.375], above + [0.68, 0.52, 0.48], # near SP2 [0.667, 0.50, 0.417], above + [0.60, 0.48, 0.37], # near SP5 [0.583, 0.50, 0.392], below + [0.72, 0.46, 0.55], # near SP6 [0.733, 0.50, 0.500], above + ], dtype=np.float64) + contact_surface_ids = np.zeros(len(contacts), dtype=int) # all target the single surface + + exported_contacts = _eval_3d(contacts) + macro_values_at_contacts = exported_contacts.scalar_field + contact_gx = exported_contacts.gx_field + contact_gy = exported_contacts.gy_field + contact_gz = exported_contacts.gz_field + contact_gradients = np.column_stack([contact_gx, contact_gy, contact_gz]) + + target_values_at_contacts = np.array([target_per_surface[sid] for sid in contact_surface_ids]) + contact_residuals = target_values_at_contacts - macro_values_at_contacts + + for i in range(len(contacts)): + print(f" contact {i}: target={target_values_at_contacts[i]:.3f} " + f"macro={macro_values_at_contacts[i]:.3f} residual={contact_residuals[i]:.3f}") + + # --- build micro constraint system --- + # When preserve_macro_points=True: contacts + macro SP as zero-residual constraints. + # When False: contacts only. + micro = options.evaluation_options.micro_anisotropic + micro.preserve_macro_points = False # this need to be false no question + preserve = micro.preserve_macro_points + + if preserve: + constraint_points = np.vstack([contacts, macro_sp_coords]) + constraint_gradients = np.vstack([contact_gradients, macro_sp_gradients]) + constraint_residuals = np.concatenate([ + contact_residuals, + np.zeros(len(macro_sp_coords)), + ]) + else: + constraint_points = contacts + constraint_gradients = contact_gradients + constraint_residuals = contact_residuals + + n_contacts = len(contacts) + n_macro = len(macro_sp_coords) if preserve else 0 + n_constraints = len(constraint_points) + + micro_kernel_range = .1 # larger than 2D case because macro domain is ~0.5 + A = compute_anisotropy_matrices_from_gradients( + constraint_points, constraint_gradients, r_vertical=.4, r_lateral=.7, + ) + all_weights = solve_micro_weights(constraint_points, constraint_residuals, A, + kernel_range=micro_kernel_range, nugget=1e-6, + kernel_type="exponential") + + print(f" micro weights: contacts {np.array2string(all_weights[:n_contacts], precision=3)}" + f"{', macro ' + np.array2string(all_weights[n_contacts:], precision=3) if preserve else ''}") + + # --- grid evaluation (small 3D dense grid) --- + grid_xyz = _build_grid_3d( + (0.25, 0.75), (0.45, 0.55), (0.3, 0.6), 16, 8, 16, + ) + + options.evaluation_options.compute_scalar_gradient = False + macro_fields = _eval_3d(grid_xyz) + + micro = options.evaluation_options.micro_anisotropic + micro.enabled = True + micro.points = constraint_points + micro.weights = all_weights + micro.anisotropy_matrices = A + micro.kernel_range = micro_kernel_range + + micro_fields = _eval_3d(grid_xyz) + + macro_field = macro_fields.scalar_field + micro_field = micro_fields.scalar_field + diff_field = micro_field - macro_field + + assert np.all(np.isfinite(micro_field)), "3D micro field is not finite" + assert np.all(np.isfinite(diff_field)), "3D diff field is not finite" + max_abs_diff = np.max(np.abs(diff_field)) + assert max_abs_diff > 1e-8, f"3D micro correction should produce nonzero change, got max abs diff = {max_abs_diff}" + + # --- contact compliance --- + micro.enabled = True + micro_exported = _eval_3d(contacts) + micro.enabled = False + corrected_contacts = micro_exported.scalar_field + rms_before = np.sqrt(np.mean(contact_residuals ** 2)) + rms_after = np.sqrt(np.mean((target_values_at_contacts - corrected_contacts) ** 2)) + assert rms_after < rms_before, ( + f"3D micro correction should reduce contact RMS error. " + f"Before: {rms_before:.6f}, After: {rms_after:.6f}" + ) + + # --- macro point preservation (only when macro SP are constraints) --- + if preserve: + options.evaluation_options.compute_scalar_gradient = False + micro.enabled = True + macro_after_exported = _eval_3d(macro_sp_coords) + micro.enabled = False + macro_after_sp = macro_after_exported.scalar_field + macro_drift = np.abs(macro_after_sp - macro_at_sp) + max_macro_drift = np.max(macro_drift) + mean_macro_drift = np.mean(macro_drift) + assert max_macro_drift < 2.0, ( + f"3D macro points shifted too much by micro correction. " + f"Max drift: {max_macro_drift:.4f}, Mean: {mean_macro_drift:.4f}" + ) + else: + max_macro_drift = 0.0 + mean_macro_drift = 0.0 + + print(f"3D RMS before: {rms_before:.6f}, RMS after: {rms_after:.6f}") + print(f"3D macro point drift — max: {max_macro_drift:.4f}, mean: {mean_macro_drift:.4f}") + + if PLOT: + _plot_3d_results( + grid_xyz, macro_field, micro_field, diff_field, + contacts, macro_values_at_contacts, corrected_contacts, target_values_at_contacts, + macro_sp_coords, n_per_surface, A, target_per_surface, n_contacts, + micro_kernel_range=micro_kernel_range, + macro_before=macro_at_sp, + macro_after=macro_after_sp if preserve else None, + ) + _plot_3d_pyvista( + grid_xyz, macro_field, micro_field, + contacts, macro_sp_coords, target_per_surface, + ) + + +# ---------------------------------------------------------------- +# 3D plotting (slice-based) +# ---------------------------------------------------------------- +def _plot_3d_results(grid_xyz, macro_field, micro_field, diff, + contacts, macro_vals_before, macro_vals_after, target_vals, + macro_sp_coords, n_per_surface, A_matrices, target_per_surface, + n_contacts, micro_kernel_range=1.0, macro_before=None, macro_after=None): + import matplotlib.pyplot as plt + from matplotlib.patches import Ellipse + + # extract grid shape + x = grid_xyz[:, 0] + y = grid_xyz[:, 1] + z = grid_xyz[:, 2] + nx = len(np.unique(x)) + ny = len(np.unique(y)) + nz = len(np.unique(z)) + shape = (nx, ny, nz) + + macro_3d = macro_field.reshape(shape) + micro_3d = micro_field.reshape(shape) + diff_3d = diff.reshape(shape) + + # take a mid-y slice + y_idx = ny // 2 + x_grid = np.unique(x) + z_grid = np.unique(z) + + macro_slice = macro_3d[:, y_idx, :].T # shape (nz, nx) + micro_slice = micro_3d[:, y_idx, :].T + diff_slice = diff_3d[:, y_idx, :].T + + vmin = min(macro_slice.min(), micro_slice.min()) + vmax = max(macro_slice.max(), micro_slice.max()) + levels = np.linspace(vmin, vmax, 20) + xlim = (x_grid[0], x_grid[-1]) + zlim = (z_grid[0], z_grid[-1]) + + fig, axes = plt.subplots(2, 2, figsize=(14, 10)) + + # --- macro slice --- + ax = axes[0, 0] + ax.set_title("3D Macro scalar field (mid-y slice)") + ax.contourf(x_grid, z_grid, macro_slice, levels=levels, cmap="viridis", extend="both") + ax.plot(macro_sp_coords[:, 0], macro_sp_coords[:, 2], "ks", + markersize=6, markeredgecolor="white", label="macro SP (projected)") + for tv in target_per_surface: + ax.contour(x_grid, z_grid, macro_slice, levels=[tv], colors="white", + linewidths=1.5, linestyles="-") + ax.plot(contacts[:, 0], contacts[:, 2], "ro", markersize=6, + markeredgecolor="black", label="contacts (projected)") + for i in range(n_contacts): + ax.annotate(f"{macro_vals_before[i]:.2f}", (contacts[i, 0], contacts[i, 2]), + textcoords="offset points", xytext=(4, 4), fontsize=6, color="red") + _draw_3d_anisotropy_disks(ax, contacts, A_matrices[:n_contacts], micro_kernel_range) + ax.set_xlim(*xlim) + ax.set_ylim(*zlim) + ax.set_aspect("equal") + ax.legend(fontsize=6) + + # --- micro slice --- + ax = axes[0, 1] + ax.set_title("3D Micro-adjusted scalar field (mid-y slice)") + ax.contourf(x_grid, z_grid, micro_slice, levels=levels, cmap="viridis", extend="both") + for tv in target_per_surface: + ax.contour(x_grid, z_grid, micro_slice, levels=[tv], colors="white", + linewidths=1.5, linestyles="-") + ax.plot(macro_sp_coords[:, 0], macro_sp_coords[:, 2], "ks", + markersize=6, markeredgecolor="white", label="macro SP") + ax.plot(contacts[:, 0], contacts[:, 2], "ro", markersize=6, + markeredgecolor="black", label="contacts") + for i in range(n_contacts): + ax.annotate(f"{macro_vals_after[i]:.2f}", (contacts[i, 0], contacts[i, 2]), + textcoords="offset points", xytext=(4, 4), fontsize=6, color="red") + _draw_3d_anisotropy_disks(ax, contacts, A_matrices[:n_contacts], micro_kernel_range) + ax.set_xlim(*xlim) + ax.set_ylim(*zlim) + ax.set_aspect("equal") + ax.legend(fontsize=6) + + # --- diff slice --- + ax = axes[1, 0] + c = ax.contourf(x_grid, z_grid, diff_slice, cmap="RdBu_r", levels=16, extend="both") + plt.colorbar(c, ax=ax, shrink=0.9) + ax.set_title("3D Micro - Macro difference (mid-y slice)") + ax.plot(macro_sp_coords[:, 0], macro_sp_coords[:, 2], "ks", + markersize=6, markeredgecolor="white", label="macro SP") + ax.plot(contacts[:, 0], contacts[:, 2], "ro", + markersize=6, markeredgecolor="black", label="contacts") + ax.set_xlim(*xlim) + ax.set_ylim(*zlim) + ax.set_aspect("equal") + ax.legend(fontsize=6) + + # --- residual bar chart --- + ax = axes[1, 1] + before_abs = np.abs(macro_vals_before - target_vals) + after_abs = np.abs(macro_vals_after - target_vals) + n = n_contacts + x_idx = np.arange(n) + width = 0.35 + ax.bar(x_idx - width / 2, before_abs, width, color="#00bfff", label="|macro - target|") + ax.bar(x_idx + width / 2, after_abs, width, color="#0099cc", label="|corrected - target|") + ax.set_xticks(x_idx) + ax.set_xticklabels([f"c{i}" for i in range(n)]) + ax.set_title("3D Contact residual error (abs)") + ax.legend(loc="upper left", fontsize=7) + + if macro_before is not None and macro_after is not None: + drift_text = (f"macro pt drift:\n" + f" max: {np.max(np.abs(macro_after - macro_before)):.4f}\n" + f" mean: {np.mean(np.abs(macro_after - macro_before)):.4f}") + ax.text(0.95, 0.85, drift_text, transform=ax.transAxes, + fontsize=7, verticalalignment="top", horizontalalignment="right", + bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.7)) + + plt.tight_layout() + plt.show() + + +def _draw_3d_anisotropy_disks(ax, points, A_matrices, kernel_range, color="yellow"): + """Project 3D anisotropy ellipsoids onto the xz plane as ellipses. + + Draws three contours per contact: + - filled core (0.15 * range) — where correction is strong + - solid e^-1 (1.00 * range) — half-decay contour + - dashed (3.00 * range) — where correction is ~5 % + + A_i is (3,3). The projected 2x2 metric on the xz plane is: + M = (A_i[:, [0,2]])^T @ A_i[:, [0,2]] + The ellipse is {v : v^T M v = radius^2}. + """ + from matplotlib.patches import Ellipse + + for i, p in enumerate(points): + A = A_matrices[i] + B = A[:, [0, 2]] + M = B.T @ B + eigvals, eigvecs = np.linalg.eigh(M) + eigvals = np.maximum(eigvals, 1e-12) + angle = np.degrees(np.arctan2(eigvecs[1, 0], eigvecs[0, 0])) + + def _ellipse(radius, **kwargs): + semi = radius / np.sqrt(eigvals) + return Ellipse(xy=(p[0], p[2]), width=2 * semi[0], height=2 * semi[1], + angle=angle, **kwargs) + + ax.add_patch(_ellipse(kernel_range * 0.15, + facecolor=color, edgecolor="none", alpha=0.35, zorder=2)) + ax.add_patch(_ellipse(kernel_range, + fill=False, edgecolor=color, linewidth=1.2, + alpha=0.7, zorder=2, label="e⁻¹")) + ax.add_patch(_ellipse(kernel_range * 3.0, + fill=False, edgecolor=color, linewidth=0.6, linestyle="--", + alpha=0.3, zorder=2, label="e⁻³")) + + # deduplicate legend entries + handles, labels = ax.get_legend_handles_labels() + seen = set() + unique = [(h, l) for h, l in zip(handles, labels) if l not in seen and not seen.add(l)] + ax.legend(handles=[h for h, _ in unique], labels=[l for _, l in unique], + fontsize=6, loc="upper right") + + +def _plot_3d_pyvista(grid_xyz, macro_field, micro_field, + contacts, macro_sp_coords, target_per_surface): + """3D PyVista visualization of the adjusted scalar field. + + Shows: semi-transparent scalar volume, target isosurfaces, + macro surface points (black), and micro contacts (red). + """ + try: + import pyvista as pv + except ImportError: + print("pyvista not installed, skipping 3D plot") + return + + x = np.unique(grid_xyz[:, 0]) + y = np.unique(grid_xyz[:, 1]) + z = np.unique(grid_xyz[:, 2]) + nx, ny, nz = len(x), len(y), len(z) + + xyz_reshaped = grid_xyz.reshape(nx, ny, nz, 3) + xv = xyz_reshaped[..., 0] + yv = xyz_reshaped[..., 1] + zv = xyz_reshaped[..., 2] + + grid = pv.StructuredGrid(xv, yv, zv) + + macro_3d = macro_field.reshape(nx, ny, nz) + micro_3d = micro_field.reshape(nx, ny, nz) + grid["macro_scalar"] = macro_3d.ravel(order="F") + grid["micro_scalar"] = micro_3d.ravel(order="F") + + p = pv.Plotter() + p.add_text("Micro-adjusted scalar field — volume + isosurfaces", font_size=10) + + p.add_mesh(grid, scalars="micro_scalar", opacity=0.3, cmap="viridis", + show_edges=False, show_scalar_bar=False) + + for target in target_per_surface: + contour = grid.contour(isosurfaces=[target], scalars="micro_scalar") + p.add_mesh(contour, color="white", opacity=0.85, show_edges=False, + label=f"target isosurface ({target:.3f})") + + p.add_mesh(pv.PolyData(macro_sp_coords), color="black", + point_size=12, render_points_as_spheres=True, + label="macro SP") + + p.add_mesh(pv.PolyData(contacts), color="red", + point_size=16, render_points_as_spheres=True, + label="micro contacts") + + p.add_axes() + p.add_legend() + p.show()