Merge remote-tracking branch 'origin/main'

Author: javoha

docs/developers_notes/dev_log/log_2024-06.md

@@ -1,7 +1,53 @@
 # Refactoring grid to add transform
 
 ## TODO:
-- [ ] Rotation of pyvista grid
-- [ ] Rotation of pyvista meshes
+- [x] Rotation of pyvista grid
+- [x] Rotation of pyvista meshes
+- [ ] Run something with octrees again
+- [ ] **We do need settings presets**
+
 
+## Notes:
+### Optimizing interpolation June 2024
+
+- Create_grad_kernel is call 27 times
+
+### Wholes issue
+- Spremberg with 1.5 std of extra refinement:
+      Number of voxels marked by stats: 424 of torch.Size([512]).
+      Number of voxels marked total by corners: 281
+      Dense Grid would be 512 voxels
+      Number of voxels marked by stats: 2535 of torch.Size([3568]).
+      Number of voxels marked total by corners: 1492
+      Dense Grid would be 4096 voxels
+      Number of voxels marked by stats: 12148 of torch.Size([21152]).
+      Number of voxels marked total by corners: 7614
+      Dense Grid would be 32768 voxels
+      Number of voxels marked by stats: 69619 of torch.Size([101432]).
+      Number of voxels marked total by corners: 37423
+      Dense Grid would be 262144 voxels
+      Chunking done: 12 chunks
+- Spremberg with 1
+        Number of voxels marked by stats: 408 of torch.Size([512]).
+        Number of voxels marked total by corners: 281
+        Dense Grid would be 512 voxels
+        Number of voxels marked by stats: 1953 of torch.Size([3504]).
+        Number of voxels marked total by corners: 1492
+        Dense Grid would be 4096 voxels
+        Number of voxels marked by stats: 8632 of torch.Size([16856]).
+        Number of voxels marked total by corners: 7612
+        Dense Grid would be 32768 voxels
+        Chunking done: 13 chunks
+        Number of voxels marked by stats: 47653 of torch.Size([79080]).
+        Number of voxels marked total by corners: 37408
+        Dense Grid would be 262144 voxels
+        Chunking done: 10 chunks
+
+### Other:
+{'_cached_pivot': [200, 200, 0],
+'_is_default_transform': False,
+'position': array([0., 0., 0.]),
+'rotation': array([  0.,  -0., -45.]),
+'scale': array([1., 1., 1.])}
+
 ## Questions:

examples/tutorials/z_other_tutorials/a1_fold.py

@@ -31,7 +31,7 @@
 geo_model = gp.create_geomodel(
     project_name='tutorial_model',
     extent=[0, 2500, 0, 1000, 0, 1110],
-    refinement=6,
+    refinement=4,
     importer_helper=gp.data.ImporterHelper(
         path_to_orientations=path_to_data + "tutorial_model_orientations.csv",
         path_to_surface_points=path_to_data + "tutorial_model_surface_points.csv"
@@ -62,17 +62,23 @@
 
 interpolation_options: gp.data.InterpolationOptions = geo_model.interpolation_options
 
-
 interpolation_options.kernel_options.range = 1.7
 interpolation_options.evaluation_options.number_octree_levels_surface = 4
-interpolation_options.evaluation_options.curvature_threshold = 0.1
+interpolation_options.evaluation_options.compute_scalar_gradient = False
+interpolation_options.evaluation_options.curvature_threshold = 1
+interpolation_options.evaluation_options.min_octree_level = 1
+
+interpolation_options.evaluation_options.verbose = True
 
 # %%
 # Compute the geological model
 # We use the specified backend (in this case, PyTorch) to compute the model.
 gp.compute_model(
     gempy_model=geo_model,
-    engine_config=gp.data.GemPyEngineConfig(backend=gp.data.AvailableBackends.PYTORCH)
+    engine_config=gp.data.GemPyEngineConfig(
+        backend=gp.data.AvailableBackends.numpy,
+        dtype="float64"
+    )
 )
 
 # %%

gempy/API/compute_API.py

@@ -10,6 +10,7 @@
 from .grid_API import set_custom_grid
 from ..core.data.gempy_engine_config import GemPyEngineConfig
 from ..core.data.geo_model import GeoModel
+from ..modules.data_manipulation.engine_factory import interpolation_input_from_structural_frame
 from ..optional_dependencies import require_gempy_legacy
 
 
@@ -39,7 +40,7 @@ def compute_model(gempy_model: GeoModel, engine_config: Optional[GemPyEngineConf
             )
 
             # TODO: To decide what to do with this.
-            interpolation_input = gempy_model.interpolation_input_copy
+            interpolation_input = interpolation_input_from_structural_frame(gempy_model)
             gempy_model.taped_interpolation_input = interpolation_input  # * This is used for gradient tape
 
             gempy_model.solutions = gempy_engine.compute_model(

gempy/API/grid_API.py

@@ -12,12 +12,11 @@
 def set_section_grid(grid: Grid, section_dict: dict):
     if grid.sections is None:
         grid.sections = Sections(regular_grid=grid.regular_grid, section_dict=section_dict)
-        grid.sections
     else:
         grid.sections.set_sections(section_dict,
                                    regular_grid=grid.regular_grid)
 
-    set_active_grid(grid, [GridTypes.SECTIONS])
+    set_active_grid(grid, [Grid.GridTypes.SECTIONS])
     return grid.sections
 
 
@@ -38,19 +37,19 @@ def set_topography_from_random(grid: Grid, fractal_dimension: float = 2.0, d_z:
         values_2d=random_topography
     )
 
-    set_active_grid(grid, [GridTypes.TOPOGRAPHY])
+    set_active_grid(grid, [Grid.GridTypes.TOPOGRAPHY])
     return grid.topography
 
 
 def set_topography_from_subsurface_structured_grid(grid: Grid, struct: "subsurface.StructuredData"):
     grid.topography = Topography.from_subsurface_structured_data(struct, grid.regular_grid)
-    set_active_grid(grid, [GridTypes.TOPOGRAPHY])
+    set_active_grid(grid, [Grid.GridTypes.TOPOGRAPHY])
     return grid.topography
 
 
 def set_topography_from_arrays(grid: Grid,  xyz_vertices: np.ndarray):
     grid.topography = Topography.from_unstructured_mesh(grid.regular_grid, xyz_vertices)
-    set_active_grid(grid, [GridTypes.TOPOGRAPHY])
+    set_active_grid(grid, [Grid.GridTypes.TOPOGRAPHY])
     return grid.topography
 
 
@@ -67,7 +66,7 @@ def set_custom_grid(grid: Grid, xyz_coord: np.ndarray):
     custom_grid = CustomGrid(xyx_coords=xyz_coord)
     grid.custom_grid = custom_grid
 
-    set_active_grid(grid, [GridTypes.CUSTOM])
+    set_active_grid(grid, [Grid.GridTypes.CUSTOM])
     return grid.custom_grid
 
 
@@ -79,7 +78,7 @@ def set_centered_grid(grid: Grid, centers: np.ndarray, resolution: Sequence[floa
         radius=radius
     )
     grid.centered_grid = centered_grid
-    set_active_grid(grid, [GridTypes.CENTERED])
+    set_active_grid(grid, [Grid.GridTypes.CENTERED])
     return grid.centered_grid
 
 
@@ -93,7 +92,7 @@ def set_topography_from_array():
 
 def set_active_grid(grid: Grid, grid_type: list[Grid.GridTypes], reset: bool = False):
     if reset is True:
-        grid.active_grids = GridTypes.NONE
+        grid.active_grids = Grid.GridTypes.NONE
     for grid_type in grid_type:
         grid.active_grids |= grid_type
 

gempy/API/initialization_API.py

@@ -6,6 +6,7 @@
 
 from gempy.API.io_API import read_surface_points, read_orientations
 from gempy_engine.core.data import InterpolationOptions
+from ..core.data.grid_modules import RegularGrid
 from ..optional_dependencies import require_subsurface
 from ..core.data import StructuralElement
 from ..core.data.geo_model import GeoModel
@@ -44,10 +45,12 @@ def create_geomodel(
     """
     # init resolutions well
     if resolution is None:
-        grid: Grid = Grid(extent=extent)
-        grid.octree_levels = refinement
+        grid: Grid = Grid.init_octree_grid(
+            extent=extent,
+            octree_levels=refinement
+        )
     else:
-        grid = Grid(
+        grid: Grid = Grid.init_dense_grid(
             extent=extent,
             resolution=resolution
         )
@@ -92,18 +95,18 @@ def structural_elements_from_borehole_set(
     Returns:
         list[StructuralElement]: A list of StructuralElement instances.
     """
-    
+
     ss = require_subsurface()
     borehole_set: ss.core.geological_formats.BoreholeSet
-    
+
     elements = []
     component_lith: dict[Hashable, np.ndarray] = borehole_set.get_bottom_coords_for_each_lith()
-    
+
     for name, properties in elements_dict.items():
         top_coordinates = component_lith.get(properties['id'])
         if top_coordinates is None:
             raise ValueError(f"Top lithology {properties['id']} not found in borehole set.")
-        
+
         element = StructuralElement(
             name=name,
             id=properties['id'],
@@ -119,7 +122,7 @@ def structural_elements_from_borehole_set(
         )
         elements.append(element)
     # Reverse the list to have the oldest rocks at the bottom
-    
+
     return elements
 
 

gempy/core/data/geo_model.py

@@ -233,7 +233,7 @@ def project_bounds(self) -> np.ndarray:
         return self.grid.bounding_box
 
     @property
-    def extent_transformed(self) -> np.ndarray:
+    def extent_transformed_transformed_by_input(self) -> np.ndarray:
         transformed = self.input_transform.apply(self.project_bounds)  # ! grid already has the grid transform applied
         new_extents = np.array([transformed[:, 0].min(), transformed[:, 0].max(),
                                 transformed[:, 1].min(), transformed[:, 1].max(),