Add micro-point support in interpolation process

MICRO_ANISOTROPIC_FIELD_DEFORMATION.md

@@ -1,6 +1,9 @@
 # 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.
+This document records the historical prototype and design exploration. For the
+current authored-contact API and limitations, see [docs/micro_points.md](docs/micro_points.md).
+The historical covariance solve and gradient-derived metrics below now live in
+`tests/test_common/test_modules/test_evaluator/micro_reference.py`, not the production API.
 
 ## Goal
 
@@ -333,7 +336,9 @@ micro_weights = solve_micro_weights(
 )
 ```
 
-10. Store on `options.evaluation_options.micro_anisotropic` and evaluate final field.
+10. Prototype-only: apply fitted weights explicitly; the current authored API
+    uses `options.micro_options` for settings and a stack-local `MicroCorrection`
+    for fitted state (see [docs/micro_points.md](docs/micro_points.md)).
 
 ### 3D Frame Construction
 

docs/micro_points.md

@@ -0,0 +1,66 @@
+# Authored micro points
+
+Micro contacts add a local anisotropic scalar-field correction after the ordinary
+macro cokriging solve. They do not enter the macro covariance system. The feature
+is opt-in: authored points are ignored with a warning while
+`options.micro_options.enabled` is `False`.
+
+For a GemPy `GeoModel`, assign contacts to the owning structural element. Each
+4-by-4 support transform has its contact position in the last column and its
+three support axes in the upper-left 3-by-3 block:
+
+```python
+import numpy as np
+import gempy as gp
+
+support = np.eye(4)[None].copy()
+support[0, :3, 3] = [0.5, 0.5, 0.52]
+support[0, :3, :3] = np.diag([0.5, 0.5, 0.5])
+element = model.structural_frame.structural_groups[0].elements[0]
+element.micro_points = gp.data.MicroPointsTable.from_transforms(
+    support, names=element.name
+)
+model.interpolation_options.micro_options.enabled = True
+gp.compute_model(model)
+```
+
+Alternatively, engine callers pass `MicroPoints(points, anisotropy_matrices,
+nuggets, surface_indices)` to `InterpolationInput(..., micro_points=...)`.
+`surface_indices` are zero-based **global** surface indices in structural-frame
+order. Inputs have shapes `(N, 3)`, `(N, 3, 3)`, `(N,)`, and `(N,)` respectively.
+Internal stack subsets rebase surface indices to the owning stack. Model preparation
+validates ownership and supported stack types once, then shares contact queries
+across stacks; the authored arrays are not converted or mutated in place.
+The matrices transform coordinate differences into each contact's local metric;
+GemPy builds them from the inverse support axes after coordinate transforms.
+
+Configure `options.micro_options.kernel_type` (`"exponential"`, `"matern_3_2"`,
+or `"matern_5_2"`), `kernel_range`, `nugget`, `strength`, and
+`preserve_macro_points` as needed. With macro preservation enabled (the default),
+surface points become additional zero-correction constraints. Micro weights are
+fitted separately per stack at each interpolation level; neither inputs nor
+options store fitted state. Set `enabled = False` to deactivate without removing
+authored contacts. PyTorch fits preserve gradients through coordinates, matrices,
+and contact nuggets; nearest metric selection for macro constraints is discrete.
+
+The final scalar field is the faulted macro field plus the stack-local micro
+deformation. Faults do not mask, displace, or split the micro kernel. Contacts
+on fault-affected stratigraphic stacks fit residuals against the macro field
+including upstream ordinary or finite faults. When any contact's stack is enabled,
+every stack evaluates the same coordinate layout: grid, all macro surface points,
+then all authored micro contacts (including contacts on disabled stacks). These
+last coordinates are evaluation queries only, not macro interpolation constraints
+or public grid cells. The fit reads macro values directly from that evaluation;
+there is no separate macro kernel call for the contacts. Upstream fault rows use
+the same suffix, including activation and finite-fault projection. Fault-row
+minimum and finite-fault normal selection use only the grid and macro surface
+point prefix, so appended contacts cannot change the reference frame. Enabled
+authored micro points on fault stacks are rejected, whether or not the fault has
+upstream dependencies. Disabled authored points on fault stacks are ignored.
+The micro fit uses a dense NumPy or PyTorch solve (not PyKeOps); macro PyKeOps
+acceleration remains available and micro evaluation is tensor-native (the dense
+fit and correction are not PyKeOps-accelerated). Flat stacks
+retain fused macro scalar evaluation when possible; gradients use the non-fused
+macro evaluator. Nearest-center metric assignment for preserved macro points is
+discrete. Finite-fault projection uses NumPy and is not end-to-end differentiable.
+Enabled authored micro contacts on external-function stacks remain unsupported.

gempy_engine/API/interp_single/_aux_faults_ops.py

@@ -11,12 +11,10 @@
 def _grab_stack_fault_data(_all_stack_values_block, _interpolation_input_i, _stack_structure, grid_size:int) -> FaultsData:
     fault_data = _interpolation_input_i.fault_values or FaultsData()
     fault_data.fault_values_everywhere = _all_stack_values_block[_stack_structure.active_faults_relations]
-    fv_on_all_sp = fault_data.fault_values_everywhere[:, grid_size:]
+    fv_on_all_sp = fault_data.fault_values_everywhere[:, grid_size:_interpolation_input_i.macro_reference_size]
     fault_data.fault_values_on_sp = fv_on_all_sp[:, _interpolation_input_i.slice_feature]
     return fault_data
 
-
-
 def _options_with_finite_fault_gradients(
         options: InterpolationOptions,
         fault_input: FaultsData | None,
@@ -35,7 +33,8 @@ def _modify_faults_values_output(
         xyz_to_interpolate: np.ndarray,
 ) -> np.ndarray:
     values_on_all_xyz = output.values_on_all_xyz
-    val_min = BackendTensor.t.min(values_on_all_xyz, axis=1).reshape(-1, 1)  # ? Is this as good as it gets?
+    reference_size = output.exported_fields._macro_reference_size or len(xyz_to_interpolate)
+    val_min = BackendTensor.t.min(values_on_all_xyz[:, :reference_size], axis=1).reshape(-1, 1)
     shifted_vals = (values_on_all_xyz - val_min)  # * Shift values between 0 and 1... hopefully
     if not fault_input.finite_fault_defined:
         return shifted_vals
@@ -65,12 +64,12 @@ def _modify_faults_values_output(
     )
 
     gradient_matrix = np.stack(gradients_np, axis=-1)
-    valid_gradient = np.linalg.norm(gradient_matrix, axis=1) > 1e-12
+    valid_gradient = np.linalg.norm(gradient_matrix[:reference_size], axis=1) > 1e-12
     if not np.any(valid_gradient):
         raise ValueError("Cannot determine finite-fault frame from near-zero gradients")
 
     center = np.asarray(fault_input.finite_fault.center)
-    distances_to_center = np.linalg.norm(projected_points - center, axis=1)
+    distances_to_center = np.linalg.norm(projected_points[:reference_size] - center, axis=1)
     center_index = np.argmin(np.where(valid_gradient, distances_to_center, np.inf))
     finite_fault_scalar_np = fault_input.finite_fault.calculate_slip(
         points=projected_points,

gempy_engine/API/interp_single/_interp_scalar_field.py

@@ -102,7 +102,7 @@ def _solve_interpolation_result(
 
 
 def _evaluate_sys_eq(eval_input: Union[SolverInput, EvaluatorInput], weights: np.ndarray, options: InterpolationOptions,
-                     grid: EngineGrid | None = None) -> ExportedFields:
+                      grid: EngineGrid | None = None) -> ExportedFields:
     inverse = None
     if options.evaluation_options.deduplicate_octree_corners:
         eval_input, inverse = _deduplicate_corners(eval_input, grid)

gempy_engine/API/interp_single/_interp_single_feature.py

@@ -16,8 +16,10 @@
 from ...core.data.kernel_classes.faults import FaultsData
 from ...core.data.options import InterpolationOptions
 from ...core.data.scalar_field_output import ScalarFieldOutput
+from ...core.data.stack_relation_type import StackRelationType
 from ...modules.activator import activator_interface
 from ...modules.data_preprocess import data_preprocess_interface
+from ...modules.evaluator.micro_correction import apply_micro_correction, fit_micro_correction
 
 
 def interpolate_feature_with_cokrig(interpolation_input: InterpolationInput,
@@ -36,11 +38,15 @@ def interpolate_feature_with_cokrig(interpolation_input: InterpolationInput,
 
     weights = compute_weights(solver_input, stack_number, options)
     exported_fields: ExportedFields = _evaluate_sys_eq(solver_input, weights, options, grid=grid)
+    correction = fit_micro_correction(interpolation_input, exported_fields.scalar_field_everywhere, options,
+                                      data_shape.number_of_points_per_surface)
+    apply_micro_correction(exported_fields, xyz, correction)
 
     exported_fields.set_structure_values(
         reference_sp_position=data_shape.reference_sp_position,
         slice_feature=interpolation_input.slice_feature,
-        grid_size=interpolation_input.grid.len_all_grids
+        grid_size=interpolation_input.grid.len_all_grids,
+        macro_reference_size=interpolation_input.macro_reference_size
     )
 
     exported_fields.debug = solver_input.debug
@@ -64,7 +70,11 @@ def interpolate_feature_with_external_function(interpolation_input: Interpolatio
         grid = interpolation_input.grid
 
     # region Interpolate scalar field
-    xyz = grid.values
+    shared_queries = (interpolation_input.evaluation_micro_points is not None or
+                      interpolation_input.stack_relation is StackRelationType.FAULT)
+    xyz = (data_preprocess_interface.prepare_grid(grid.values, interpolation_input.all_surface_points,
+                                                   interpolation_input.evaluation_micro_points)
+           if shared_queries else grid.values)
 
     exported_fields: ExportedFields = _interpolate_external_function(
         interp_funct=external_interp_funct,
@@ -74,7 +84,8 @@ def interpolate_feature_with_external_function(interpolation_input: Interpolatio
     exported_fields.set_structure_values(
         reference_sp_position=None,
         slice_feature=None,
-        grid_size=xyz.shape[0]
+        grid_size=grid.len_all_grids,
+        macro_reference_size=interpolation_input.macro_reference_size if shared_queries else None
     )
     output = _segment(exported_fields, external_segment_funct, grid, interpolation_input, options, xyz)
 
@@ -92,7 +103,8 @@ def input_preprocess(data_shape: TensorsStructure, interpolation_input: Interpol
     # * We need to interpolate in ALL the surface points not only the surface points of the stack
     grid_internal: np.ndarray = data_preprocess_interface.prepare_grid(
         grid=grid.values,
-        surface_points=interpolation_input.all_surface_points
+        surface_points=interpolation_input.all_surface_points,
+        micro_points=interpolation_input.evaluation_micro_points
     )
 
     fault_values: FaultsData = interpolation_input.fault_values
@@ -149,11 +161,14 @@ def _scalar_field_segmentation(exported_fields: ExportedFields, external_segment
 
 
 def _interpolate_external_function(interp_funct, xyz):
+    def field(func):
+        return BackendTensor.t.array(func(xyz), dtype=BackendTensor.dtype) if func is not None else None
+
     exported_fields = ExportedFields(
-        _scalar_field=interp_funct.implicit_function(xyz),
-        _gx_field=interp_funct.gx_function(xyz) if interp_funct.gx_function is not None else None,
-        _gy_field=interp_funct.gy_function(xyz) if interp_funct.gy_function is not None else None,
-        _gz_field=interp_funct.gz_function(xyz) if interp_funct.gz_function is not None else None,
+        _scalar_field=field(interp_funct.implicit_function),
+        _gx_field=field(interp_funct.gx_function),
+        _gy_field=field(interp_funct.gy_function),
+        _gz_field=field(interp_funct.gz_function),
         _scalar_field_at_surface_points=interp_funct.scalar_field_at_surface_points
     )
     return exported_fields

gempy_engine/API/interp_single/_multi_scalar_field_manager.py

@@ -19,6 +19,7 @@
 from ...core.data.scalar_field_output import ScalarFieldOutput
 from ...core.data.stack_relation_type import StackRelationType
 from ...core.data.stacks_structure import StacksStructure
+from ...modules.data_preprocess.data_preprocess_interface import prepare_grid
 
 
 # @off
@@ -33,7 +34,7 @@ def interpolate_all_fields(interpolation_input: InterpolationInput, options: Int
     )
 
     if (os.getenv("GEMPY_FLAT_STACKS", "False").lower() in ("true", "1", "t", "y", "yes") and
-            BackendTensor.use_pykeops 
+            BackendTensor.use_pykeops
             # and not has_external_functions
     ):
         all_scalar_fields_outputs: List[ScalarFieldOutput] = _interpolate_stack_flat(data_descriptor, interpolation_input, options)
@@ -57,7 +58,7 @@ def _interpolate_stack_flat(root_data_descriptor: InputDataDescriptor, root_inte
                             options: InterpolationOptions) -> ScalarFieldOutput | List[ScalarFieldOutput]:
     stack_structure = root_data_descriptor.stack_structure
 
-    xyz_to_interpolate_size: int = root_interpolation_input.grid.len_all_grids + root_interpolation_input.surface_points.n_points
+    xyz_to_interpolate_size: int = root_interpolation_input.evaluation_size
     all_stack_values_block: np.ndarray = BackendTensor.t.zeros(
         (stack_structure.n_stacks, xyz_to_interpolate_size),
         dtype=BackendTensor.dtype_obj
@@ -95,7 +96,7 @@ def _interpolate_stack(root_data_descriptor: InputDataDescriptor, root_interpola
 
     all_scalar_fields_outputs: List[ScalarFieldOutput | None] = [None] * stack_structure.n_stacks
 
-    xyz_to_interpolate_size: int = root_interpolation_input.grid.len_all_grids + root_interpolation_input.surface_points.n_points
+    xyz_to_interpolate_size: int = root_interpolation_input.evaluation_size
     all_stack_values_block: np.ndarray = BackendTensor.t.zeros(
         (stack_structure.n_stacks, xyz_to_interpolate_size),
         dtype=BackendTensor.dtype_obj)  # * Used for faults
@@ -113,6 +114,7 @@ def _interpolate_stack(root_data_descriptor: InputDataDescriptor, root_interpola
                 all_interpolation_input=root_interpolation_input,
                 stack_structure=stack_structure
             )
+            fault_input = interpolation_input_i.fault_values
 
             output: ScalarFieldOutput = interpolate_feature_with_external_function(
                 interpolation_input=interpolation_input_i,
@@ -154,7 +156,10 @@ def _interpolate_stack(root_data_descriptor: InputDataDescriptor, root_interpola
             values_output = _modify_faults_values_output(  # ! This is all_STACK_values_block (not all_scalar_fields_outputs)
                 fault_input=fault_input,
                 output=output,
-                xyz_to_interpolate=solver_input.xyz_to_interpolate
+                xyz_to_interpolate=(solver_input.xyz_to_interpolate if stack_structure.interp_function is None else
+                                    prepare_grid(interpolation_input_i.grid.values,
+                                                 interpolation_input_i.all_surface_points,
+                                                 interpolation_input_i.evaluation_micro_points))
             )
             all_stack_values_block[i, :] = values_output
 

gempy_engine/API/interp_single/_stack_ops.py

@@ -22,6 +22,7 @@
 
 from ...modules.activator import activator_interface
 from ...modules.data_preprocess import data_preprocess_interface
+from ...modules.evaluator.micro_correction import apply_micro_correction, fit_micro_correction
 
 
 @dataclass
@@ -91,7 +92,9 @@ def process_chunk(state: InterpolationState, chunk: list[int]):
         state.solver_inputs[i] = chunk_solver_inputs[idx]
 
     # Evaluate this chunk
-    if any(fault_input.finite_fault_defined for fault_input in (item.fault_values for item in chunk_interpolation_inputs)):
+    # The fused PyKeOps evaluator cannot concatenate scalar and gradient lazy kernels.
+    if (any(item.fault_values.finite_fault_defined for item in chunk_interpolation_inputs)
+            or any(opt.compute_scalar_gradient for opt in chunk_options)):
         chunk_eval_inputs, chunk_exported_fields = _evaluate(
             interpolation_inputs=chunk_interpolation_inputs,
             options=state.options,
@@ -114,6 +117,10 @@ def process_chunk(state: InterpolationState, chunk: list[int]):
 
     for idx, i in enumerate(chunk):
         state.eval_inputs[i] = chunk_eval_inputs[idx]
+        correction = fit_micro_correction(chunk_interpolation_inputs[idx],
+                                          chunk_exported_fields[idx].scalar_field_everywhere,
+                                          chunk_options[idx], chunk_tensor_structs[idx].number_of_points_per_surface)
+        apply_micro_correction(chunk_exported_fields[idx], chunk_eval_inputs[idx].xyz_to_interpolate, correction)
 
     # Segment this chunk
     chunk_outputs = _segment(
@@ -238,8 +245,8 @@ def _evaluate(interpolation_inputs: list[InterpolationInput], options: Interpola
 
 
 def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options: InterpolationOptions, solver_inputs, stack_structure: StacksStructure,
-                         tensor_structs: list[TensorsStructure], stack_indices: list[int] | None = None,
-                         options_per_stack: list[InterpolationOptions] | None = None) -> tuple[list[EvaluatorInput], list[ExportedFields]]:
+                          tensor_structs: list[TensorsStructure], stack_indices: list[int] | None = None,
+                          options_per_stack: list[InterpolationOptions] | None = None) -> tuple[list[EvaluatorInput], list[ExportedFields]]:
     from gempy_engine.modules.evaluator.symbolic_evaluator import symbolic_evaluator_optimized_stacked
 
     eval_inputs: list[EvaluatorInput] = []
@@ -270,7 +277,7 @@ def _evaluate_optimized(interpolation_inputs: list[InterpolationInput], options:
     exported_fields_list: list[ExportedFields] = symbolic_evaluator_optimized_stacked(
         eval_inputs=reduced_inputs,
         weights_list=weights_list,
-        options_list=options_list
+        options_list=options_list,
     )
 
     for idx, exported_fields in enumerate(exported_fields_list):
@@ -379,7 +386,6 @@ def _process_external_chunk(state: InterpolationState, chunk: list[int]):
         all_interpolation_input=state.root_interpolation_input,
         stack_structure=state.stack_structure
     )
-
     output: ScalarFieldOutput = interpolate_feature_with_external_function(
         interpolation_input=interpolation_input_i,
         options=options_i,
@@ -388,3 +394,9 @@ def _process_external_chunk(state: InterpolationState, chunk: list[int]):
     )
 
     state.all_scalar_fields_outputs[i] = output
+    if interpolation_input_i.stack_relation is StackRelationType.FAULT:
+        xyz = data_preprocess_interface.prepare_grid(
+            interpolation_input_i.grid.values, interpolation_input_i.all_surface_points,
+            interpolation_input_i.evaluation_micro_points)
+        state.all_stack_values_block[i, :] = _modify_faults_values_output(
+            interpolation_input_i.fault_values, output, xyz)

gempy_engine/API/model/model_api.py

@@ -20,6 +20,7 @@
 from ...core.utils import gempy_profiler_decorator
 from ...core.exceptions import GemPyEngineInputError
 from ...core.data.options.temp_interpolation_values import TempInterpolationValues
+from ...modules.data_preprocess.micro_points import prepare_micro_points
 from ...modules.geophysics.fw_gravity import compute_gravity
 from ...modules.geophysics.fw_magnetic import compute_tmi
 from ...modules.weights_cache.weights_cache_interface import WeightCache
@@ -43,6 +44,8 @@ def compute_model(interpolation_input: InterpolationInput, options: Interpolatio
         # Check input is valid
         _check_input_validity(interpolation_input, options, data_descriptor)  # TODO
 
+        prepare_micro_points(interpolation_input, options, data_descriptor.stack_structure)
+
         output: list[OctreeLevel] = interpolate_n_octree_levels(
             interpolation_input=interpolation_input,
             options=options,