Merge remote-tracking branch 'origin/main'

Author: javoha

.gitignore

@@ -157,3 +157,9 @@ examples/integrations/*.out
 # Approval tests
 *.received.txt
 *.Identifier
+/.obsidian/app.json
+/.obsidian/appearance.json
+/.obsidian/core-plugins.json
+/.obsidian/core-plugins-migration.json
+/.obsidian/hotkeys.json
+/.obsidian/workspace.json

DevelopersGuide.md

@@ -27,15 +27,6 @@ New
     twine upload dist/*
 ```
 
-Old
-```
-    #  First create the dist
-    python3 setup.py sdist bdist_wheel
-
-    # Second upload the distributions
-    twine upload dist/*
-
-
 ### Type of commits:
 
 - ENH: Enhancement, new functionality

gempy/NOTES.md renamed to docs/developers_notes/dev_log/2022-2023.md

@@ -0,0 +1,93 @@
+## Optimizing solver
+
+### Debugging Vector model 1
+
+#### Questions
+
+- [x] How many times do the system of equations gets solved?
+    - Without caching is 1 per scalar field x 2 (centers and corners) x n_octree levels
+
+- [x] How many iterations we need to reach the solution?
+  - Model 1 takes 5000 iterations (actually starts to diverge after 4k)
+   
+- [ ] Has Pykeops become slower with the latest additions to the gradient?
+
+- [x] Is the Jacobi pre-conditioner bad for all cases?
+  - Model 1:
+    - Numpy Jacobi: Converge in
+      - 1602 iterations and reach 5000
+      - time 150 sec
+      - More or less good results: Mesh closed
+    - Numpy No jacobi
+      - 1018 iterations and reach 5000
+      - time 150 sec too
+      - Results similar? I would say maybe a bit worse
+    - Pykeosp No jacobi
+      - 92 iterations and reach 5000
+      - time 324 (all on running the solver)
+- [ ] Can be that pykeops is faster than invert the matrix but still slower that any other CG
+  - Since the matrix is already constructed, it makes sense
+    - [ ] Then which cases is worth using pykeops? The number of iterations will need to be rather small
+     Will it make sense on edits reusing the previous weights? This is going to need a hell of orchestrator
+- [ ] How much memory takes numpy direct, numpy cg and pykeops cg
+  - numpy direct: 1.5 or so
+  - numpy cg: 1.25
+  - pykeops cg: 1.1
+  - **Note** Now the exporting takes quite a lot of memory
+
+
+
+### Experiments conclusions
+- Using pykeops CG to solve a ill posed problem is not realistic 
+  - if the model fits in memory probably the best thing to do is make the first run in torch an use the weights as initial point
+  - if the model does not fit in memory we should try to construct the sparse version of it and try to run the incomplete cholesky decomposition
+
+- I didn't find any improvement using the Jacobi pre-conditioner. Mostly the opposite
+- With ichol we get different results with numpy or pykeops
+	- Pykeops
+- ichol does not give super good results either at least before the convergence happens
+  - Model 1 converge **without** ichol on 1290 and 8518
+  - Model 1 converge **With** ichol 496 (506 in pykeops) 5080 (4995 in pykeops)
+	
+
+### TODO:
+
+- [x] Clean solver interface
+- [ ] Add citation for pymatting
+
+### Things to try:
+- [x] Pykeops as linear operator to use exact solver
+  - Does not look good. Probably is not even possible mathematically
+- [x] Falkon
+
+### Optimal solver
+
+- Tools we have:
+  - keops
+  - CG
+  - Incomplete cholesky decomposition
+  - passing the weights as initial point
+- Variables to choose what to use:
+  - Size of the covariance matrix
+  - Number of CG iterations we are going to need. This depends on:
+    - Conditional number
+    - Do we have weights
+
+- **Strategy**:
+  - If the Number of CG iterations is a handful -> **use keops**
+  - If the number of CG iterations is large:
+    - Use **direct solver** if the matrix fits in memory
+    - Then for changes hopefully we can use the weights as initial point and use **keops**
+  - If the matrix does not fit in memory:
+    - Construct the sparse version of the matrix
+    - Use **incomplete cholesky decomposition**
+      - Once we have the decomposition is better to use **keops** to solve the system (except if the number of iterations has to be very large) 
+  - Optimizing nugget effects still seems the best way to do this
+  - If memory is not an issue. **Direct solver with torch and keops for exporting seems to be the best option**
+
+
+#### Notes
+- [ ] I need to make similar default env file for gempy itself and using the default when we call the engine or whatever
+
+
+

docs/developers_notes/dev_log/log_2024-05.md

@@ -0,0 +1,7 @@
+# Refactoring grid to add transform
+
+## TODO:
+- [ ] Rotation of pyvista grid
+- [ ] Rotation of pyvista meshes
+ 
+## Questions:

docs/developers_notes/dev_log/log_2024-06.md

@@ -8,6 +8,7 @@
 
 import gempy as gp
 from gempy.core.data import Grid
+from gempy.core.data.grid_modules import RegularGrid
 
 np.random.seed(55500)
 
@@ -103,11 +104,15 @@
 # activate them.
 # 
 
+
 # %% 
-help(grid.create_regular_grid)
+help(RegularGrid)
 
 # %% 
-grid.create_regular_grid(extent=[0, 100, 0, 100, -100, 0], resolution=[20, 20, 20])
+extent = np.array([0, 100, 0, 100, -100, 0])
+resolution = np.array([20, 20, 20])
+grid.regular_grid = RegularGrid(extent, resolution)
+grid.regular_grid
 
 # %%
 # Now the regular grid object composed on ``Grid`` has been filled:

examples/tutorials/ch1_fundamentals/ch1_3a_grids.py

@@ -47,7 +47,7 @@
 geo_model.structural_frame
 
 # %%
-geo_model.transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
+geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
 gp.add_structural_group(
     model=geo_model,
     group_index=0,

examples/tutorials/ch1_fundamentals/ch1_4_onlap_relations.py

@@ -63,7 +63,7 @@
 
 # %% 
 
-geo_model.transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
+geo_model.input_transform.apply_anisotropy(gp.data.GlobalAnisotropy.NONE)
 if False:
     gp.compute_model(geo_model)
     # %%

examples/tutorials/ch1_fundamentals/ch1_5_fault_relations.py

@@ -93,7 +93,7 @@
 # creating a domain object from the gempy solution, a defined domain conditioning data
 domain = kriging.Domain(
     model_solutions=sol,
-    transform=geo_data.transform,
+    transform=geo_data.input_transform,
     domain=[2],
     data=cond_data
 )

examples/tutorials/ch3-Interpolations/ch3_1_kriging_interpolation_and_simulation.py

@@ -2,17 +2,17 @@
 
 import numpy as np
 
-from gempy.core.data import Grid
-from gempy.core.data.grid import GridTypes
-from gempy.core.data.grid_modules import CustomGrid
-from gempy.core.data.grid_modules.topography import Topography
-from gempy.modules.grids.create_topography import create_random_topography
-from gempy.optional_dependencies import require_subsurface
+from ..core.data import Grid
+from ..core.data.grid_modules import CustomGrid, Sections
+from ..core.data.grid_modules.topography import Topography
+from ..modules.grids.create_topography import create_random_topography
+from ..optional_dependencies import require_subsurface
 
 
 def set_section_grid(grid: Grid, section_dict: dict):
     if grid.sections is None:
-        grid.create_section_grid(section_dict=section_dict)
+        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)
@@ -91,12 +91,12 @@ def set_topography_from_array():
     raise NotImplementedError("This is not implemented yet")
 
 
-def set_active_grid(grid: Grid, grid_type: list[GridTypes], reset: bool = False):
+def set_active_grid(grid: Grid, grid_type: list[Grid.GridTypes], reset: bool = False):
     if reset is True:
-        grid.deactivate_all_grids()
+        grid.active_grids = GridTypes.NONE
     for grid_type in grid_type:
-        grid.active_grids_bool[grid_type.value] = True
+        grid.active_grids |= grid_type
 
-    print(f'Active grids: {grid.grid_types[grid.active_grids_bool]}')
+    print(f'Active grids: {grid.active_grids}')
 
     return grid