Source code for gempy.API.compute_API

from typing import Optional

import numpy as np

import gempy_engine
from gempy_engine.core.backend_tensor import BackendTensor
from gempy.API.gp2_gp3_compatibility.gp3_to_gp2_input import gempy3_to_gempy2
from gempy_engine.config import AvailableBackends
from gempy_engine.core.data import Solutions
from gempy_engine.core.data.interpolation_input import InterpolationInput
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


[docs] def compute_model(gempy_model: GeoModel, engine_config: Optional[GemPyEngineConfig] = None) -> Solutions: """ Compute the geological model given the provided GemPy model. Args: gempy_model (GeoModel): The GemPy model to compute. engine_config (Optional[GemPyEngineConfig]): Configuration for the computational engine. Defaults to None, in which case a default configuration will be used. Raises: ValueError: If the provided backend in the engine_config is not supported. Returns: Solutions: The computed geological model. """ engine_config = engine_config or GemPyEngineConfig(use_gpu=False) match engine_config.backend: case AvailableBackends.numpy | AvailableBackends.PYTORCH: BackendTensor.change_backend_gempy( engine_backend=engine_config.backend, use_gpu=engine_config.use_gpu, dtype=engine_config.dtype ) # TODO: To decide what to do with this. 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( interpolation_input=interpolation_input, options=gempy_model.interpolation_options, data_descriptor=gempy_model.input_data_descriptor, geophysics_input=gempy_model.geophysics_input, ) case AvailableBackends.aesara | AvailableBackends.legacy: gempy_model.legacy_model = _legacy_compute_model(gempy_model) case _: raise ValueError(f'Backend {engine_config} not supported') return gempy_model.solutions
def compute_model_at(gempy_model: GeoModel, at: np.ndarray, engine_config: Optional[GemPyEngineConfig] = None) -> np.ndarray: """ Compute the geological model at specific coordinates. Note: This function sets a custom grid and computes the model so be wary of side effects. Args: gempy_model (GeoModel): The GemPy model to compute. at (np.ndarray): The coordinates at which to compute the model. engine_config (Optional[GemPyEngineConfig], optional): Configuration for the computational engine. Defaults to None, in which case a default configuration will be used. Returns: np.ndarray: The computed geological model at the specified coordinates. """ set_custom_grid( grid=gempy_model.grid, xyz_coord=at ) sol = compute_model(gempy_model, engine_config) return sol.raw_arrays.custom def optimize_and_compute(geo_model: GeoModel, engine_config: GemPyEngineConfig, max_epochs: int = 10, convergence_criteria: float = 1e5): if engine_config.backend != AvailableBackends.PYTORCH: raise ValueError(f'Only PyTorch backend is supported for optimization. Received {engine_config.backend}') BackendTensor.change_backend_gempy( engine_backend=engine_config.backend, use_gpu=engine_config.use_gpu, dtype=engine_config.dtype ) import torch from gempy_engine.core.data.continue_epoch import ContinueEpoch interpolation_input: InterpolationInput = interpolation_input_from_structural_frame(geo_model) geo_model.taped_interpolation_input = interpolation_input nugget_effect_scalar: torch.Tensor = geo_model.taped_interpolation_input.surface_points.nugget_effect_scalar optimizer = torch.optim.Adam( params=[nugget_effect_scalar], lr=0.01, ) # Optimization loop geo_model.interpolation_options.kernel_options.optimizing_condition_number = True def _check_convergence_criterion(conditional_number: float, condition_number_old: float, conditional_number_target: float = 1e5): reached_conditional_target = conditional_number < conditional_number_target if reached_conditional_target == False and epoch > 10: condition_number_change = torch.abs(conditional_number - condition_number_old) / condition_number_old if condition_number_change < 0.01: reached_conditional_target = True return reached_conditional_target previous_condition_number = 0 for epoch in range(max_epochs): optimizer.zero_grad() try: # geo_model.taped_interpolation_input.grid = geo_model.interpolation_input_copy.grid gempy_engine.compute_model( interpolation_input=geo_model.taped_interpolation_input, options=geo_model.interpolation_options, data_descriptor=geo_model.input_data_descriptor, geophysics_input=geo_model.geophysics_input, ) except ContinueEpoch: # Get absolute values of gradients grad_magnitudes = torch.abs(nugget_effect_scalar.grad) # Get indices of the 10 largest gradients grad_magnitudes.size # * This ignores 90 percent of the gradients # To int n_values = int(grad_magnitudes.size()[0] * 0.9) _, indices = torch.topk(grad_magnitudes, n_values, largest=False) # Zero out gradients that are not in the top 10 mask = torch.ones_like(nugget_effect_scalar.grad) mask[indices] = 0 nugget_effect_scalar.grad *= mask # Update the vector optimizer.step() nugget_effect_scalar.data = nugget_effect_scalar.data.clamp_(min=1e-7) # Replace negative values with 0 # optimizer.zero_grad() # Monitor progress if epoch % 1 == 0: # print(f"Epoch {epoch}: Condition Number = {condition_number.item()}") print(f"Epoch {epoch}") if _check_convergence_criterion( conditional_number=geo_model.interpolation_options.kernel_options.condition_number, condition_number_old=previous_condition_number, conditional_number_target=convergence_criteria, ): break previous_condition_number = geo_model.interpolation_options.kernel_options.condition_number continue geo_model.interpolation_options.kernel_options.optimizing_condition_number = False geo_model.solutions = gempy_engine.compute_model( interpolation_input=geo_model.taped_interpolation_input, options=geo_model.interpolation_options, data_descriptor=geo_model.input_data_descriptor, geophysics_input=geo_model.geophysics_input, ) return geo_model.solutions def _legacy_compute_model(gempy_model: GeoModel) -> 'gempy_legacy.Project': gpl = require_gempy_legacy() legacy_model: gpl.Project = gempy3_to_gempy2(gempy_model) gpl.set_interpolator(legacy_model) gpl.compute_model(legacy_model) return legacy_model