diff --git a/docs/octree_refinement.md b/docs/octree_refinement.md new file mode 100644 index 00000000..4c965159 --- /dev/null +++ b/docs/octree_refinement.md @@ -0,0 +1,79 @@ +# Octree refinement support modes + +```python +from gempy_engine.core.data import OctreeRefinementMode + +options.evaluation_options.number_octree_levels_surface = 4 +options.evaluation_options.octree_refinement_mode = OctreeRefinementMode.PRECISE +``` + +The mode also accepts the strings `"fast"`, `"balanced"`, and `"precise"`. + +| Mode | Support stencil | Isolated interior parent count | +| --- | --- | --- | +| `fast` (default) | No additional support; existing selector | 1 | +| `balanced` | Six face neighbors | 7 | +| `precise` | 26 face, edge, and corner neighbors | 27 | + +This setting does not change octree depth, minimum-level refinement, curvature +thresholds, or categorical surface selection. One dilation of their combined +selection is applied per transition, only while generating surface extraction +levels and only when mesh extraction is enabled. Support cells never seed another +dilation in the same transition. Each selected parent generates eight children. + +`balanced` is a cost/coverage compromise, **not 2:1 octree balancing**. It does not +provide diagonal support around primal edges. `precise` supplies the full touching +neighborhood among existing sparse cells, clipped to the physical domain. Neither +mode guarantees discovery of unsampled features or watertight surfaces. A surface +newly detected at the outer edge of an earlier support band can still request an +absent branch; this produces a warning rather than silently claiming closure. +Physical extent capping and extraction-mask stitching are separate concerns. + +## Diagnostics + +With `options.debug = True` (or evaluation `verbose`), generated octree grids expose +`refinement_debug`: primary surface, additional, and support-only parent masks; +their counts; final count; closure multiplier; generated children; and unique +missing in-domain support requests. These masks index the **previous** generation, +like `active_cells`. Counts for surface and additional selection may overlap. + +For non-fast modes, or when `options.debug` is true, each extracted mesh exposes +`support_report`. This CPU diagnostic enumerates unique sampled sign-changing +primal edges before extraction masking. Missing incident cells are classified as +physical, mask-removed, or never generated. Records include edge direction and +integer coordinate, missing coordinates, and the last existing ancestor level. +An edge can have more than one boundary classification. Internal refinement +failures emit a warning; tests can assert the corresponding count is zero. + +The report checks true sign changes, not the extractor's slightly extrapolated +edge-intersection tolerance. It does not audit later fault/overlap triangle +removal or all final mesh edge incidences. Debug reports and coordinate sets add +CPU memory/time overhead, especially for many surfaces. + +Integer extraction coordinates now use signed int64 and theoretical domain +bounds rather than byte packing and maximum active coordinates. Normal fast-mode +selection is unchanged; meshes affected by the old coordinate/bounds defects can +change even in fast mode. + +## Performance + +For an interior planar one-cell selection, the full halo approaches 3x the parent +count; isolated selections can reach 27x (7x in balanced mode). Each selected +parent contributes eight centers and 64 stored corner rows at the next level. +No corner deduplication or interpolation caching is introduced here. + +A NumPy float64 lookup-only smoke benchmark on a dense `32 x 32 x 32` lattice +gave the following counts (not an end-to-end interpolation benchmark): + +| Selection | Primary | Balanced final | Precise final | +| --- | ---: | ---: | ---: | +| Plane `ix == 16` | 1,024 | 3,072 (3x) | 3,072 (3x) | +| Sphere band `abs(norm(coord - 15.5) - 9) < 0.5` | 1,032 | 2,696 (2.61x) | 4,024 (3.90x) | +| Isolated seeds `all(coord % 4 == 2)` | 512 | 3,584 (7x) | 13,824 (27x) | + +Single-call lookup times on the development machine were 0.36–0.63 ms for +balanced and 1.05–1.68 ms for precise; these are indicative, not performance +thresholds or measurements of interpolation/reporting overhead. + +The default remains fast. Representative curved, multi-stack, faulted, and GPU +time/memory benchmarks are still needed before recommending a different default. diff --git a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py index 35955806..654a0cf5 100644 --- a/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py +++ b/gempy_engine/API/dual_contouring/multi_scalar_dual_contouring.py @@ -1,4 +1,5 @@ import copy +import warnings from typing import List, Any import numpy as np @@ -22,6 +23,8 @@ from ...modules.dual_contouring.dual_contouring_interface import (find_intersection_on_edge, get_triangulation_codes, get_masked_codes, mask_generation) from ...modules.dual_contouring.overlapping import average_overlapping_vertices, remove_fault_overlap_triangles +from ...modules.dual_contouring._support_report import mesh_support_report +from ...core.data.options.evaluation_options import OctreeRefinementMode @gempy_profiler_decorator @@ -106,6 +109,7 @@ def dual_contouring_multi_scalar( surface_to_stack = [] # track which stack each surface belongs to # Generate meshes for each scalar field dc_data_per_surface_all = [] + support_reports = [] stack_relations = data_descriptor.stack_structure.masking_descriptor for n_scalar_field in range(data_descriptor.stack_structure.n_stacks): if stack_relations[n_scalar_field] is StackRelationType.NULL_SPACE: @@ -114,6 +118,24 @@ def dual_contouring_multi_scalar( mask = all_mask_arrays[n_scalar_field] n_surfaces_to_export = output.scalar_field_at_sp.shape[0] for surface_i in range(n_surfaces_to_export): + report = None + if options.debug or options.evaluation_options.octree_refinement_mode != OctreeRefinementMode.FAST: + report = mesh_support_report( + left_right_codes, + output.exported_fields.scalar_field[output.grid.corners_grid_slice], + output.scalar_field_at_sp[surface_i], base_number, mask, + surface_index=surface_i, + ancestor_coordinates=[level.grid.octree_grid.integer_coordinates for level in octree_list[:-1]] + ) + report['stack_index'] = n_scalar_field + if report['internal_refinement_boundary_edge_count']: + warnings.warn( + f"Stack {n_scalar_field}, surface {surface_i}: " + f"{report['internal_refinement_boundary_edge_count']} crossing edges lack " + "in-domain refinement support; see mesh.support_report.", + RuntimeWarning, stacklevel=2 + ) + support_reports.append(report) valid_edges = all_valid_edges[n_scalar_field] valid_edges_per_surface = valid_edges.reshape((n_surfaces_to_export, -1, 12)) slice_object = _surface_slicer(surface_i, valid_edges_per_surface) @@ -153,6 +175,8 @@ def dual_contouring_multi_scalar( dc_data_list=dc_data_per_surface_all, max_workers=None ) + for mesh, report in zip(all_meshes, support_reports): + mesh.support_report = report # endregion # Save differentiable vertices before in-place overlap modifications, # then replace mesh.vertices with a detached copy so averaging doesn't diff --git a/gempy_engine/API/interp_single/interp_features.py b/gempy_engine/API/interp_single/interp_features.py index 7fb08aa9..48d83899 100644 --- a/gempy_engine/API/interp_single/interp_features.py +++ b/gempy_engine/API/interp_single/interp_features.py @@ -28,7 +28,8 @@ def interpolate_n_octree_levels(interpolation_input: InterpolationInput, options grid_1_centers: EngineGrid = get_next_octree_grid( prev_octree=next_octree, evaluation_options=options.evaluation_options, - current_octree_level=i + current_octree_level=i, + debug=options.debug ) interpolation_input.set_temp_grid(grid_1_centers) octree_list.append(next_octree) diff --git a/gempy_engine/core/data/__init__.py b/gempy_engine/core/data/__init__.py index 3bfa6adc..2715672e 100644 --- a/gempy_engine/core/data/__init__.py +++ b/gempy_engine/core/data/__init__.py @@ -3,4 +3,5 @@ from .kernel_classes.orientations import Orientations, OrientationsInternals from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals from .options.interpolation_options import InterpolationOptions +from .options.evaluation_options import OctreeRefinementMode from .solutions import Solutions diff --git a/gempy_engine/core/data/dual_contouring_mesh.py b/gempy_engine/core/data/dual_contouring_mesh.py index e2d26eb3..692ad997 100644 --- a/gempy_engine/core/data/dual_contouring_mesh.py +++ b/gempy_engine/core/data/dual_contouring_mesh.py @@ -10,6 +10,7 @@ class DualContouringMesh: vertices: np.ndarray edges: np.ndarray dc_data: Optional[DualContouringData] = None # * In principle we need this just for testing + support_report: Optional[dict] = None def __repr__(self): return f"DualContouringMesh({self.vertices.shape[0]} vertices, {self.edges.shape[0]} edges)" diff --git a/gempy_engine/core/data/options/evaluation_options.py b/gempy_engine/core/data/options/evaluation_options.py index 67e3de64..b302004c 100644 --- a/gempy_engine/core/data/options/evaluation_options.py +++ b/gempy_engine/core/data/options/evaluation_options.py @@ -8,6 +8,13 @@ from .micro_anisotropic_options import MicroAnisotropicOptions +class OctreeRefinementMode(str, enum.Enum): + """Support stencil, not 2:1 balancing or a guarantee of watertightness.""" + FAST = "fast" + BALANCED = "balanced" + PRECISE = "precise" + + class MeshExtractionMaskingOptions(enum.Enum): NOTHING = enum.auto() # * This is only for testing DISJOINT = enum.auto() @@ -22,6 +29,7 @@ class EvaluationOptions: octree_curvature_threshold: float = -1. #: Threshold to do octree refinement due to curvature to deal with angular geometries. This curvature assumes that 1 is the maximum curvature of any voxel octree_error_threshold: float = 1. #: Number of standard deviations to consider a voxel as candidate to refine octree_min_level: int = 2 + octree_refinement_mode: OctreeRefinementMode = OctreeRefinementMode.FAST mesh_extraction: bool = True mesh_extraction_masking_options: MeshExtractionMaskingOptions = MeshExtractionMaskingOptions.INTERSECT diff --git a/gempy_engine/core/data/regular_grid.py b/gempy_engine/core/data/regular_grid.py index 2f72de7c..80cabd71 100644 --- a/gempy_engine/core/data/regular_grid.py +++ b/gempy_engine/core/data/regular_grid.py @@ -16,6 +16,8 @@ class RegularGrid: _active_cells: np.ndarray = field(default=None, repr=False, init=False) left_right: np.ndarray = field(default=None, repr=False, init=False) + _integer_coordinates: np.ndarray = field(default=None, repr=False, init=False) + refinement_debug: dict = field(default=None, repr=False, init=False) values: np.ndarray = field(default=None, repr=False, init=False) original_values: np.ndarray = field(default=None, repr=False, init=False) #: When the regular grid is representing a octree level, only active cells are stored in values. This is the original values of the regular grid. @@ -68,6 +70,10 @@ def from_octree_level(cls, xyz_coords_octree: np.ndarray, previous_regular_grid: regular_grid_for_octree_level.values = xyz_coords_octree # ! Overwrite the common values regular_grid_for_octree_level._active_cells = active_cells regular_grid_for_octree_level.left_right = left_right + regular_grid_for_octree_level._integer_coordinates = ( + 2 * BackendTensor.tfnp.repeat(previous_regular_grid.integer_coordinates[active_cells], 8, axis=0) + + BackendTensor.t.array(left_right, dtype='int64') + ) return regular_grid_for_octree_level @@ -86,6 +92,16 @@ def from_schema(cls, schema: GridSchema): left_right=None ) + @property + def integer_coordinates(self): + """Signed lattice coordinates in sparse row order (Z fastest at root).""" + if self._integer_coordinates is None: + axes = [BackendTensor.arange(int(n), dtype='int64') for n in self.regular_grid_shape] + self._integer_coordinates = BackendTensor.t.stack( + BackendTensor.t.meshgrid(*axes, indexing='ij'), axis=-1 + ).reshape(-1, 3) + return self._integer_coordinates + @property def active_cells(self) -> np.ndarray: if self._active_cells is not None: diff --git a/gempy_engine/modules/dual_contouring/_support_report.py b/gempy_engine/modules/dual_contouring/_support_report.py new file mode 100644 index 00000000..08a32bac --- /dev/null +++ b/gempy_engine/modules/dual_contouring/_support_report.py @@ -0,0 +1,69 @@ +"""CPU-only topology diagnostics, independent of triangle-candidate filtering.""" +from itertools import product + +import numpy as np + +from ...core.backend_tensor import BackendTensor + + +def mesh_support_report(coordinates, scalar_corners, isovalue, domain_shape, + mask=None, surface_index=0, ancestor_coordinates=()): + """Classify missing incident cells for unique sampled sign-changing edges. + + This diagnoses sampled crossings, not unsampled components or later triangle + removal by overlap/fault processing. Counts of missing cells are incidences. + """ + to_numpy = BackendTensor.t.to_numpy + coords = np.asarray(to_numpy(coordinates), dtype=np.int64) + scalar = np.asarray(to_numpy(scalar_corners)).reshape(-1, 8) + iso = float(isovalue) + retained = np.ones(len(coords), dtype=bool) if mask is None else np.asarray(to_numpy(mask), dtype=bool) + generated = set(map(tuple, coords)) + kept = set(map(tuple, coords[retained])) + ancestors = [set(map(tuple, to_numpy(c))) for c in ancestor_coordinates] + bounds = tuple(int(n) for n in domain_shape) + corners = np.array(list(product((0, 1), repeat=3)), dtype=np.int64) + edges = set() + for direction, pairs in enumerate(( + ((0, 4), (1, 5), (2, 6), (3, 7)), + ((0, 2), (1, 3), (4, 6), (5, 7)), + ((0, 1), (2, 3), (4, 5), (6, 7)), + )): + for a, b in pairs: + crossing = (scalar[:, a] >= iso) != (scalar[:, b] >= iso) + edges.update((direction, *p) for p in coords[crossing] + corners[a]) + report = dict(surface_index=surface_index, crossing_edge_count=len(edges), + missing_incident_cell_count=0, physical_boundary_edge_count=0, + mask_boundary_edge_count=0, internal_refinement_boundary_edge_count=0, + violations=[]) + for edge in sorted(edges): + direction, *origin = edge + transverse = [i for i in range(3) if i != direction] + missing = [] + kinds = set() + for offsets in product((-1, 0), repeat=2): + cell = list(origin) + for axis, offset in zip(transverse, offsets): + cell[axis] += offset + cell = tuple(cell) + if cell in kept: + continue + outside = any(c < 0 or c >= n for c, n in zip(cell, bounds)) + kind = 'physical' if outside else 'mask' if cell in generated else 'refinement' + kinds.add(kind) + stopped = None + if kind == 'refinement': + for level, existing in enumerate(ancestors): + shift = len(ancestors) - level + if tuple(c >> shift for c in cell) in existing: + stopped = level + missing.append(dict(coordinate=cell, kind=kind, outside_extent=outside, + ancestor_stop_level=stopped)) + if missing: + report['missing_incident_cell_count'] += len(missing) + for kind, key in (('physical', 'physical_boundary_edge_count'), + ('mask', 'mask_boundary_edge_count'), + ('refinement', 'internal_refinement_boundary_edge_count')): + report[key] += kind in kinds + report['violations'].append(dict(direction=direction, coordinate=tuple(origin), missing=missing)) + return report diff --git a/gempy_engine/modules/dual_contouring/fancy_triangulation.py b/gempy_engine/modules/dual_contouring/fancy_triangulation.py index ab6ed9b3..638af625 100644 --- a/gempy_engine/modules/dual_contouring/fancy_triangulation.py +++ b/gempy_engine/modules/dual_contouring/fancy_triangulation.py @@ -4,100 +4,15 @@ def get_left_right_array(octree_list: list[OctreeLevel]): - dtype = bool - match BackendTensor.engine_backend: - case BackendTensor.engine_backend.PYTORCH: - dtype = BackendTensor.tfnp.bool - case BackendTensor.engine_backend.numpy: - dtype = bool - case _: - raise ValueError("Unsupported backend") - - # === Local function === - def _compute_voxel_binary_code(root_bits_list, dir_idx: int, left_right_all, voxel_select_all): - - # Calculate the voxels from root - processed_root_bits = [] - for bit_array in root_bits_list: - idx_curr = bit_array - for active_voxels_per_lvl in voxel_select_all: - idx_curr = BackendTensor.tfnp.repeat(idx_curr[active_voxels_per_lvl], 8, axis=0) - processed_root_bits.append(idx_curr) - - left_right_list = [] - voxel_select_op = list(voxel_select_all[1:]) - voxel_select_op.append(BackendTensor.tfnp.ones( - left_right_all[-1].shape[0], - dtype=dtype - ) - ) - left_right_all = left_right_all[::-1] - voxel_select_op = voxel_select_op[::-1] - - for e, left_right_per_lvl in enumerate(left_right_all): # size is equal to the depth of the tree (except root) - left_right_per_lvl_dir = left_right_per_lvl[:, dir_idx] - for n_rep in range(e): - inner = left_right_per_lvl_dir[voxel_select_op[e - n_rep]] - left_right_per_lvl_dir = BackendTensor.tfnp.repeat(inner, 8, axis=0) - left_right_list.append(left_right_per_lvl_dir) - - # Combine refinement bits (LSB->MSB) with root bits (LSB->MSB) - final_list = left_right_list + processed_root_bits - binary_code = BackendTensor.tfnp.stack(final_list) - return binary_code - - # === Local function === - - if len(octree_list) == 1: - # * Not only that, the current implementation only works with pure octree starting at [2,2,2] - raise ValueError("Octree list must have more than one level") - - voxel_select_all = [octree_iter.grid.octree_grid.active_cells for octree_iter in octree_list[1:]] - left_right_all = [octree_iter.grid.octree_grid.left_right for octree_iter in octree_list[1:]] - - # Dynamic generation of root indices - import numpy as np - root_res = octree_list[0].grid.octree_grid_shape - nx, ny, nz = int(root_res[0]), int(root_res[1]), int(root_res[2]) - - # Generate coordinate grids (Order: Z fast, Y, X slow) - x_indices = np.repeat(np.arange(nx), ny * nz) - y_indices = np.tile(np.repeat(np.arange(ny), nz), nx) - z_indices = np.tile(np.arange(nz), nx * ny) - - def get_root_bits_list(indices): - max_val = max(nx, ny, nz) - # Calculate needed bits (at least 1) - n_bits = int(max_val - 1).bit_length() if max_val > 1 else 1 - - bits_list = [] - for i in range(n_bits): - # Extract bit i (LSB to MSB) - bit_val = (indices >> i) & 1 - bits_list.append(BackendTensor.tfnp.array(bit_val, dtype=dtype)) - return bits_list - - binary_x = _compute_voxel_binary_code(get_root_bits_list(x_indices), 0, left_right_all, voxel_select_all) - binary_y = _compute_voxel_binary_code(get_root_bits_list(y_indices), 1, left_right_all, voxel_select_all) - binary_z = _compute_voxel_binary_code(get_root_bits_list(z_indices), 2, left_right_all, voxel_select_all) - - bool_to_int_x = BackendTensor.tfnp.packbits(binary_x, axis=0, bitorder="little") - bool_to_int_y = BackendTensor.tfnp.packbits(binary_y, axis=0, bitorder="little") - bool_to_int_z = BackendTensor.tfnp.packbits(binary_z, axis=0, bitorder="little") - left_right_array = BackendTensor.tfnp.vstack([bool_to_int_x, bool_to_int_y, bool_to_int_z]).T - - if left_right_array.shape[0] == 0: - base_x = base_y = base_z = 0 - else: - base_x = bool_to_int_x.max() + 1 - base_y = bool_to_int_y.max() + 1 - base_z = bool_to_int_z.max() + 1 - - return left_right_array, (base_x, base_y, base_z) + """Return signed sparse coordinates and theoretical (not active) bounds.""" + grid = octree_list[-1].grid.octree_grid + return grid.integer_coordinates, tuple(int(n) for n in grid.regular_grid_shape) def _get_pack_factors(base_x, base_y, base_z): """Generates [base_y * base_z, base_z, 1] for packing 3D coordinates.""" + if int(base_x) * int(base_y) * int(base_z) > 2 ** 63 - 1: + raise OverflowError("Octree domain exceeds signed int64 coordinate-key capacity") # Ensure we use int64 for packing to avoid overflow bx = BackendTensor.tfnp.array(base_x, dtype='int64') by = BackendTensor.tfnp.array(base_y, dtype='int64') @@ -346,7 +261,10 @@ def _get_indices_via_searchsorted(voxel_code, compressed_0, compressed_1, compre """ Memory-efficient replacement for broadcasting, preserving exact original indices. """ - vc_1d = voxel_code.squeeze() + vc_1d = voxel_code.reshape(-1) + if len(vc_1d) == 0: + empty = BackendTensor.tfnp.zeros(0, dtype='int64') + return empty, empty, empty # 1. Get sorting indices to map back to original positions later sort_indices = BackendTensor.tfnp.argsort(vc_1d) diff --git a/gempy_engine/modules/octrees_topology/_neighbor_closure.py b/gempy_engine/modules/octrees_topology/_neighbor_closure.py new file mode 100644 index 00000000..1b76816f --- /dev/null +++ b/gempy_engine/modules/octrees_topology/_neighbor_closure.py @@ -0,0 +1,39 @@ +"""One-step sparse lattice support; never seeds a second halo from support cells.""" +from itertools import product +import math + +from ...core.backend_tensor import BackendTensor +from ...core.data.options.evaluation_options import OctreeRefinementMode + + +def close_refinement_mask(coordinates, primary_mask, domain_shape, mode): + """Return support-only mask and number of unique absent in-domain cells.""" + t = BackendTensor.tfnp + mode = OctreeRefinementMode(mode) + support = t.zeros(primary_mask.shape, dtype=bool) + if mode == OctreeRefinementMode.FAST or not t.any(primary_mask): + return support, 0 + shape = tuple(int(n) for n in domain_shape) + if math.prod(shape) > 2 ** 63 - 1: + raise OverflowError("Octree domain exceeds signed int64 coordinate-key capacity") + bounds = t.array(shape, dtype='int64') + factors = t.array([shape[1] * shape[2], shape[2], 1], dtype='int64') + keys = (coordinates * factors).sum(axis=1) + order = t.argsort(keys) + sorted_keys = keys[order] + seeds = coordinates[primary_mask] + missing = [] + for offset in product((-1, 0, 1), repeat=3): + distance = sum(abs(v) for v in offset) + if distance == 0 or (mode == OctreeRefinementMode.BALANCED and distance != 1): + continue + neighbors = seeds + t.array(offset, dtype='int64') + inside = ((neighbors >= 0) & (neighbors < bounds)).all(axis=1) + requested = (neighbors[inside] * factors).sum(axis=1) + positions = t.searchsorted(sorted_keys, requested) + safe = t.clip(positions, 0, len(sorted_keys) - 1) + found = (positions < len(sorted_keys)) & (sorted_keys[safe] == requested) + support[order[safe[found]]] = True + missing.append(requested[~found]) + missing_count = int(t.unique(t.concatenate(missing)).shape[0]) + return support & ~primary_mask, missing_count diff --git a/gempy_engine/modules/octrees_topology/_octree_internals.py b/gempy_engine/modules/octrees_topology/_octree_internals.py index a9c85f9e..735f2bab 100644 --- a/gempy_engine/modules/octrees_topology/_octree_internals.py +++ b/gempy_engine/modules/octrees_topology/_octree_internals.py @@ -1,19 +1,22 @@ import numpy as np +import warnings from typing import List from ._curvature_analysis import mark_highest_curvature_voxels from ._octree_common import _generate_next_level_centers +from ._neighbor_closure import close_refinement_mask +from ...config import DEBUG_MODE from ...core.backend_tensor import BackendTensor from ...core.data.engine_grid import EngineGrid from ...core.data.exported_fields import ExportedFields from ...core.data.interp_output import InterpOutput from ...core.data.octree_level import OctreeLevel -from ...core.data.options.evaluation_options import EvaluationOptions +from ...core.data.options.evaluation_options import EvaluationOptions, OctreeRefinementMode from ...core.data.regular_grid import RegularGrid def compute_next_octree_locations(prev_octree: OctreeLevel, evaluation_options: EvaluationOptions, - current_octree_level: int) -> EngineGrid: + current_octree_level: int, debug: bool = False) -> EngineGrid: ids = prev_octree.litho_faults_ids_corners_grid uv_8 = ids.reshape((-1, 8)) @@ -27,7 +30,26 @@ def compute_next_octree_locations(prev_octree: OctreeLevel, evaluation_options: prev_octree=prev_octree ) - voxel_select = voxel_select | additional_voxel_selected_to_refinement + surface_mask = voxel_select + primary_mask = surface_mask | additional_voxel_selected_to_refinement + mode = OctreeRefinementMode(evaluation_options.octree_refinement_mode) + closure_enabled = ( + mode != OctreeRefinementMode.FAST + and evaluation_options.mesh_extraction + and current_octree_level < evaluation_options.number_octree_levels_surface - 1 + ) + support_mask, missing_count = close_refinement_mask( + prev_octree.grid.octree_grid.integer_coordinates if closure_enabled else None, + primary_mask, prev_octree.grid.octree_grid.regular_grid_shape, + mode if closure_enabled else OctreeRefinementMode.FAST + ) + voxel_select = primary_mask | support_mask + if missing_count: + warnings.warn( + f"Octree level {current_octree_level}: {missing_count} in-domain support cells " + "are absent from the sparse generation; neighbor closure cannot restore them.", + RuntimeWarning, stacklevel=2 + ) if compute_topology := False: # TODO: Fix topology function raise NotImplementedError @@ -50,6 +72,33 @@ def compute_next_octree_locations(prev_octree: OctreeLevel, evaluation_options: ), ) + if debug or DEBUG_MODE or evaluation_options.verbose: + primary_count = int(primary_mask.sum()) + minimum_level = current_octree_level < evaluation_options.octree_min_level + child_grid = grid_next_centers.octree_grid + corner_offsets = BackendTensor.t.array( + [[x, y, z] for x in (0, 1) for y in (0, 1) for z in (0, 1)], dtype='int64' + ) + # Diagnostics only: count stored rows representing the same lattice point. + corner_coordinates = (child_grid.integer_coordinates[:, None, :] + corner_offsets).reshape(-1, 3) + unique_corners = np.unique(BackendTensor.t.to_numpy(corner_coordinates), axis=0).shape[0] + grid_next_centers.octree_grid.refinement_debug = dict( + primary_surface_mask=surface_mask, + additional_refinement_mask=additional_voxel_selected_to_refinement, + support_only_mask=support_mask, + current_cell_count=len(primary_mask), + primary_surface_count=int(surface_mask.sum()), + additional_refinement_count=int(additional_voxel_selected_to_refinement.sum()), + minimum_level_count=len(primary_mask) if minimum_level else 0, + curvature_count=0 if minimum_level else int(additional_voxel_selected_to_refinement.sum()), + support_only_count=int(support_mask.sum()), + final_refinement_count=int(voxel_select.sum()), + closure_multiplier=int(voxel_select.sum()) / primary_count if primary_count else 1.0, + missing_current_neighbor_count=missing_count, + generated_child_count=len(xyz_coords), + duplicate_corner_count=len(corner_coordinates) - unique_corners + ) + if True: grid_next_centers.debug_vals = (xyz_coords, xyz_anchor, shift_select_xyz, bool_idx, voxel_select, grid_next_centers) return grid_next_centers # TODO: This is going to break the tests that were using this @@ -126,4 +175,4 @@ def _additional_refinement_tests(voxel_select, current_octree_level, evaluation_ curvature_threshold=evaluation_options.octree_curvature_threshold # * This curvature assumes that 1 is the maximum curvature of any voxel ) - return BackendTensor.t.array(additional_voxel_selected_to_refinement) \ No newline at end of file + return BackendTensor.t.array(additional_voxel_selected_to_refinement) diff --git a/gempy_engine/modules/octrees_topology/octrees_topology_interface.py b/gempy_engine/modules/octrees_topology/octrees_topology_interface.py index 79b90d7c..6f84b865 100644 --- a/gempy_engine/modules/octrees_topology/octrees_topology_interface.py +++ b/gempy_engine/modules/octrees_topology/octrees_topology_interface.py @@ -19,11 +19,12 @@ def get_next_octree_grid(prev_octree: OctreeLevel, evaluation_options: EvaluationOptions, - current_octree_level: int = 9999) -> EngineGrid: + current_octree_level: int = 9999, debug: bool = False) -> EngineGrid: octree_from_output: EngineGrid = compute_next_octree_locations( prev_octree=prev_octree, evaluation_options=evaluation_options, - current_octree_level=current_octree_level + current_octree_level=current_octree_level, + debug=debug ) return octree_from_output diff --git a/tests/test_common/test_modules/test_neighbor_closure.py b/tests/test_common/test_modules/test_neighbor_closure.py new file mode 100644 index 00000000..ce4172f4 --- /dev/null +++ b/tests/test_common/test_modules/test_neighbor_closure.py @@ -0,0 +1,172 @@ +from itertools import product +from types import SimpleNamespace + +import numpy as np +import pytest + +from gempy_engine.config import AvailableBackends +from gempy_engine.core.backend_tensor import BackendTensor +from gempy_engine.core.data.regular_grid import RegularGrid +from gempy_engine.core.data.options.evaluation_options import EvaluationOptions +from gempy_engine.modules.octrees_topology._neighbor_closure import close_refinement_mask +from gempy_engine.modules.octrees_topology._octree_common import _generate_next_level_centers +from gempy_engine.modules.octrees_topology._octree_internals import compute_next_octree_locations +from gempy_engine.modules.dual_contouring.fancy_triangulation import get_left_right_array +from gempy_engine.modules.dual_contouring.fancy_triangulation import triangulate +from gempy_engine.modules.dual_contouring._support_report import mesh_support_report + + +@pytest.fixture(params=['numpy', 'PYTORCH']) +def backend(request): + if request.param == 'PYTORCH': + pytest.importorskip('torch') + old = BackendTensor.engine_backend + old_gpu = BackendTensor.use_gpu + old_dtype = BackendTensor.dtype + old_keops = BackendTensor.use_pykeops + BackendTensor._change_backend(engine_backend=getattr(AvailableBackends, request.param), use_gpu=False, dtype='float64') + yield + BackendTensor._change_backend(engine_backend=old, use_gpu=old_gpu, dtype=old_dtype, use_pykeops=old_keops) + + +@pytest.mark.parametrize('mode,seed,expected', [ + ('fast', (2, 2, 2), 1), ('balanced', (2, 2, 2), 7), + ('precise', (2, 2, 2), 27), ('precise', (0, 2, 2), 18), + ('precise', (0, 0, 2), 12), ('precise', (0, 0, 0), 8), +]) +def test_stencil(backend, mode, seed, expected): + t = BackendTensor.t + coords = t.array(list(product(range(5), range(6), range(7))), dtype='int64') + primary = (coords == t.array(seed)).all(axis=1) + support, missing = close_refinement_mask(coords, primary, (5, 6, 7), mode) + assert int((support | primary).sum()) == expected + assert not bool((support & primary).any()) + assert missing == 0 + + +def test_sheet_overlap_and_empty(backend): + t = BackendTensor.t + coords = t.array(list(product(range(7), repeat=3)), dtype='int64') + primary = coords[:, 0] == 3 + support, missing = close_refinement_mask(coords, primary, (7, 7, 7), 'precise') + np.testing.assert_array_equal(t.to_numpy(support | primary), t.to_numpy(abs(coords[:, 0] - 3) <= 1)) + assert missing == 0 + support, missing = close_refinement_mask(coords[:0], primary[:0], (7, 7, 7), 'precise') + assert len(support) == missing == 0 + + +def test_missing_and_deep_coordinates(backend): + t = BackendTensor.t + coords = t.array([[255, 1, 1], [256, 1, 1], [257, 1, 1]], dtype='int64') + primary = t.array([True, False, False], dtype=bool) + support, missing = close_refinement_mask(coords, primary, (1024, 3, 3), 'balanced') + np.testing.assert_array_equal(t.to_numpy(support), [False, True, False]) + assert missing == 5 + + +def test_lineage(backend): + t = BackendTensor.t + root = RegularGrid([0, 300, 0, 2, 0, 2], [300, 2, 2]) + selected = root.integer_coordinates[:, 0] == 255 + xyz, bits = _generate_next_level_centers(root.values[selected], root.dxdydz) + child = RegularGrid.from_octree_level(xyz, root, selected, bits) + expected = 2 * np.repeat(t.to_numpy(root.integer_coordinates[selected]), 8, axis=0) + t.to_numpy(bits) + np.testing.assert_array_equal(t.to_numpy(child.integer_coordinates), expected) + assert len(np.unique(expected, axis=0)) == len(expected) + codes, bounds = get_left_right_array([SimpleNamespace(grid=SimpleNamespace(octree_grid=child))]) + assert bounds == (600, 4, 4) + np.testing.assert_array_equal(t.to_numpy(codes), expected) + geometric = (t.to_numpy(xyz) - t.to_numpy(root.orthogonal_extent)[::2]) / np.array([float(v) for v in child.dxdydz]) - 0.5 + np.testing.assert_allclose(geometric, expected, atol=1e-6) + + +@pytest.mark.parametrize('mode,count', [('fast', 1), ('balanced', 7), ('precise', 27)]) +def test_refinement_and_surface_depth(backend, mode, count): + t = BackendTensor.t + grid = RegularGrid([0, 5, 0, 5, 0, 5], [5, 5, 5]) + ids = t.zeros((125, 8)) + ids[62, 0] = 1 + prev = SimpleNamespace(litho_faults_ids_corners_grid=ids.reshape(-1), + dxdydz=grid.dxdydz, grid=SimpleNamespace(octree_grid=grid), outputs=[]) + options = EvaluationOptions(_number_octree_levels=4, _number_octree_levels_surface=2, + octree_min_level=0, octree_refinement_mode=mode, verbose=True) + result = compute_next_octree_locations(prev, options, 0).octree_grid + assert int(result.active_cells.sum()) == count + assert len(result.values) == count * 8 + assert result.refinement_debug['primary_surface_count'] == 1 + assert result.refinement_debug['support_only_count'] == count - 1 + result = compute_next_octree_locations(prev, options, 1).octree_grid + assert int(result.active_cells.sum()) == 1 + + +def test_support_classification(backend): + t = BackendTensor.t + coords = np.array(list(product(range(3), repeat=3))) + corners = np.array(list(product((0, 1), repeat=3))) + scalar = (coords[:, None, :] + corners)[..., 0].astype(float) + report = mesh_support_report(t.array(coords), t.array(scalar), 1.5, (3, 3, 3)) + assert report['crossing_edge_count'] == 16 + assert report['physical_boundary_edge_count'] == 12 + assert report['internal_refinement_boundary_edge_count'] == 0 + retained = ~np.all(coords == [1, 1, 1], axis=1) + report = mesh_support_report(t.array(coords), t.array(scalar), 1.5, (3, 3, 3), t.array(retained)) + assert report['mask_boundary_edge_count'] == 4 + assert report['internal_refinement_boundary_edge_count'] == 0 + report = mesh_support_report(t.array(coords[retained]), t.array(scalar[retained]), 1.5, (3, 3, 3)) + assert report['internal_refinement_boundary_edge_count'] == 4 + + +def test_closed_sphere_connectivity(backend): + t = BackendTensor.t + corners = t.array(list(product((0, 1), repeat=3)), dtype='int64') + grid = RegularGrid([0, 8, 0, 8, 0, 8], [8, 8, 8]) + options = EvaluationOptions(_number_octree_levels=3, _number_octree_levels_surface=3, + octree_min_level=0, octree_refinement_mode='precise') + for level in range(3): + points = (grid.integer_coordinates[:, None, :] + corners) / (2 ** level) + scalar = ((points - 4.13) ** 2).sum(axis=2) - 2.37 ** 2 + if level < 2: + prev = SimpleNamespace(litho_faults_ids_corners_grid=t.array(scalar >= 0, dtype='int64').reshape(-1), + dxdydz=grid.dxdydz, grid=SimpleNamespace(octree_grid=grid), outputs=[]) + grid = compute_next_octree_locations(prev, options, level).octree_grid + report = mesh_support_report(grid.integer_coordinates, scalar, 0., grid.regular_grid_shape) + assert report['internal_refinement_boundary_edge_count'] == 0 + assert report['physical_boundary_edge_count'] == 0 + pairs = [(0, 4), (1, 5), (2, 6), (3, 7), (0, 2), (1, 3), (4, 6), (5, 7), + (0, 1), (2, 3), (4, 5), (6, 7)] + valid = t.stack([(scalar[:, a] >= 0) != (scalar[:, b] >= 0) for a, b in pairs], axis=1) + active = valid.any(axis=1) + coords = grid.integer_coordinates[active] + vertices = t.array(coords, dtype='float64') + 0.5 + triangles = triangulate(coords, valid[active], 3, t.ones((len(coords), 12, 3)), + vertices, tuple(int(n) for n in grid.regular_grid_shape)) + triangles = t.to_numpy(triangles) + edges = np.sort(np.concatenate([triangles[:, [0, 1]], triangles[:, [1, 2]], triangles[:, [2, 0]]]), axis=1) + _, counts = np.unique(edges, axis=0, return_counts=True) + assert len(counts) > 0 + assert np.all(counts == 2) + + +@pytest.mark.parametrize('mode', ['balanced', 'precise']) +def test_fault_model_mode_integration(graben_fault_model, mode): + from gempy_engine.API.model.model_api import compute_model + + interpolation_input, structure, options = graben_fault_model + options.evaluation_options.number_octree_levels = 4 + options.evaluation_options.octree_refinement_mode = mode + options.debug = True + solutions = compute_model(interpolation_input, options, structure) + assert solutions.dc_meshes + for mesh in solutions.dc_meshes: + assert mesh.support_report is not None + assert 'stack_index' in mesh.support_report + + +def test_mode_serialization(): + from gempy_engine.core.data import InterpolationOptions, OctreeRefinementMode + + options = InterpolationOptions.from_args(range=1., c_o=1.) + assert options.evaluation_options.octree_refinement_mode == OctreeRefinementMode.FAST + options.evaluation_options.octree_refinement_mode = OctreeRefinementMode.PRECISE + restored = InterpolationOptions.model_validate_json(options.model_dump_json()) + assert restored.evaluation_options.octree_refinement_mode == OctreeRefinementMode.PRECISE