From 48accccde1975d9cc813d4fbedf287c6067b929a Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 3 Mar 2026 16:03:34 -0800 Subject: [PATCH 01/24] Begin refactoring --- simpeg_drivers/driver.py | 13 +- simpeg_drivers/electricals/base_2d.py | 296 ++++++++++++++++++ .../pseudo_three_dimensions/driver.py | 36 --- .../pseudo_three_dimensions/options.py | 92 ------ .../direct_current/two_dimensions/driver.py | 2 +- .../direct_current/two_dimensions/options.py | 4 +- .../pseudo_three_dimensions/driver.py | 40 --- .../pseudo_three_dimensions/options.py | 92 ------ .../two_dimensions/driver.py | 2 +- simpeg_drivers/utils/surveys.py | 2 +- simpeg_drivers/utils/synthetics/driver.py | 4 +- .../{meshes/octrees.py => meshes.py} | 36 ++- .../utils/synthetics/meshes/__init__.py | 9 - .../utils/synthetics/meshes/factory.py | 44 --- .../utils/synthetics/meshes/tensors.py | 55 ---- simpeg_drivers/utils/synthetics/options.py | 8 +- simpeg_drivers/utils/synthetics/topography.py | 12 +- simpeg_drivers/utils/utils.py | 57 +--- tests/run_tests/driver_dc_b2d_test.py | 40 ++- 19 files changed, 394 insertions(+), 450 deletions(-) create mode 100644 simpeg_drivers/electricals/base_2d.py delete mode 100644 simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/driver.py delete mode 100644 simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/options.py delete mode 100644 simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/driver.py delete mode 100644 simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/options.py rename simpeg_drivers/utils/synthetics/{meshes/octrees.py => meshes.py} (76%) delete mode 100644 simpeg_drivers/utils/synthetics/meshes/__init__.py delete mode 100644 simpeg_drivers/utils/synthetics/meshes/factory.py delete mode 100644 simpeg_drivers/utils/synthetics/meshes/tensors.py diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index fc5a24e1..05d75bed 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -823,7 +823,18 @@ def get_tiles(self) -> dict[str, list[np.ndarray]]: tiles = [[tile] for tile in np.array_split(indices, n_chunks)] elif "2d" in self.params.inversion_type: - tiles = [[indices]] + tiles = [ + [ + [ + np.where( + self.params.line_selection.line_object.values == line_id + )[0] + ] + for line_id in np.unique( + self.params.line_selection.line_object.values + ) + ] + ] else: tiles = tile_locations( diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py new file mode 100644 index 00000000..7eb6f41d --- /dev/null +++ b/simpeg_drivers/electricals/base_2d.py @@ -0,0 +1,296 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + + +from __future__ import annotations + +import numpy as np +from geoapps_utils.utils.locations import get_locations +from geoapps_utils.utils.numerical import weighted_average +from geoh5py.data import Data, IntegerData +from geoh5py.groups import PropertyGroup +from geoh5py.objects import DrapeModel +from geoh5py.shared.merging.drape_model import DrapeModelMerger +from geoh5py.ui_json.ui_json import fetch_active_workspace +from geoh5py.workspace import Workspace + +from simpeg_drivers.components.data import InversionData +from simpeg_drivers.components.meshes import InversionMesh +from simpeg_drivers.components.topography import InversionTopography +from simpeg_drivers.components.windows import InversionWindow +from simpeg_drivers.driver import InversionDriver +from simpeg_drivers.line_sweep.driver import LineSweepDriver +from simpeg_drivers.options import ( + BaseForwardOptions, + BaseInversionOptions, + DrapeModelOptions, + LineSelectionOptions, +) +from simpeg_drivers.utils.surveys import extract_dcip_survey +from simpeg_drivers.utils.utils import get_drape_model + + +class Base2DDriver(InversionDriver): + """ + Base class for 2D DC and IP forward and inversion drivers. + + Survey lines are inverted independently and internally stacked as a single + long survey. The inversion mesh is created as a drape mesh over the survey lines. + """ + + @property + def inversion_mesh(self) -> InversionMesh: + """Inversion mesh""" + if getattr(self, "_inversion_mesh", None) is None: + with fetch_active_workspace(self.workspace, mode="r+"): + entity = self.params.mesh + if entity is None: + entity = create_mesh_by_line_id( + self.workspace, + self.params.line_selection.line_object, + self.params.drape_model, + parent=self.out_group, + ) + + self._inversion_mesh = InversionMesh( + self.workspace, self.params, entity=entity + ) + + return self._inversion_mesh + + +def create_mesh_by_line_id( + workspace: Workspace, + line_ids: IntegerData, + drape_options: DrapeModelOptions, + **object_kwargs, +) -> DrapeModel: + """ + Create a drape mesh for the dc resistivity survey lines. + + :param workspace: Workspace to create the drape mesh in. + :param line_ids: IntegerData object containing the line IDs for each vertex. + :param drape_options: DrapeModelOptions containing the parameters for the drape mesh + :param object_kwargs: Additional keyword arguments to pass to the DrapeModelMerger.create_object method. + + :return: A DrapeModel object containing the merged drape mesh for all survey lines. + """ + drape_models = [] + temp_work = Workspace() + for line_id in np.unique(line_ids.values): + poles = get_poles_by_line_id(line_ids, line_id) + drape_model = get_drape_model( + temp_work, + poles, + [ + drape_options.u_cell_size, + drape_options.v_cell_size, + ], + drape_options.depth_core, + [0.0] * 2 + [drape_options.vertical_padding, 1], + drape_options.expansion_factor, + ) + drape_models.append(drape_model) + + entity = DrapeModelMerger.create_object(workspace, drape_models, **object_kwargs) + + return entity + + +def get_poles_by_line_id(line_ids: IntegerData, uid: int) -> np.ndarray: + """Get the vertices associated with a given line ID.""" + mn_mask = line_ids.values == uid + + unique_tx = np.unique(line_ids.parent.ab_cell_id.values[mn_mask]) + + ab_mask = np.isin(line_ids.parent.complement.ab_cell_id.values, unique_tx) + + return np.vstack( + [ + line_ids.parent.vertices[line_ids.parent.cells[mn_mask].flatten()], + line_ids.parent.current_electrodes.vertices[ + line_ids.parent.current_electrodes.cells[ab_mask].flatten() + ], + ] + ) + + +# class BaseBatch2DDriver(LineSweepDriver): +# """Base class for batch 2D DC and IP forward and inversion drivers.""" +# +# _params_class: type[BaseForwardOptions | BaseInversionOptions] +# _params_2d_class: type[BaseForwardOptions | BaseInversionOptions] +# +# _model_list: list[str] = [] +# +# def __init__(self, params): +# super().__init__(params) +# if params.file_control.files_only: +# sys.exit("Files written") +# +# def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID | float]: +# """ +# Transfer models from the input parameters to the output drape mesh. +# +# :param mesh: Destination DrapeModel object. +# """ +# models = {"starting_model": self.batch2d_params.models.starting_model} +# +# for model in self._model_list: +# models[model] = getattr(self.batch2d_params, model, None) +# +# if not self.batch2d_params.forward_only: +# for model in ["reference_model", "lower_bound", "upper_bound"]: +# models[model] = getattr(self.batch2d_params.models, model) +# +# if self.batch2d_params.models.gradient_rotation is not None: +# group_properties = {} +# for prop in self.batch2d_params.models.gradient_rotation.properties: +# model = self.batch2d_params.mesh.get_data(prop)[0] +# group_properties[model.name] = model +# +# models.update(group_properties) +# +# if self.batch2d_params.mesh is not None: +# xyz_in = get_locations(self.workspace, self.batch2d_params.mesh) +# xyz_out = mesh.centroids +# +# for name, model in models.items(): +# if model is None: +# continue +# elif isinstance(model, Data): +# model_values = weighted_average( +# xyz_in, xyz_out, [model.values], n=1 +# )[0] +# else: +# model_values = model * np.ones(len(xyz_out)) +# +# model_object = mesh.add_data({name: {"values": model_values}}) +# models[name] = model_object +# +# if ( +# not self.batch2d_params.forward_only +# and self.batch2d_params.models.gradient_rotation is not None +# ): +# pg = PropertyGroup( +# mesh, +# properties=[models[prop] for prop in group_properties], +# property_group_type=self.batch2d_params.models.gradient_rotation.property_group_type, +# ) +# models["gradient_rotation"] = pg +# del models["azimuth"] +# del models["dip"] +# +# return models +# +# def write_files(self, lookup): +# """Write ui.geoh5 and ui.json files for sweep trials.""" +# +# kwargs_2d = {} +# with fetch_active_workspace(self.workspace, mode="r+"): +# for uid, trial in lookup.items(): +# if trial["status"] != "pending": +# continue +# +# filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" +# +# if filepath.exists(): +# warnings.warn( +# f"File {filepath} already exists but 'status' marked as 'pending'. " +# "Over-writing file." +# ) +# filepath.unlink() +# +# with Workspace.create(filepath) as iter_workspace: +# cell_mask: np.ndarray = ( +# self.batch2d_params.line_selection.line_object.values +# == trial["line_id"] +# ) +# +# if not np.any(cell_mask): +# continue +# +# receiver_entity = extract_dcip_survey( +# iter_workspace, self.inversion_data.entity, cell_mask +# ) +# current_entity = receiver_entity.current_electrodes +# receiver_locs = np.vstack( +# [receiver_entity.vertices, current_entity.vertices] +# ) +# +# mesh = get_drape_model( +# iter_workspace, +# "Models", +# receiver_locs, +# [ +# self.batch2d_params.drape_model.u_cell_size, +# self.batch2d_params.drape_model.v_cell_size, +# ], +# self.batch2d_params.drape_model.depth_core, +# [self.batch2d_params.drape_model.horizontal_padding] * 2 +# + [self.batch2d_params.drape_model.vertical_padding, 1], +# self.batch2d_params.drape_model.expansion_factor, +# )[0] +# +# model_parameters = self.transfer_models(mesh) +# +# for key in self._params_2d_class.model_fields: +# param = getattr(self.batch2d_params, key, None) +# if key not in ["title", "inversion_type"]: +# kwargs_2d[key] = param +# +# self.batch2d_params.active_cells.topography_object.copy( +# parent=iter_workspace, copy_children=True +# ) +# +# kwargs_2d.update( +# dict( +# **{ +# "geoh5": iter_workspace, +# "mesh": mesh, +# "data_object": receiver_entity, +# "line_selection": LineSelectionOptions( +# line_object=receiver_entity.get_data( +# self.batch2d_params.line_selection.line_object.name +# )[0], +# line_id=trial["line_id"], +# ), +# "out_group": None, +# }, +# **model_parameters, +# ) +# ) +# +# params = self._params_2d_class(**kwargs_2d) +# params.write_ui_json(Path(self.working_directory) / f"{uid}.ui.json") +# +# lookup[uid]["status"] = "written" +# +# _ = self.update_lookup(lookup) # pylint: disable=no-member +# +# @property +# def inversion_data(self) -> InversionData: +# """Inversion data""" +# if getattr(self, "_inversion_data", None) is None: +# with fetch_active_workspace(self.workspace, mode="r+"): +# self._inversion_data = InversionData( +# self.workspace, self.batch2d_params +# ) +# +# return self._inversion_data +# +# @property +# def inversion_topography(self): +# """Inversion topography""" +# if getattr(self, "_inversion_topography", None) is None: +# self._inversion_topography = InversionTopography( +# self.workspace, self.batch2d_params +# ) +# return self._inversion_topography diff --git a/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/driver.py b/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/driver.py deleted file mode 100644 index 69d80ead..00000000 --- a/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/driver.py +++ /dev/null @@ -1,36 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.options import ( - DCBatch2DForwardOptions, - DCBatch2DInversionOptions, -) -from simpeg_drivers.electricals.direct_current.two_dimensions.options import ( - DC2DForwardOptions, - DC2DInversionOptions, -) -from simpeg_drivers.electricals.driver import BaseBatch2DDriver - - -class DCBatch2DForwardDriver(BaseBatch2DDriver): - """Direct Current batch 2D forward driver.""" - - _params_class = DCBatch2DForwardOptions - _params_2d_class = DC2DForwardOptions - - -class DCBatch2DInversionDriver(BaseBatch2DDriver): - """Direct Current batch 2D inversion driver.""" - - _params_class = DCBatch2DInversionOptions - _params_2d_class = DC2DInversionOptions diff --git a/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/options.py b/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/options.py deleted file mode 100644 index 59083593..00000000 --- a/simpeg_drivers/electricals/direct_current/pseudo_three_dimensions/options.py +++ /dev/null @@ -1,92 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from pathlib import Path -from typing import ClassVar - -from geoh5py.data import FloatData -from geoh5py.objects import Octree, PotentialElectrode - -from simpeg_drivers import assets_path -from simpeg_drivers.electricals.options import ( - FileControlOptions, -) -from simpeg_drivers.options import ( - BaseForwardOptions, - BaseInversionOptions, - ConductivityModelOptions, - DrapeModelOptions, - LineSelectionOptions, -) - - -class DCBatch2DForwardOptions(BaseForwardOptions): - """ - Direct Current batch 2D forward options. - - :param data_object: DC survey object. - :param potential_channel_bool: Potential channel boolean. - :param line_selection: Line selection parameters. - :param mesh: Optional mesh object if providing a heterogeneous model. - :param drape_model: Drape model parameters common to all 2D simulations. - :param file_control: File control parameters. - """ - - name: ClassVar[str] = "Direct Current Pseudo 3D Forward" - default_ui_json: ClassVar[Path] = ( - assets_path() / "uijson/direct_current_batch2d_forward.ui.json" - ) - - title: str = "Direct Current (DC) 2D Batch Forward" - physical_property: str = "conductivity" - inversion_type: str = "direct current pseudo 3d" - - data_object: PotentialElectrode - potential_channel_bool: bool = True - line_selection: LineSelectionOptions - mesh: Octree | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() - file_control: FileControlOptions = FileControlOptions() - models: ConductivityModelOptions - - -class DCBatch2DInversionOptions(BaseInversionOptions): - """ - Direct Current batch 2D Inversion options. - - :param data_object: DC survey object. - :param potential_channel: Potential data channel. - :param potential_uncertainty: Potential data uncertainty channel. - :param line_selection: Line selection parameters. - :param mesh: Optional mesh object if providing a heterogeneous model. - :param drape_model: Drape model parameters. - :param file_control: File control parameters. - """ - - name: ClassVar[str] = "Direct Current Pseudo 3D Inversion" - default_ui_json: ClassVar[Path] = ( - assets_path() / "uijson/direct_current_batch2d_inversion.ui.json" - ) - - title: str = "Direct Current (DC) 2D Batch Inversion" - physical_property: str = "conductivity" - inversion_type: str = "direct current pseudo 3d" - - data_object: PotentialElectrode - potential_channel: FloatData - potential_uncertainty: float | FloatData - line_selection: LineSelectionOptions - mesh: Octree | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() - file_control: FileControlOptions = FileControlOptions() - models: ConductivityModelOptions diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py index 0f5b2417..5642d81a 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py @@ -11,7 +11,7 @@ from __future__ import annotations -from simpeg_drivers.electricals.driver import Base2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver from .options import DC2DForwardOptions, DC2DInversionOptions diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py index 53425542..89292cb8 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py @@ -49,7 +49,7 @@ class DC2DForwardOptions(BaseForwardOptions): potential_channel_bool: bool = True line_selection: LineSelectionOptions mesh: DrapeModel | None = None - drape_model: DrapeModelOptions + drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions @@ -77,5 +77,5 @@ class DC2DInversionOptions(BaseInversionOptions): potential_uncertainty: float | FloatData | None = None line_selection: LineSelectionOptions mesh: DrapeModel | None = None - drape_model: DrapeModelOptions + drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions diff --git a/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/driver.py b/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/driver.py deleted file mode 100644 index 394586e8..00000000 --- a/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/driver.py +++ /dev/null @@ -1,40 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from simpeg_drivers.electricals.driver import BaseBatch2DDriver -from simpeg_drivers.electricals.induced_polarization.pseudo_three_dimensions.options import ( - IPBatch2DForwardOptions, - IPBatch2DInversionOptions, -) -from simpeg_drivers.electricals.induced_polarization.two_dimensions.options import ( - IP2DForwardOptions, - IP2DInversionOptions, -) - - -class IPBatch2DForwardDriver(BaseBatch2DDriver): - """Induced Polarization batch 2D forward driver.""" - - _params_class = IPBatch2DForwardOptions - _params_2d_class = IP2DForwardOptions - - _model_list = ["conductivity_model"] - - -class IPBatch2DInversionDriver(BaseBatch2DDriver): - """Induced Polarization batch 2D inversion driver.""" - - _params_class = IPBatch2DInversionOptions - _params_2d_class = IP2DInversionOptions - - _model_list = ["conductivity_model"] diff --git a/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/options.py b/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/options.py deleted file mode 100644 index 726caa2f..00000000 --- a/simpeg_drivers/electricals/induced_polarization/pseudo_three_dimensions/options.py +++ /dev/null @@ -1,92 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from pathlib import Path -from typing import ClassVar - -from geoh5py.data import FloatData -from geoh5py.objects import Octree, PotentialElectrode - -from simpeg_drivers import assets_path -from simpeg_drivers.electricals.options import ( - FileControlOptions, - IPModelOptions, -) -from simpeg_drivers.options import ( - BaseForwardOptions, - BaseInversionOptions, - DrapeModelOptions, - LineSelectionOptions, -) - - -class IPBatch2DForwardOptions(BaseForwardOptions): - """ - Induced Polarization batch 2D forward options. - - :param data_object: DC/IP survey object. - :param chargeability_channel_bool: Chargeability channel boolean. - :param line_selection: Line selection parameters. - :param mesh: Optional mesh object if providing a heterogeneous model. - :param drape_model: Drape model parameters common to all 2D simulations. - :param file_control: File control parameters. - """ - - name: ClassVar[str] = "Induced Polarization Pseudo 3D Forward" - default_ui_json: ClassVar[Path] = ( - assets_path() / "uijson/induced_polarization_batch2d_forward.ui.json" - ) - - title: str = "Induced Polarization (IP) 2D Batch Forward" - physical_property: str = "chargeability" - inversion_type: str = "induced polarization pseudo 3d" - - data_object: PotentialElectrode - chargeability_channel_bool: bool = True - line_selection: LineSelectionOptions - mesh: Octree | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() - file_control: FileControlOptions = FileControlOptions() - models: IPModelOptions - - -class IPBatch2DInversionOptions(BaseInversionOptions): - """ - Induced Polarization batch 2D inversion options. - - :param data_object: DC/IP survey object. - :param chargeability_channel: Chargeability data channel. - :param chargeability_uncertainty: Chargeability data uncertainty channel. - :param line_selection: Line selection parameters. - :param mesh: Optional mesh object if providing a heterogeneous model. - :param drape_model: Drape model parameters common to all 2D simulations. - :param file_control: File control parameters. - """ - - name: ClassVar[str] = "Induced Polarization Pseudo 3D Inversion" - default_ui_json: ClassVar[Path] = ( - assets_path() / "uijson/induced_polarization_batch2d_inversion.ui.json" - ) - - title: str = "Induced Polarization (IP) 2D Batch Inversion" - physical_property: str = "chargeability" - inversion_type: str = "induced polarization pseudo 3d" - - data_object: PotentialElectrode - chargeability_channel: FloatData - chargeability_uncertainty: float | FloatData - line_selection: LineSelectionOptions - mesh: Octree | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() - file_control: FileControlOptions = FileControlOptions() - models: IPModelOptions diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py index 6b646606..b1322e92 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py @@ -11,7 +11,7 @@ from __future__ import annotations -from simpeg_drivers.electricals.driver import Base2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver from .options import ( IP2DForwardOptions, diff --git a/simpeg_drivers/utils/surveys.py b/simpeg_drivers/utils/surveys.py index 21b69a08..78da3ff8 100644 --- a/simpeg_drivers/utils/surveys.py +++ b/simpeg_drivers/utils/surveys.py @@ -76,7 +76,7 @@ def compute_alongline_distance(points: np.ndarray, ordered: bool = True): distances = np.cumsum( np.r_[0, np.linalg.norm(np.diff(points[:, :2], axis=0), axis=1)] ) - if points.shape[1] == 3: + if points.shape[1] > 2: distances = np.c_[distances, points[:, 2:]] return distances diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index 9b6c848f..ec921484 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -13,7 +13,7 @@ from geoh5py.data import BooleanData, FloatData from geoh5py.objects import DrapeModel, ObjectBase, Octree, Surface -from simpeg_drivers.utils.synthetics.meshes.factory import get_mesh +from simpeg_drivers.utils.synthetics.meshes import get_mesh from simpeg_drivers.utils.synthetics.models import get_model from simpeg_drivers.utils.synthetics.options import SyntheticsComponentsOptions from simpeg_drivers.utils.synthetics.surveys.factory import get_survey @@ -71,7 +71,7 @@ def survey(self): @property def mesh(self): - entity = self.geoh5.get_entity(self.options.mesh.name)[0] + entity = self.geoh5.get_entity("mesh")[0] assert isinstance(entity, Octree | DrapeModel | type(None)) if entity is None: assert self.options is not None diff --git a/simpeg_drivers/utils/synthetics/meshes/octrees.py b/simpeg_drivers/utils/synthetics/meshes.py similarity index 76% rename from simpeg_drivers/utils/synthetics/meshes/octrees.py rename to simpeg_drivers/utils/synthetics/meshes.py index 795cf952..db351883 100644 --- a/simpeg_drivers/utils/synthetics/meshes/octrees.py +++ b/simpeg_drivers/utils/synthetics/meshes.py @@ -9,12 +9,44 @@ # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' import numpy as np -from discretize import TreeMesh +from discretize import TensorMesh, TreeMesh from discretize.utils import mesh_builder_xyz -from geoh5py.objects import Octree, Points, Surface +from geoh5py.objects import CellObject, DrapeModel, Octree, Points, Surface from grid_apps.octree_creation.driver import OctreeDriver from grid_apps.utils import treemesh_2_octree +from simpeg_drivers.electricals.base_2d import create_mesh_by_line_id +from simpeg_drivers.options import DrapeModelOptions +from simpeg_drivers.utils.synthetics.options import MeshOptions + + +def get_mesh( + method: str, + survey: Points, + topography: Surface, + options: MeshOptions | DrapeModelOptions, +) -> DrapeModel | Octree: + """Factory for mesh creation with behaviour modified by the provided method.""" + + if "2d" in method: + line_data = survey.get_entity("line_ids")[0] + + return create_mesh_by_line_id( + survey.workspace, + line_data, + options, + name="mesh", + ) + + return get_octree_mesh( + survey=survey, + topography=topography, + cell_size=options.cell_size, + refinement=options.refinement, + padding_distance=options.padding_distance, + name=options.name, + ) + def get_base_octree( survey: Points, diff --git a/simpeg_drivers/utils/synthetics/meshes/__init__.py b/simpeg_drivers/utils/synthetics/meshes/__init__.py deleted file mode 100644 index df32b204..00000000 --- a/simpeg_drivers/utils/synthetics/meshes/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/meshes/factory.py b/simpeg_drivers/utils/synthetics/meshes/factory.py deleted file mode 100644 index 8ee89802..00000000 --- a/simpeg_drivers/utils/synthetics/meshes/factory.py +++ /dev/null @@ -1,44 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -from discretize import TensorMesh, TreeMesh -from geoh5py.objects import CellObject, DrapeModel, Octree, Points, Surface - -from simpeg_drivers.utils.synthetics.options import MeshOptions - -from .octrees import get_octree_mesh -from .tensors import get_tensor_mesh - - -def get_mesh( - method: str, - survey: Points, - topography: Surface, - options: MeshOptions, -) -> DrapeModel | Octree: - """Factory for mesh creation with behaviour modified by the provided method.""" - - if "2d" in method: - assert isinstance(survey, CellObject) - return get_tensor_mesh( - survey=survey, - cell_size=options.cell_size, - padding_distance=options.padding_distance, - name=options.name, - ) - - return get_octree_mesh( - survey=survey, - topography=topography, - cell_size=options.cell_size, - refinement=options.refinement, - padding_distance=options.padding_distance, - name=options.name, - ) diff --git a/simpeg_drivers/utils/synthetics/meshes/tensors.py b/simpeg_drivers/utils/synthetics/meshes/tensors.py deleted file mode 100644 index 39ff84f4..00000000 --- a/simpeg_drivers/utils/synthetics/meshes/tensors.py +++ /dev/null @@ -1,55 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -import numpy as np -from geoh5py.data import IntegerData -from geoh5py.objects import CellObject, DrapeModel - -from simpeg_drivers.utils.utils import get_drape_model - - -def get_tensor_mesh( - survey: CellObject, - cell_size: tuple[float, float, float], - padding_distance: float, - line_id: int = 101, - name: str = "mesh", -) -> DrapeModel: - """ - Generate a tensor mesh and the colocated DrapeModel. - - :param survey: Survey object with vertices that define the core of the - tensor mesh. - :param cell_size: Tuple defining the cell size in all directions. - :param padding_distance: Distance to pad the mesh in all directions. - :param line_id: Chooses line from the survey to define the drape model. - - :return entity: The DrapeModel object that shares cells with the discretize - tensor mesh and which stores the results of computations. - :return mesh: The discretize tensor mesh object for computations. - """ - - line_data = survey.get_entity("line_ids")[0] - assert isinstance(line_data, IntegerData) - lines = line_data.values - entity, _mesh, _ = get_drape_model( # pylint: disable=unbalanced-tuple-unpacking - survey.workspace, - name, - survey.vertices[np.unique(survey.cells[lines == line_id, :]), :], - [cell_size[0], cell_size[2]], - 100.0, - [padding_distance] * 2 + [padding_distance] * 2, - 1.1, - parent=None, - return_colocated_mesh=True, - return_sorting=True, - ) - - return entity diff --git a/simpeg_drivers/utils/synthetics/options.py b/simpeg_drivers/utils/synthetics/options.py index 87fd1da0..564f8beb 100644 --- a/simpeg_drivers/utils/synthetics/options.py +++ b/simpeg_drivers/utils/synthetics/options.py @@ -14,6 +14,8 @@ from geoapps_utils.utils.locations import gaussian from pydantic import BaseModel, ConfigDict +from simpeg_drivers.options import DrapeModelOptions + class SurveyOptions(BaseModel): center: tuple[float, float] = (0.0, 0.0) @@ -40,6 +42,10 @@ class MeshOptions(BaseModel): refinement: tuple = (4, 6) padding_distance: float = 100.0 name: str = "mesh" + depth_core: float = 100.0 + expansion_factor: float = 1.1 + horizontal_padding: float = 1000.0 + vertical_padding: float = 1000.0 class ModelOptions(BaseModel): @@ -63,6 +69,6 @@ class SyntheticsComponentsOptions(BaseModel): method: str = "gravity" survey: SurveyOptions = SurveyOptions() - mesh: MeshOptions = MeshOptions() + mesh: MeshOptions | DrapeModelOptions = MeshOptions() model: ModelOptions = ModelOptions() active: ActiveCellsOptions = ActiveCellsOptions() diff --git a/simpeg_drivers/utils/synthetics/topography.py b/simpeg_drivers/utils/synthetics/topography.py index 82444000..039e66f5 100644 --- a/simpeg_drivers/utils/synthetics/topography.py +++ b/simpeg_drivers/utils/synthetics/topography.py @@ -13,6 +13,7 @@ from geoh5py.objects import DrapeModel, Octree, Surface from scipy.spatial import Delaunay +from simpeg_drivers.options import DrapeModelOptions from simpeg_drivers.utils.synthetics.options import ( MeshOptions, SurveyOptions, @@ -58,7 +59,7 @@ def get_topography_surface( def compute_mesh_extents( - survey_options: SurveyOptions, mesh_options: MeshOptions + survey_options: SurveyOptions, mesh_options: MeshOptions | DrapeModelOptions ) -> tuple[float, float]: """ Estimates the extent of the mesh from survey and mesh options. @@ -71,8 +72,13 @@ def compute_mesh_extents( """ width = survey_options.width height = survey_options.height - cell_size = mesh_options.cell_size - padding = mesh_options.padding_distance + + if isinstance(mesh_options, DrapeModelOptions): + cell_size = (mesh_options.u_cell_size, mesh_options.u_cell_size) + padding = mesh_options.horizontal_padding + else: + cell_size = mesh_options.cell_size + padding = mesh_options.padding_distance def next_pow2_cells(span, cell_size, padding): return 2 ** np.ceil(np.log2((span + 2 * padding) / cell_size)) diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 37eafb8c..942c8578 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -345,33 +345,28 @@ def floating_active(mesh: TensorMesh | TreeMesh, active: np.ndarray): def get_drape_model( workspace: Workspace, - name: str, locations: np.ndarray, h: list, depth_core: float, pads: list, expansion_factor: float, - parent: Group | None = None, return_colocated_mesh: bool = False, - return_sorting: bool = False, -) -> tuple: + **object_kwargs, +) -> DrapeModel | tuple[DrapeModel, TensorMesh]: """ Create a BlockModel object from parameters. :param workspace: Workspace. - :param parent: Group to contain the result. - :param name: Block model name. :param locations: Location points. :param h: Cell size(s) for the core mesh. :param depth_core: Depth of core mesh below locs. :param pads: len(4) Padding distances [W, E, Down, Up] :param expansion_factor: Expansion factor for padding cells. :param return_colocated_mesh: If true return TensorMesh. - :param return_sorting: If true, return the indices required to map - values stored in the TensorMesh to the drape model. + :param object_kwargs: Extra arguments to pass to the DrapeModel.create() method. + :return object_out: Output block model. """ - locations = truncate_locs_depths(locations, depth_core) order = traveling_salesman(locations) # Smooth the locations @@ -387,10 +382,8 @@ def get_drape_model( ] ) distances = compute_alongline_distance(xy_smooth) - distances[:, -1] += locations[:, 2].max() - distances[:, -1].max() + h[1] x_interp = interp1d(distances[:, 0], xy_smooth[:, 0], fill_value="extrapolate") y_interp = interp1d(distances[:, 0], xy_smooth[:, 1], fill_value="extrapolate") - mesh = mesh_utils.mesh_builder_xyz( distances, h, @@ -407,29 +400,21 @@ def get_drape_model( locations_top = np.c_[ x_interp(mesh.cell_centers_x), y_interp(mesh.cell_centers_x), top ] - model = xyz_2_drape_model(workspace, locations_top, hz, name, parent) - val = [model] + drape_model = xyz_2_drape_model(workspace, locations_top, hz, **object_kwargs) + if return_colocated_mesh: - val.append(mesh) - if return_sorting: - sorting = np.arange(mesh.n_cells) - sorting = sorting.reshape(mesh.shape_cells[1], mesh.shape_cells[0], order="C") - sorting = sorting[::-1].T.flatten() - val.append(sorting) - return val + return drape_model, mesh + return drape_model -def xyz_2_drape_model( - workspace, locations, depths, name=None, parent=None -) -> DrapeModel: +def xyz_2_drape_model(workspace, locations, depths, **object_kwargs) -> DrapeModel: """ Convert a list of cell tops and layer depths to a DrapeModel object. :param workspace: Workspace object :param locations: n x 3 array of cell centers [x, y, z_top] :param depths: n x 1 array of layer depths - :param name: Name of the new DrapeModel object - :param parent: Parent group for the new DrapeModel object + :param object_kwargs: Additional keyword arguments to pass to DrapeModel.create() :returns: DrapeModel object """ @@ -450,9 +435,7 @@ def xyz_2_drape_model( prisms = np.vstack(prisms) layers = np.vstack(layers) - model = DrapeModel.create( - workspace, layers=layers, name=name, prisms=prisms, parent=parent - ) + model = DrapeModel.create(workspace, layers=layers, prisms=prisms, **object_kwargs) model.add_data( { "indices": { @@ -558,24 +541,6 @@ def active_from_xyz( return mask_under_horizon(locations, horizon=topo, triangulation=triangulation) -def truncate_locs_depths(locs: np.ndarray, depth_core: float) -> np.ndarray: - """ - Sets locations below core to core bottom. - - :param locs: Location points. - :param depth_core: Depth of core mesh below locs. - - :return locs: locs with depths truncated. - """ - zmax = locs[:, -1].max() # top of locs - below_core_ind = (zmax - locs[:, -1]) > depth_core - core_bottom_elev = zmax - depth_core - locs[below_core_ind, -1] = ( - core_bottom_elev # sets locations below core to core bottom - ) - return locs - - def get_neighbouring_cells(mesh: TreeMesh, indices: list | np.ndarray) -> tuple: """ Get the indices of neighbouring cells along a given axis for a given list of diff --git a/tests/run_tests/driver_dc_b2d_test.py b/tests/run_tests/driver_dc_b2d_test.py index 82998c3a..3a3822d2 100644 --- a/tests/run_tests/driver_dc_b2d_test.py +++ b/tests/run_tests/driver_dc_b2d_test.py @@ -17,13 +17,13 @@ from geoh5py.groups import SimPEGGroup from geoh5py.workspace import Workspace -from simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.driver import ( - DCBatch2DForwardDriver, - DCBatch2DInversionDriver, +from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( + DC2DForwardDriver, + DC2DInversionDriver, ) -from simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.options import ( - DCBatch2DForwardOptions, - DCBatch2DInversionOptions, +from simpeg_drivers.electricals.direct_current.two_dimensions.options import ( + DC2DForwardOptions, + DC2DInversionOptions, ) from simpeg_drivers.electricals.options import ( FileControlOptions, @@ -54,28 +54,25 @@ def test_dc_p3d_fwr_run( tmp_path: Path, n_electrodes=10, n_lines=3, - refinement=(4, 6), ): # Run the forward opts = SyntheticsComponentsOptions( - method="direct current pseudo 3d", + method="direct current 2d", survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), + mesh=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=50.0, + expansion_factor=1.1, + vertical_padding=100.0, + ), model=ModelOptions(background=0.01, anomaly=10.0), ) with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: components = SyntheticsComponents(geoh5=geoh5, options=opts) - params = DCBatch2DForwardOptions.build( + params = DC2DForwardOptions.build( geoh5=geoh5, mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), topography_object=components.topography, data_object=components.survey, starting_model=components.model, @@ -83,7 +80,7 @@ def test_dc_p3d_fwr_run( line_object=components.survey.get_data("line_ids")[0] ), ) - fwr_driver = DCBatch2DForwardDriver(params) + fwr_driver = DC2DForwardDriver(params) fwr_driver.run() @@ -103,7 +100,7 @@ def test_dc_p3d_run( potential = survey.get_data("Iteration_0_potential")[0] # Run the inverse - params = DCBatch2DInversionOptions.build( + params = DC2DInversionOptions.build( geoh5=geoh5, mesh=components.mesh, drape_model=DrapeModelOptions( @@ -136,7 +133,7 @@ def test_dc_p3d_run( ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - DCBatch2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) + DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) basepath = workpath.parent with open(basepath / "lookup.json", encoding="utf8") as f: @@ -163,7 +160,6 @@ def test_dc_p3d_run( Path("./"), n_electrodes=20, n_lines=3, - refinement=(4, 4), ) test_dc_p3d_run( Path("./"), From 0b110ad2523f084aa723953307a4bb726c779f1a Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 4 Mar 2026 22:39:02 -0800 Subject: [PATCH 02/24] Continue refacotring --- .../uijson/direct_current_2d_forward.ui.json | 2 + .../direct_current_2d_inversion.ui.json | 2 + .../direct_current_batch2d_forward.ui.json | 226 ------- .../direct_current_batch2d_inversion.ui.json | 570 ----------------- .../induced_polarization_2d_forward.ui.json | 2 + .../induced_polarization_2d_inversion.ui.json | 2 + ...duced_polarization_batch2d_forward.ui.json | 237 ------- ...ced_polarization_batch2d_inversion.ui.json | 581 ------------------ simpeg_drivers/components/data.py | 18 +- simpeg_drivers/driver.py | 16 +- simpeg_drivers/electricals/base_2d.py | 37 +- simpeg_drivers/options.py | 4 +- simpeg_drivers/utils/nested.py | 63 +- simpeg_drivers/utils/utils.py | 81 ++- 14 files changed, 176 insertions(+), 1665 deletions(-) delete mode 100644 simpeg_drivers-assets/uijson/direct_current_batch2d_forward.ui.json delete mode 100644 simpeg_drivers-assets/uijson/direct_current_batch2d_inversion.ui.json delete mode 100644 simpeg_drivers-assets/uijson/induced_polarization_batch2d_forward.ui.json delete mode 100644 simpeg_drivers-assets/uijson/induced_polarization_batch2d_inversion.ui.json diff --git a/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json b/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json index 4facc425..cc08960c 100644 --- a/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json +++ b/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json @@ -33,6 +33,8 @@ "min": 1, "label": "Line number", "value": 1, + "optional": true, + "enabled": false, "tooltip": "Selects the line of data to be processed." }, "potential_channel_bool": true, diff --git a/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json b/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json index b919c9e7..3a471d24 100644 --- a/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json +++ b/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json @@ -33,6 +33,8 @@ "min": 1, "label": "Line number", "value": 1, + "optional": true, + "enabled": false, "tooltip": "Selects the line of data to be processed." }, "potential_channel": { diff --git a/simpeg_drivers-assets/uijson/direct_current_batch2d_forward.ui.json b/simpeg_drivers-assets/uijson/direct_current_batch2d_forward.ui.json deleted file mode 100644 index f481914e..00000000 --- a/simpeg_drivers-assets/uijson/direct_current_batch2d_forward.ui.json +++ /dev/null @@ -1,226 +0,0 @@ -{ - "version": "0.4.0", - "title": "Direct Current (DC) 2D Batch Forward", - "icon": "PotentialElectrode", - "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/latest/", - "conda_environment": "simpeg_drivers", - "run_command": "simpeg_drivers.driver", - "geoh5": "", - "monitoring_directory": "", - "inversion_type": "direct current pseudo 3d", - "physical_property": "conductivity", - "forward_only": true, - "data_object": { - "main": true, - "group": "Survey", - "label": "Object", - "meshType": "{275ecee9-9c24-4378-bf94-65f3c5fbe163}", - "value": "" - }, - "line_object": { - "association": "Cell", - "dataType": "Referenced", - "group": "Survey", - "main": true, - "label": "Line ID", - "parent": "data_object", - "value": "", - "tooltip": "Selects the data representing the different lines in the survey" - }, - "line_id": 1, - "potential_channel_bool": true, - "u_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Easting core cell size (m)", - "value": 25.0 - }, - "v_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Northing core cell size (m)", - "value": 25.0 - }, - "depth_core": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Depth of core (m)", - "value": 500.0 - }, - "horizontal_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Horizontal padding (m)", - "value": 1000.0 - }, - "vertical_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "dependencyType": "disabled", - "label": "Vertical padding (m)", - "value": 1000.0 - }, - "expansion_factor": { - "main": true, - "group": "Mesh and models", - "label": "Expansion factor", - "value": 1.1 - }, - "mesh": { - "group": "Mesh and models", - "main": true, - "label": "Mesh", - "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", - "optional": true, - "enabled": false, - "value": "", - "visible": true - }, - "model_type": { - "choiceList": [ - "Conductivity (S/m)", - "Resistivity (Ohm-m)" - ], - "main": true, - "group": "Mesh and models", - "label": "Model units", - "tooltip": "Select the units of the model", - "value": "Conductivity (S/m)" - }, - "starting_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Value(s)", - "property": "", - "value": 0.001 - }, - "topography_object": { - "main": true, - "group": "Topography", - "label": "Topography", - "meshType": [ - "{202c5db1-a56d-4004-9cad-baafd8899406}", - "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", - "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", - "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", - "{b020a277-90e2-4cd7-84d6-612ee3f25051}" - ], - "value": "", - "optional": true, - "enabled": true, - "tooltip": "Select a topography object to define the active cells for inversion" - }, - "topography": { - "association": [ - "Vertex", - "Cell" - ], - "dataType": "Float", - "group": "Topography", - "main": true, - "optional": true, - "enabled": false, - "label": "Elevation channel", - "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object", - "parent": "topography_object", - "dependency": "topography_object", - "dependencyType": "enabled", - "value": "", - "verbose": 2 - }, - "active_model": { - "association": "Cell", - "dataType": [ - "Referenced", - "Boolean", - "Integer" - ], - "group": "Topography", - "main": true, - "enabled": false, - "dependency": "topography_object", - "dependencyType": "disabled", - "label": "Active model", - "tooltip": "Provide the active cell Boolean model directly if topography not set", - "parent": "mesh", - "value": "" - }, - "n_cpu": { - "min": 1, - "group": "Compute", - "optional": true, - "enabled": false, - "label": "Number of CPUs", - "value": 1, - "visible": false - }, - "solver_type": { - "choiceList": [ - "Pardiso", - "Mumps" - ], - "group": "Compute", - "label": "Direct solver", - "tooltip": "Direct solver to use for the forward calculations", - "value": "Pardiso" - }, - "max_chunk_size": { - "min": 0, - "group": "Compute", - "optional": true, - "enabled": true, - "label": "Maximum chunk size (Mb)", - "value": 128, - "verbose": 3, - "visible": false, - "tooltip": "Limit the chunk size used by Dask for distributed computation" - }, - "out_group": { - "label": "SimPEG group", - "value": "", - "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", - "group": "Drag-and-drop options", - "visible": true, - "optional": true, - "enabled": false, - "tooltip": "Optionally set the SimPEG group to which results will be saved" - }, - "generate_sweep": { - "label": "Generate sweep file", - "group": "Python run preferences", - "main": true, - "value": false, - "visible": false, - "tooltip": "Generates a file for sweeping parameters instead of running the application" - }, - "files_only": { - "label": "Generate files only", - "group": "Python run preferences", - "main": true, - "value": false - }, - "cleanup": { - "main": true, - "group": "Python run preferences", - "label": "Clean directory", - "value": false - }, - "n_workers": null, - "n_threads": null, - "max_ram": "", - "performance_report": false, - "distributed_workers": "" -} diff --git a/simpeg_drivers-assets/uijson/direct_current_batch2d_inversion.ui.json b/simpeg_drivers-assets/uijson/direct_current_batch2d_inversion.ui.json deleted file mode 100644 index eea2fd30..00000000 --- a/simpeg_drivers-assets/uijson/direct_current_batch2d_inversion.ui.json +++ /dev/null @@ -1,570 +0,0 @@ -{ - "version": "0.4.0", - "title": "Direct Current (DC) 2D Batch Inversion", - "icon": "PotentialElectrode", - "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/latest/", - "conda_environment": "simpeg_drivers", - "run_command": "simpeg_drivers.driver", - "geoh5": "", - "monitoring_directory": "", - "inversion_type": "direct current pseudo 3d", - "physical_property": "conductivity", - "forward_only": false, - "data_object": { - "main": true, - "group": "Data", - "label": "Object", - "meshType": "{275ecee9-9c24-4378-bf94-65f3c5fbe163}", - "value": "" - }, - "line_object": { - "association": "Cell", - "dataType": "Referenced", - "group": "Data", - "main": true, - "label": "Line ID", - "parent": "data_object", - "value": "", - "tooltip": "Selects the data representing the different lines in the survey" - }, - "line_id": 1, - "potential_channel": { - "association": "Cell", - "dataType": "Float", - "group": "Data", - "main": true, - "label": "Potential (V/I)", - "parent": "data_object", - "value": "" - }, - "potential_uncertainty": { - "association": "Cell", - "dataType": "Float", - "group": "Data", - "main": true, - "isValue": true, - "label": "Uncertainty", - "parent": "data_object", - "property": "", - "value": 1.0 - }, - "u_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Easting core cell size (m)", - "value": 25.0 - }, - "v_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Northing core cell size (m)", - "value": 25.0 - }, - "depth_core": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Depth of core (m)", - "value": 500.0 - }, - "horizontal_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Horizontal padding (m)", - "value": 1000.0 - }, - "vertical_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "dependencyType": "disabled", - "label": "Vertical padding (m)", - "value": 1000.0 - }, - "expansion_factor": { - "main": true, - "group": "Mesh and models", - "label": "Expansion factor", - "value": 1.1 - }, - "mesh": { - "group": "Mesh and models", - "main": true, - "label": "Mesh", - "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", - "optional": true, - "enabled": false, - "value": "", - "visible": true - }, - "model_type": { - "choiceList": [ - "Conductivity (S/m)", - "Resistivity (Ohm-m)" - ], - "main": true, - "group": "Mesh and models", - "label": "Model units", - "tooltip": "Select the units of the model", - "value": "Conductivity (S/m)" - }, - "starting_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Initial Conductivity (S/m)", - "property": "", - "value": 0.001 - }, - "reference_model": { - "association": "Cell", - "dataType": "Float", - "main": true, - "group": "Mesh and models", - "isValue": true, - "optional": true, - "enabled": false, - "parent": "mesh", - "label": "Reference Conductivity (S/m)", - "property": "", - "value": 0.001 - }, - "lower_bound": { - "association": "Cell", - "main": true, - "dataType": "Float", - "group": "Mesh and models", - "isValue": true, - "parent": "mesh", - "label": "Lower bound", - "property": "", - "optional": true, - "value": 1e-08, - "enabled": false - }, - "upper_bound": { - "association": "Cell", - "main": true, - "dataType": "Float", - "group": "Mesh and models", - "isValue": true, - "parent": "mesh", - "label": "Upper bound", - "property": "", - "optional": true, - "value": 100.0, - "enabled": false - }, - "topography_object": { - "main": true, - "group": "Topography", - "label": "Topography", - "meshType": [ - "{202c5db1-a56d-4004-9cad-baafd8899406}", - "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", - "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", - "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", - "{b020a277-90e2-4cd7-84d6-612ee3f25051}" - ], - "value": "", - "optional": true, - "enabled": true, - "tooltip": "Select a topography object to define the active cells for inversion" - }, - "topography": { - "association": [ - "Vertex", - "Cell" - ], - "dataType": "Float", - "group": "Topography", - "main": true, - "optional": true, - "enabled": false, - "label": "Elevation channel", - "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object", - "parent": "topography_object", - "dependency": "topography_object", - "dependencyType": "enabled", - "value": "", - "verbose": 2 - }, - "active_model": { - "association": "Cell", - "dataType": [ - "Referenced", - "Boolean", - "Integer" - ], - "group": "Topography", - "main": true, - "enabled": false, - "dependency": "topography_object", - "dependencyType": "disabled", - "label": "Active model", - "tooltip": "Provide the active cell Boolean model directly if topography not set", - "parent": "mesh", - "value": "" - }, - "inversion_style": "voxel", - "alpha_s": { - "min": 0.0, - "group": "Regularization", - "label": "Reference weight", - "value": 1.0, - "tooltip": "Constant ratio compared to other weights. Larger values result in models that remain close to the reference model", - "dependency": "reference_model", - "dependencyType": "enabled", - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "length_scale_x": { - "min": 0.0, - "group": "Regularization", - "label": "X-smoothness weight", - "tooltip": "Larger values relative to other smoothness weights will result in x biased smoothness", - "value": 1.0, - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "length_scale_z": { - "min": 0.0, - "group": "Regularization", - "label": "Z-smoothness weight", - "tooltip": "Larger values relative to other smoothness weights will result in z biased smoothess", - "value": 1.0, - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "gradient_rotation": { - "group": "Regularization", - "association": "Cell", - "dataType": "Float", - "dataGroupType": [ - "Strike & dip", - "Dip direction & dip", - "3D vector" - ], - "label": "Gradient rotation", - "optional": true, - "enabled": false, - "parent": "mesh", - "value": "" - }, - "s_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "Smallness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 0.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": true, - "enabled": true, - "dependency": "reference_model", - "dependencyType": "enabled", - "tooltip": "Lp-norm used in the smallness term of the objective function" - }, - "x_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "X-smoothness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 2.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": false, - "enabled": true, - "tooltip": "Lp-norm used in the x-smoothness term of the objective function" - }, - "z_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "Z-smoothness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 2.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": false, - "enabled": true, - "tooltip": "Lp-norm used in the z-smoothness term of the objective function" - }, - "max_irls_iterations": { - "min": 0, - "group": "Sparse/blocky model", - "label": "Maximum IRLS iterations", - "tooltip": "Incomplete Re-weighted Least Squares iterations for non-L2 problems", - "value": 25, - "enabled": true, - "verbose": 2 - }, - "starting_chi_factor": { - "group": "Sparse/blocky model", - "label": "IRLS start chi factor", - "enabled": true, - "value": 1.0, - "tooltip": "This chi factor will be used to determine the misfit threshold after which IRLS iterations begin", - "verbose": 3 - }, - "beta_tol": { - "group": "Update IRLS directive", - "label": "Beta tolerance", - "value": 0.5, - "min": 0.0001, - "verbose": 3, - "visible": false - }, - "percentile": { - "group": "Update IRLS directive", - "label": "Percentile", - "value": 95, - "max": 100, - "min": 5, - "verbose": 3, - "visible": false - }, - "chi_factor": { - "min": 0.1, - "max": 20.0, - "precision": 1, - "lineEdit": false, - "group": "Cooling schedule/target", - "label": "Chi factor", - "value": 1.0, - "enabled": true, - "tooltip": "The global target data misfit value" - }, - "auto_scale_tiles": { - "group": "Cooling schedule/target", - "label": "Auto-scale tiles", - "value": false, - "verbose": 3, - "visible": true, - "tooltip": "Whether to auto-scale the misfit function of tiles based on chi-factor." - }, - "initial_beta_ratio": { - "min": 0.0, - "precision": 2, - "group": "Cooling schedule/target", - "optional": true, - "enabled": true, - "label": "Initial beta ratio", - "value": 100.0, - "verbose": 2, - "tooltip": "Estimate the trade-off parameter by scaling the ratio between the largest derivatives in the objective function gradients" - }, - "initial_beta": { - "min": 0.0, - "group": "Cooling schedule/target", - "optional": true, - "enabled": false, - "dependency": "initial_beta_ratio", - "dependencyType": "disabled", - "label": "Initial beta", - "value": 1.0, - "verbose": 2, - "tooltip": "Trade-off parameter between data misfit and regularization" - }, - "cooling_factor": { - "group": "Cooling schedule/target", - "label": "Beta cooling factor", - "tooltip": "Each beta cooling step will be calculated by dividing the current beta by this factor", - "value": 2.0, - "min": 1.1, - "max": 100, - "precision": 1, - "lineEdit": false, - "verbose": 2 - }, - "cooling_rate": { - "group": "Optimization", - "label": "Iterations per beta", - "value": 1, - "min": 1, - "LineEdit": false, - "max": 10, - "precision": 1, - "verbose": 2, - "enabled": true, - "tooltip": "Set the number of iterations per beta value. Use higher values for more non-linear optimization problems" - }, - "epsilon_cooling_factor": 1.2, - "max_global_iterations": { - "min": 1, - "lineEdit": false, - "group": "Optimization", - "label": "Maximum iterations", - "tooltip": "Number of L2 and IRLS iterations combined", - "value": 50, - "enabled": true - }, - "max_line_search_iterations": { - "group": "Optimization", - "label": "Maximum number of line searches", - "value": 20, - "min": 1, - "enabled": true, - "verbose": 3, - "tooltip": "Perform an Armijo backtracking line search for the provided number of iterations" - }, - "max_cg_iterations": { - "min": 0, - "group": "Optimization", - "label": "Maximum CG iterations", - "value": 30, - "enabled": true, - "verbose": 2 - }, - "tol_cg": { - "min": 0, - "group": "Optimization", - "label": "Conjugate gradient tolerance", - "value": 0.0001, - "enabled": true, - "verbose": 3 - }, - "f_min_change": { - "group": "Optimization", - "label": "Minimum change in objective function", - "value": 0.01, - "min": 1e-06, - "verbose": 3, - "enabled": true, - "tooltip": "Minimum decrease in regularization beyond which the IRLS procedure is deemed to have completed" - }, - "sens_wts_threshold": { - "group": "Update sensitivity weights directive", - "tooltip": "Update sensitivity weight threshold", - "label": "Threshold (%)", - "value": 1.0, - "max": 100.0, - "min": 0.0, - "precision": 3, - "enabled": true, - "verbose": 2 - }, - "every_iteration_bool": { - "group": "Update sensitivity weights directive", - "tooltip": "Update weights at every iteration", - "label": "Every iteration", - "value": true, - "verbose": 2, - "enabled": true - }, - "save_sensitivities": { - "group": "Update sensitivity weights directive", - "label": "Save sensitivities", - "tooltip": "Save the summed square row sensitivities to geoh5", - "value": false - }, - "n_cpu": { - "min": 1, - "group": "Compute", - "optional": true, - "enabled": false, - "label": "Number of CPUs", - "value": 1, - "visible": false - }, - "solver_type": { - "choiceList": [ - "Pardiso", - "Mumps" - ], - "group": "Compute", - "label": "Direct solver", - "tooltip": "Direct solver to use for the forward calculations", - "value": "Pardiso" - }, - "tile_spatial": 1, - "store_sensitivities": { - "choiceList": [ - "disk", - "ram" - ], - "group": "Compute", - "label": "Storage device", - "tooltip": "Use disk on a fast local SSD, and RAM elsewhere", - "value": "ram" - }, - "max_chunk_size": { - "min": 0, - "group": "Compute", - "optional": true, - "enabled": true, - "label": "Maximum chunk size (Mb)", - "value": 128, - "verbose": 3, - "visible": false, - "tooltip": "Limit the chunk size used by Dask for distributed computation" - }, - "out_group": { - "label": "SimPEG group", - "value": "", - "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", - "group": "Drag-and-drop options", - "visible": true, - "optional": true, - "enabled": false, - "tooltip": "Optionally set the SimPEG group to which results will be saved" - }, - "generate_sweep": { - "label": "Generate sweep file", - "group": "Python run preferences", - "main": true, - "value": false, - "visible": false, - "tooltip": "Generates a file for sweeping parameters instead of running the application" - }, - "files_only": { - "label": "Generate files only", - "group": "Python run preferences", - "main": true, - "value": false - }, - "cleanup": { - "main": true, - "group": "Python run preferences", - "label": "Clean directory", - "value": true - }, - "n_workers": null, - "n_threads": null, - "max_ram": "", - "performance_report": false, - "distributed_workers": "" -} diff --git a/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json index 7c3e4123..16e9f8c5 100644 --- a/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json +++ b/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json @@ -33,6 +33,8 @@ "min": 1, "label": "Line number", "value": 1, + "optional": true, + "enabled": false, "tooltip": "Selects the line of data to be processed." }, "chargeability_channel_bool": true, diff --git a/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json index 1e679ef4..60437986 100644 --- a/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json +++ b/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json @@ -33,6 +33,8 @@ "min": 1, "label": "Line number", "value": 1, + "optional": true, + "enabled": false, "tooltip": "Selects the line of data to be processed." }, "chargeability_channel": { diff --git a/simpeg_drivers-assets/uijson/induced_polarization_batch2d_forward.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_batch2d_forward.ui.json deleted file mode 100644 index fcba512f..00000000 --- a/simpeg_drivers-assets/uijson/induced_polarization_batch2d_forward.ui.json +++ /dev/null @@ -1,237 +0,0 @@ -{ - "version": "0.4.0", - "title": "Induced Polarization (IP) 2D Batch Forward", - "icon": "PotentialElectrode", - "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/latest/", - "conda_environment": "simpeg_drivers", - "run_command": "simpeg_drivers.driver", - "geoh5": "", - "monitoring_directory": "", - "inversion_type": "induced polarization pseudo 3d", - "physical_property": "chargeability", - "forward_only": true, - "data_object": { - "main": true, - "group": "Survey", - "label": "Object", - "meshType": "{275ecee9-9c24-4378-bf94-65f3c5fbe163}", - "value": "" - }, - "line_object": { - "association": "Cell", - "dataType": "Referenced", - "group": "Survey", - "main": true, - "label": "Line ID", - "parent": "data_object", - "value": "", - "tooltip": "Selects the data representing the different lines in the survey" - }, - "line_id": 1, - "chargeability_channel_bool": true, - "u_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Easting core cell size (m)", - "value": 25.0 - }, - "v_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Northing core cell size (m)", - "value": 25.0 - }, - "depth_core": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Depth of core (m)", - "value": 500.0 - }, - "horizontal_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Horizontal padding (m)", - "value": 1000.0 - }, - "vertical_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "dependencyType": "disabled", - "label": "Vertical padding (m)", - "value": 1000.0 - }, - "expansion_factor": { - "main": true, - "group": "Mesh and models", - "label": "Expansion factor", - "value": 1.1 - }, - "mesh": { - "group": "Mesh and models", - "main": true, - "optional": true, - "enabled": false, - "label": "Mesh", - "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", - "value": "", - "visible": true - }, - "model_type": { - "choiceList": [ - "Conductivity (S/m)", - "Resistivity (Ohm-m)" - ], - "main": true, - "group": "Mesh and models", - "label": "Model units", - "tooltip": "Select the units of the model", - "value": "Conductivity (S/m)" - }, - "conductivity_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Conductivity (S/m)", - "property": "", - "value": 0.001 - }, - "starting_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Value(s)", - "property": "", - "value": 0.0001 - }, - "topography_object": { - "main": true, - "group": "Topography", - "label": "Topography", - "meshType": [ - "{202c5db1-a56d-4004-9cad-baafd8899406}", - "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", - "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", - "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", - "{b020a277-90e2-4cd7-84d6-612ee3f25051}" - ], - "value": "", - "optional": true, - "enabled": true, - "tooltip": "Select a topography object to define the active cells for inversion" - }, - "topography": { - "association": [ - "Vertex", - "Cell" - ], - "dataType": "Float", - "group": "Topography", - "main": true, - "optional": true, - "enabled": false, - "label": "Elevation channel", - "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object", - "parent": "topography_object", - "dependency": "topography_object", - "dependencyType": "enabled", - "value": "", - "verbose": 2 - }, - "active_model": { - "association": "Cell", - "dataType": [ - "Referenced", - "Boolean", - "Integer" - ], - "group": "Topography", - "main": true, - "enabled": false, - "dependency": "topography_object", - "dependencyType": "disabled", - "label": "Active model", - "tooltip": "Provide the active cell Boolean model directly if topography not set", - "parent": "mesh", - "value": "" - }, - "n_cpu": { - "min": 1, - "group": "Compute", - "optional": true, - "enabled": false, - "label": "Number of CPUs", - "value": 1, - "visible": false - }, - "solver_type": { - "choiceList": [ - "Pardiso", - "Mumps" - ], - "group": "Compute", - "label": "Direct solver", - "tooltip": "Direct solver to use for the forward calculations", - "value": "Pardiso" - }, - "max_chunk_size": { - "min": 0, - "group": "Compute", - "optional": true, - "enabled": true, - "label": "Maximum chunk size (Mb)", - "value": 128, - "verbose": 3, - "visible": false, - "tooltip": "Limit the chunk size used by Dask for distributed computation" - }, - "out_group": { - "label": "SimPEG group", - "value": "", - "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", - "group": "Drag-and-drop options", - "visible": true, - "optional": true, - "enabled": false, - "tooltip": "Optionally set the SimPEG group to which results will be saved" - }, - "generate_sweep": { - "label": "Generate sweep file", - "group": "Python run preferences", - "main": true, - "value": false, - "visible": false, - "tooltip": "Generates a file for sweeping parameters instead of running the application" - }, - "files_only": { - "label": "Generate files only", - "group": "Python run preferences", - "main": true, - "value": false - }, - "cleanup": { - "main": true, - "group": "Python run preferences", - "label": "Clean directory", - "value": false - }, - "n_workers": null, - "n_threads": null, - "max_ram": "", - "performance_report": false, - "distributed_workers": "" -} diff --git a/simpeg_drivers-assets/uijson/induced_polarization_batch2d_inversion.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_batch2d_inversion.ui.json deleted file mode 100644 index cfc61c6c..00000000 --- a/simpeg_drivers-assets/uijson/induced_polarization_batch2d_inversion.ui.json +++ /dev/null @@ -1,581 +0,0 @@ -{ - "version": "0.4.0", - "title": "Induced Polarization (IP) 2D Batch Inversion", - "icon": "PotentialElectrode", - "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/latest/", - "conda_environment": "simpeg_drivers", - "run_command": "simpeg_drivers.driver", - "geoh5": "", - "monitoring_directory": "", - "inversion_type": "induced polarization pseudo 3d", - "physical_property": "chargeability", - "forward_only": false, - "data_object": { - "main": true, - "group": "Data", - "label": "Object", - "meshType": "{275ecee9-9c24-4378-bf94-65f3c5fbe163}", - "value": "" - }, - "line_object": { - "association": "Cell", - "dataType": "Referenced", - "group": "Data", - "main": true, - "label": "Line ID", - "parent": "data_object", - "value": "", - "tooltip": "Selects the data representing the different lines in the survey" - }, - "line_id": 1, - "chargeability_channel": { - "association": "Cell", - "dataType": "Float", - "group": "Data", - "main": true, - "label": "Chargeability (V/V)", - "parent": "data_object", - "value": "" - }, - "chargeability_uncertainty": { - "association": "Cell", - "dataType": "Float", - "group": "Data", - "main": true, - "isValue": true, - "label": "Uncertainty", - "parent": "data_object", - "property": "", - "value": 1.0 - }, - "u_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Easting core cell size (m)", - "value": 25.0 - }, - "v_cell_size": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Northing core cell size (m)", - "value": 25.0 - }, - "depth_core": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Depth of core (m)", - "value": 500.0 - }, - "horizontal_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "enabled": true, - "label": "Horizontal padding (m)", - "value": 1000.0 - }, - "vertical_padding": { - "min": 0.0, - "group": "Mesh and models", - "main": true, - "dependencyType": "disabled", - "label": "Vertical padding (m)", - "value": 1000.0 - }, - "expansion_factor": { - "main": true, - "group": "Mesh and models", - "label": "Expansion factor", - "value": 1.1 - }, - "mesh": { - "group": "Mesh and models", - "main": true, - "optional": true, - "enabled": false, - "label": "Mesh", - "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", - "value": "", - "visible": true - }, - "model_type": { - "choiceList": [ - "Conductivity (S/m)", - "Resistivity (Ohm-m)" - ], - "main": true, - "group": "Mesh and models", - "label": "Model units", - "tooltip": "Select the units of the model", - "value": "Conductivity (S/m)" - }, - "conductivity_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Conductivity (S/m)", - "property": "", - "value": 0.001 - }, - "starting_model": { - "association": "Cell", - "dataType": "Float", - "group": "Mesh and models", - "main": true, - "isValue": true, - "parent": "mesh", - "label": "Initial chargeability (V/V)", - "property": "", - "value": 0.0001 - }, - "reference_model": { - "association": "Cell", - "dataType": "Float", - "main": true, - "group": "Mesh and models", - "isValue": true, - "parent": "mesh", - "label": "Reference chargeability (V/V)", - "property": "", - "optional": true, - "enabled": false, - "value": 0.001 - }, - "lower_bound": { - "association": "Cell", - "main": true, - "dataType": "Float", - "group": "Mesh and models", - "isValue": true, - "parent": "mesh", - "label": "Lower bound (V/V)", - "property": "", - "optional": true, - "value": 0.0, - "enabled": true - }, - "upper_bound": { - "association": "Cell", - "main": true, - "dataType": "Float", - "group": "Mesh and models", - "isValue": true, - "parent": "mesh", - "label": "Upper bound (V/V)", - "property": "", - "optional": true, - "value": 100.0, - "enabled": false - }, - "topography_object": { - "main": true, - "group": "Topography", - "label": "Topography", - "meshType": [ - "{202c5db1-a56d-4004-9cad-baafd8899406}", - "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", - "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", - "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", - "{b020a277-90e2-4cd7-84d6-612ee3f25051}" - ], - "value": "", - "optional": true, - "enabled": true, - "tooltip": "Select a topography object to define the active cells for inversion" - }, - "topography": { - "association": [ - "Vertex", - "Cell" - ], - "dataType": "Float", - "group": "Topography", - "main": true, - "optional": true, - "enabled": false, - "label": "Elevation channel", - "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object", - "parent": "topography_object", - "dependency": "topography_object", - "dependencyType": "enabled", - "value": "", - "verbose": 2 - }, - "active_model": { - "association": "Cell", - "dataType": [ - "Referenced", - "Boolean", - "Integer" - ], - "group": "Topography", - "main": true, - "enabled": false, - "dependency": "topography_object", - "dependencyType": "disabled", - "label": "Active model", - "tooltip": "Provide the active cell Boolean model directly if topography not set", - "parent": "mesh", - "value": "" - }, - "alpha_s": { - "min": 0.0, - "group": "Regularization", - "label": "Reference weight", - "value": 1.0, - "tooltip": "Constant ratio compared to other weights. Larger values result in models that remain close to the reference model", - "dependency": "reference_model", - "dependencyType": "enabled", - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "length_scale_x": { - "min": 0.0, - "group": "Regularization", - "label": "X-smoothness weight", - "tooltip": "Larger values relative to other smoothness weights will result in x biased smoothness", - "value": 1.0, - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "length_scale_y": "", - "length_scale_z": { - "min": 0.0, - "group": "Regularization", - "label": "Z-smoothness weight", - "tooltip": "Larger values relative to other smoothness weights will result in z biased smoothess", - "value": 1.0, - "isValue": true, - "parent": "mesh", - "association": "Cell", - "dataType": "Float", - "property": "", - "enabled": true - }, - "gradient_rotation": { - "group": "Regularization", - "association": "Cell", - "dataType": "Float", - "dataGroupType": [ - "Strike & dip", - "Dip direction & dip", - "3D vector" - ], - "label": "Gradient rotation", - "optional": true, - "enabled": false, - "parent": "mesh", - "value": "" - }, - "s_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "Smallness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 0.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": true, - "enabled": true, - "dependency": "reference_model", - "dependencyType": "enabled", - "tooltip": "Lp-norm used in the smallness term of the objective function" - }, - "x_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "X-smoothness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 2.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": false, - "enabled": true, - "tooltip": "Lp-norm used in the x-smoothness term of the objective function" - }, - "y_norm": "", - "z_norm": { - "association": "Cell", - "dataType": "Float", - "group": "Sparse/blocky model", - "label": "Z-smoothness norm", - "isValue": true, - "parent": "mesh", - "property": "", - "value": 2.0, - "min": 0.0, - "max": 2.0, - "precision": 2, - "lineEdit": false, - "enabled": true, - "tooltip": "Lp-norm used in the z-smoothness term of the objective function" - }, - "max_irls_iterations": { - "min": 0, - "group": "Sparse/blocky model", - "label": "Maximum IRLS iterations", - "tooltip": "Incomplete Re-weighted Least Squares iterations for non-L2 problems", - "value": 25, - "enabled": true, - "verbose": 2 - }, - "starting_chi_factor": { - "group": "Sparse/blocky model", - "label": "IRLS start chi factor", - "enabled": true, - "value": 1.0, - "tooltip": "This chi factor will be used to determine the misfit threshold after which IRLS iterations begin", - "verbose": 3 - }, - "beta_tol": { - "group": "Update IRLS directive", - "label": "Beta tolerance", - "value": 0.5, - "min": 0.0001, - "verbose": 3, - "visible": false - }, - "percentile": { - "group": "Update IRLS directive", - "label": "Percentile", - "value": 95, - "max": 100, - "min": 5, - "verbose": 3, - "visible": false - }, - "chi_factor": { - "min": 0.1, - "max": 20.0, - "precision": 1, - "lineEdit": false, - "group": "Cooling schedule/target", - "label": "Chi factor", - "value": 1.0, - "enabled": true, - "tooltip": "The global target data misfit value" - }, - "auto_scale_tiles": { - "group": "Cooling schedule/target", - "label": "Auto-scale tiles", - "value": false, - "verbose": 3, - "visible": true, - "tooltip": "Whether to auto-scale the misfit function of tiles based on chi-factor." - }, - "initial_beta_ratio": { - "min": 0.0, - "precision": 2, - "group": "Cooling schedule/target", - "optional": true, - "enabled": true, - "label": "Initial beta ratio", - "value": 100.0, - "verbose": 2, - "tooltip": "Estimate the trade-off parameter by scaling the ratio between the largest derivatives in the objective function gradients" - }, - "initial_beta": { - "min": 0.0, - "group": "Cooling schedule/target", - "optional": true, - "enabled": false, - "dependency": "initial_beta_ratio", - "dependencyType": "disabled", - "label": "Initial beta", - "value": 1.0, - "verbose": 2, - "tooltip": "Trade-off parameter between data misfit and regularization" - }, - "cooling_factor": { - "group": "Cooling schedule/target", - "label": "Beta cooling factor", - "tooltip": "Each beta cooling step will be calculated by dividing the current beta by this factor", - "value": 2.0, - "min": 1.1, - "max": 100, - "precision": 1, - "lineEdit": false, - "verbose": 2 - }, - "cooling_rate": { - "group": "Optimization", - "label": "Iterations per beta", - "value": 1, - "min": 1, - "LineEdit": false, - "max": 10, - "precision": 1, - "verbose": 2, - "enabled": true, - "tooltip": "Set the number of iterations per beta value. Use higher values for more non-linear optimization problems" - }, - "epsilon_cooling_factor": 1.2, - "max_global_iterations": { - "min": 1, - "lineEdit": false, - "group": "Optimization", - "label": "Maximum iterations", - "tooltip": "Number of L2 and IRLS iterations combined", - "value": 50, - "enabled": true - }, - "max_line_search_iterations": { - "group": "Optimization", - "label": "Maximum number of line searches", - "value": 20, - "min": 1, - "enabled": true, - "verbose": 3, - "tooltip": "Perform an Armijo backtracking line search for the provided number of iterations" - }, - "max_cg_iterations": { - "min": 0, - "group": "Optimization", - "label": "Maximum CG iterations", - "value": 30, - "enabled": true, - "verbose": 2 - }, - "tol_cg": { - "min": 0, - "group": "Optimization", - "label": "Conjugate gradient tolerance", - "value": 0.0001, - "enabled": true, - "verbose": 3 - }, - "f_min_change": { - "group": "Optimization", - "label": "Minimum change in objective function", - "value": 0.01, - "min": 1e-06, - "verbose": 3, - "enabled": true, - "tooltip": "Minimum decrease in regularization beyond which the IRLS procedure is deemed to have completed" - }, - "sens_wts_threshold": { - "group": "Update sensitivity weights directive", - "tooltip": "Update sensitivity weight threshold", - "label": "Threshold (%)", - "value": 1.0, - "max": 100.0, - "min": 0.0, - "precision": 3, - "enabled": true, - "verbose": 2 - }, - "every_iteration_bool": { - "group": "Update sensitivity weights directive", - "tooltip": "Update weights at every iteration", - "label": "Every iteration", - "value": true, - "verbose": 2, - "enabled": true - }, - "save_sensitivities": { - "group": "Update sensitivity weights directive", - "label": "Save sensitivities", - "tooltip": "Save the summed square row sensitivities to geoh5", - "value": false - }, - "n_cpu": { - "min": 1, - "group": "Compute", - "optional": true, - "enabled": false, - "label": "Number of CPUs", - "value": 1, - "visible": false - }, - "solver_type": { - "choiceList": [ - "Pardiso", - "Mumps" - ], - "group": "Compute", - "label": "Direct solver", - "tooltip": "Direct solver to use for the forward calculations", - "value": "Pardiso" - }, - "store_sensitivities": { - "choiceList": [ - "disk", - "ram" - ], - "group": "Compute", - "label": "Storage device", - "tooltip": "Use disk on a fast local SSD, and RAM elsewhere", - "value": "ram" - }, - "max_chunk_size": { - "min": 0, - "group": "Compute", - "optional": true, - "enabled": true, - "label": "Maximum chunk size (Mb)", - "value": 128, - "verbose": 3, - "visible": false, - "tooltip": "Limit the chunk size used by Dask for distributed computation" - }, - "out_group": { - "label": "SimPEG group", - "value": "", - "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", - "group": "Drag-and-drop options", - "visible": true, - "optional": true, - "enabled": false, - "tooltip": "Optionally set the SimPEG group to which results will be saved" - }, - "generate_sweep": { - "label": "Generate sweep file", - "group": "Python run preferences", - "main": true, - "value": false, - "visible": false, - "tooltip": "Generates a file for sweeping parameters instead of running the application" - }, - "files_only": { - "label": "Generate files only", - "group": "Python run preferences", - "main": true, - "value": false - }, - "cleanup": { - "main": true, - "group": "Python run preferences", - "label": "Clean directory", - "value": true - }, - "n_workers": null, - "n_threads": null, - "max_ram": "", - "performance_report": false, - "distributed_workers": "" -} diff --git a/simpeg_drivers/components/data.py b/simpeg_drivers/components/data.py index a82076c4..6a342c5d 100644 --- a/simpeg_drivers/components/data.py +++ b/simpeg_drivers/components/data.py @@ -94,7 +94,10 @@ def _initialize(self) -> None: self.has_tensor = InversionData.check_tensor(self.params.components) self.locations = super().get_locations(self.params.data_object) - if "2d" in self.params.inversion_type: + if ( + "2d" in self.params.inversion_type + and self.params.line_selection.line_id is not None + ): self.mask = ( self.params.line_selection.line_object.values == self.params.line_selection.line_id @@ -166,7 +169,9 @@ def drape_locations(self, locations: np.ndarray) -> np.ndarray: local_tensor = drape_2_tensor(self.params.mesh) # Interpolate distance assuming always inside the mesh trace - tree = cKDTree(self.params.mesh.prisms[:, :2]) + actives = self.params.mesh.prisms[:, -1] != 1 + prisms = self.params.mesh.prisms[actives, :] + tree = cKDTree(prisms[:, :2]) rad, ind = tree.query(locations[:, :2], k=2) distance_interp = 0.0 for ii in range(2): @@ -176,7 +181,9 @@ def drape_locations(self, locations: np.ndarray) -> np.ndarray: distance_interp /= ((rad + 1e-8) ** -1.0).sum(axis=1) - return np.c_[distance_interp, locations[:, 2:]] + # Adjust elevation relative to the origin + delta = prisms[0, 2] - prisms[ind[:, 0], 2] + return np.c_[distance_interp, locations[:, 2] + delta] def get_data(self) -> tuple[list, dict, dict]: """ @@ -352,6 +359,11 @@ def create_survey(self): if "induced polarization" in self.params.inversion_type: survey.cells = self.entity.cells + if "2d" in self.params.inversion_type: + survey.line_ids = self.entity.get_data("line_ids")[0].values[ + survey_factory.sorting + ] + return survey @property diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index 05d75bed..cf333289 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -821,21 +821,13 @@ def get_tiles(self) -> dict[str, list[np.ndarray]]: n_chunks = np.max([n_chunks, 1, len(self.workers)]) tiles = [[tile] for tile in np.array_split(indices, n_chunks)] - + # Split per line for 2D inversions elif "2d" in self.params.inversion_type: tiles = [ - [ - [ - np.where( - self.params.line_selection.line_object.values == line_id - )[0] - ] - for line_id in np.unique( - self.params.line_selection.line_object.values - ) - ] + [np.where(self.params.line_selection.line_object.values == line_id)[0]] + for line_id in np.unique(self.params.line_selection.line_object.values) ] - + # Kmeans split with subsequent splitting to optimize load else: tiles = tile_locations( self.inversion_data.locations, diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 7eb6f41d..26c6ffa6 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -16,7 +16,7 @@ from geoapps_utils.utils.numerical import weighted_average from geoh5py.data import Data, IntegerData from geoh5py.groups import PropertyGroup -from geoh5py.objects import DrapeModel +from geoh5py.objects import DrapeModel, PotentialElectrode from geoh5py.shared.merging.drape_model import DrapeModelMerger from geoh5py.ui_json.ui_json import fetch_active_workspace from geoh5py.workspace import Workspace @@ -86,6 +86,9 @@ def create_mesh_by_line_id( temp_work = Workspace() for line_id in np.unique(line_ids.values): poles = get_poles_by_line_id(line_ids, line_id) + poles = np.unique(poles, axis=0) + poles = normalize_vertically(poles, line_ids.parent, drape_options.v_cell_size) + drape_model = get_drape_model( temp_work, poles, @@ -94,7 +97,8 @@ def create_mesh_by_line_id( drape_options.v_cell_size, ], drape_options.depth_core, - [0.0] * 2 + [drape_options.vertical_padding, 1], + [drape_options.horizontal_padding] * 2 + + [drape_options.vertical_padding, 1], drape_options.expansion_factor, ) drape_models.append(drape_model) @@ -104,6 +108,35 @@ def create_mesh_by_line_id( return entity +def normalize_vertically( + poles: np.ndarray, survey: PotentialElectrode, z_cell_size +) -> np.ndarray: + """ + Given a set of pole locations, normalize the vertical component to the minimum + and maximum elevations of the survey electrodes, rounded to the nearest cell thickness. + + This ensures that the drape mesh has uniform vertical discretization across all survey lines. + + :param poles: Array of pole locations to normalize. + :param survey: PotentialElectrode object containing the survey electrode locations. + :param z_cell_size: Cell size in the vertical direction for rounding the minimum and maximum + + :return: Array of pole locations with normalized vertical component. + """ + min_elev = np.min(np.r_[survey.vertices[:, 2], survey.complement.vertices[:, 2]]) + max_elev = np.max(np.r_[survey.vertices[:, 2], survey.complement.vertices[:, 2]]) + + delta = ((max_elev - min_elev) // z_cell_size + 2) * z_cell_size + min_poles_z = poles[:, 2].min() + poles[:, 2] -= min_poles_z + poles[:, 2] *= delta / poles[:, 2].max() + + # Shift back vertically and round to the nearest cell size to align + poles[:, 2] += (min_poles_z // z_cell_size - 1) * z_cell_size + + return poles + + def get_poles_by_line_id(line_ids: IntegerData, uid: int) -> np.ndarray: """Get the vertices associated with a given line ID.""" mn_mask = line_ids.values == uid diff --git a/simpeg_drivers/options.py b/simpeg_drivers/options.py index 207e685a..ec6028a7 100644 --- a/simpeg_drivers/options.py +++ b/simpeg_drivers/options.py @@ -511,7 +511,7 @@ class LineSelectionOptions(BaseModel): model_config = ConfigDict( arbitrary_types_allowed=True, ) - line_id: int = 1 + line_id: int | None = None line_object: ReferencedData @field_validator("line_object", mode="before") @@ -523,7 +523,7 @@ def validate_cell_association(cls, value): @model_validator(mode="after") def line_id_referenced(self): - if self.line_id not in self.line_object.values: + if self.line_id is not None and self.line_id not in self.line_object.values: raise ValueError("Line id isn't referenced in the line object.") return self diff --git a/simpeg_drivers/utils/nested.py b/simpeg_drivers/utils/nested.py index c47ce602..cf67a2a4 100644 --- a/simpeg_drivers/utils/nested.py +++ b/simpeg_drivers/utils/nested.py @@ -52,13 +52,49 @@ ) -def create_mesh( +def create_nested_mesh( + survey: BaseSurvey, base_mesh: TreeMesh | TensorMesh, **mesh_kwargs +) -> TreeMesh | TensorMesh: + """ + Create a nested mesh with the same extent as the input global mesh. + """ + if isinstance(base_mesh, TreeMesh): + return create_nested_treemesh(survey, base_mesh, **mesh_kwargs) + + if base_mesh.dim == 1: + return base_mesh + + return create_nested_2dmesh(survey, base_mesh) + + +def create_nested_2dmesh( + survey: BaseSurvey, + base_mesh: TensorMesh, +) -> TensorMesh: + """ + Create a nested 2D mesh using the survey line id as reference. + """ + in_cell = np.searchsorted(base_mesh.cell_centers_x, survey.locations_a[:, 0]) + unique_parts = np.unique(base_mesh.parts[in_cell]) + cells_in_part = np.hstack( + [np.where(base_mesh.parts == part)[0] for part in unique_parts] + ) + + h_x = np.diff(base_mesh.nodes_x[cells_in_part.min() : cells_in_part.max() + 2]) + h_z = base_mesh.h[1] + + return TensorMesh( + [h_x, h_z], x0=[base_mesh.nodes_x[cells_in_part.min()], base_mesh.x0[1]] + ) + + +def create_nested_treemesh( survey: BaseSurvey, - base_mesh: TreeMesh | TensorMesh, + base_mesh: TreeMesh, padding_cells: int = 8, minimum_level: int = 4, finalize: bool = True, -) -> TreeMesh | TensorMesh: +) -> TreeMesh: """ Create a nested mesh with the same extent as the input global mesh. Refinement levels are preserved only around the input locations (local survey). @@ -189,7 +225,7 @@ def _misfit_from_indices( local_survey = create_survey( simulation.survey, indices=shared_indices, channel=channel ) - local_mesh = create_mesh( + local_mesh = create_nested_mesh( local_survey, simulation.mesh, minimum_level=3, @@ -248,7 +284,7 @@ def create_simulation( args = () else: if local_mesh is None: - local_mesh = create_mesh( + local_mesh = create_nested_mesh( local_survey, simulation.mesh, minimum_level=3, @@ -270,8 +306,21 @@ def create_simulation( actives = mapping.local_active # For DCIP-2D else: - actives = simulation.active_cells - mapping = maps.IdentityMap(nP=int(actives.sum())) + actives_2d = simulation.active_cells.reshape( + simulation.mesh.shape_cells, order="F" + ) + local_active = actives_2d[simulation.mesh.parts == tile_id, :] + actives = local_active.flatten(order="F") + + # Create a projection from the global active cells to the local active cells + n_actives = simulation.active_cells.sum() + activate_ind = np.zeros(simulation.mesh.n_cells, dtype=int) + activate_ind[np.where(simulation.active_cells)[0]] = np.arange(n_actives) + activate_ind = activate_ind.reshape(simulation.mesh.shape_cells, order="F") + local_active_ind = activate_ind[ + simulation.mesh.parts == tile_id, : + ].flatten(order="F")[actives] + mapping = maps.Projection(n_actives, local_active_ind) n_actives = int(actives.sum()) if getattr(simulation, "_chiMap", None) is not None: diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 942c8578..2520d97f 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -273,6 +273,10 @@ def drape_2_tensor(drape_model: DrapeModel, return_sorting: bool = False) -> tup """ Convert a geoh5 drape model to discretize.TensorMesh. + If ghost prisms are present in the drape model, they will be skipped and the resulting + TensorMesh will have fewer cells than the DrapeModel, assuming a continuous + TensorMesh with no ghost cells. + :param: drape_model: geoh5py.DrapeModel object. :param: return_sorting: If True then return an index array that would re-sort a model in TensorMesh order to DrapeModel order. @@ -281,48 +285,74 @@ def drape_2_tensor(drape_model: DrapeModel, return_sorting: bool = False) -> tup layers = drape_model.layers # Deal with ghost points - ghosts = prisms[:, -1] == 1 - prisms = prisms[~ghosts, :] + actives = prisms[:, -1] != 1 - nu_layers = np.unique(prisms[:, -1]) + nu_layers = np.unique(prisms[actives, -1]) if len(nu_layers) > 1: raise ValueError( "Drape model conversion to TensorMesh must have uniform number of layers." ) n_layers = nu_layers[0].astype(int) - filt_layers = ghosts[layers[:, 0].astype(int)] - layers = layers[~filt_layers, :] + n_columns = actives.sum() + + # Sorting array from DrapeModel to TensorMesh order row-wise, skipping ghost points + sorting = np.arange(n_columns * n_layers) + sorting = sorting.reshape(n_layers, n_columns, order="C") + sorting = np.argsort(sorting[::-1].T.flatten()) + filt_layers = actives[layers[:, 0].astype(int)] + layers = layers[filt_layers, :] hz = np.r_[ prisms[0, 2] - layers[0, 2], -np.diff(layers[:n_layers, 2]), ][::-1] - x = compute_alongline_distance(prisms[:, :2], ordered=False) - dx = np.diff(x) - cell_width = np.r_[dx[0], (dx[:-1] + dx[1:]) / 2.0, dx[-1]] - h = [cell_width, hz] - origin = [0, layers[:, 2].min()] + # Skip indices for ghost points + count = -1 + part = 0 + parts = [] + cell_widths = [] + section = [] + for ii, active in enumerate(actives): + if not active: + sorting[sorting > count] += 1 + count += 1 + + if section: + cell_widths.append(cell_width_from_centers(np.vstack(section))) + parts.append(np.full(len(section), part)) + section = [] + part += 1 + else: + section.append(np.c_[prisms[ii, 0], 0]) + count += n_layers + + cell_widths.append(cell_width_from_centers(np.vstack(section))) + parts.append(np.full(len(section), part)) + + h = [np.hstack(cell_widths), hz] + origin = [0, prisms[0, 2] - hz.sum()] mesh = TensorMesh(h, origin=origin) + mesh.parts = np.hstack(parts) # Assign part numbers to cells if return_sorting: - sorting = np.arange(mesh.n_cells) - sorting = sorting.reshape(mesh.shape_cells[1], mesh.shape_cells[0], order="C") - sorting = np.argsort(sorting[::-1].T.flatten()) - - # Skip indices for ghost points - count = -1 - for ghost in ghosts: - if ghost: - sorting[sorting > count] += 1 - count += 1 - else: - count += n_layers - return (mesh, sorting) - else: - return mesh + + return mesh + + +def cell_width_from_centers(centers: np.ndarray) -> np.ndarray: + """ + Compute cell widths from cell center locations. + + :param centers: n x 3 array of cell center locations + + :returns: n-1 array of cell widths + """ + x = compute_alongline_distance(centers[:, :2]) + half_dx = np.diff(x) / 2.0 + return np.r_[half_dx[0] * 2, (half_dx[:-1] + half_dx[1:]), half_dx[-1] * 2] def floating_active(mesh: TensorMesh | TreeMesh, active: np.ndarray): @@ -381,6 +411,7 @@ def get_drape_model( np.c_[locations[order[-1], :]].T, ] ) + distances = compute_alongline_distance(xy_smooth) x_interp = interp1d(distances[:, 0], xy_smooth[:, 0], fill_value="extrapolate") y_interp = interp1d(distances[:, 0], xy_smooth[:, 1], fill_value="extrapolate") From f21b36302982339a065a2158c25e060a9057844c Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 4 Mar 2026 22:50:53 -0800 Subject: [PATCH 03/24] Use max relief of all lines --- simpeg_drivers/electricals/base_2d.py | 214 +++----------------------- 1 file changed, 25 insertions(+), 189 deletions(-) diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 26c6ffa6..3f774ad5 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -84,10 +84,13 @@ def create_mesh_by_line_id( """ drape_models = [] temp_work = Workspace() + + relief = get_max_line_relief(line_ids, drape_options.v_cell_size) + for line_id in np.unique(line_ids.values): poles = get_poles_by_line_id(line_ids, line_id) poles = np.unique(poles, axis=0) - poles = normalize_vertically(poles, line_ids.parent, drape_options.v_cell_size) + poles = normalize_vertically(poles, relief) drape_model = get_drape_model( temp_work, @@ -108,31 +111,38 @@ def create_mesh_by_line_id( return entity -def normalize_vertically( - poles: np.ndarray, survey: PotentialElectrode, z_cell_size -) -> np.ndarray: +def get_max_line_relief(line_ids, z_cell_size) -> float: + """ + Get the maximum relief across all survey lines, rounded to the nearest cell thickness. + + :param line_ids: IntegerData object containing the line IDs for each vertex. + :param z_cell_size: Cell size in the vertical direction for the drape mesh. + """ + max_relief = 0 + for line_id in np.unique(line_ids.values): + poles = get_poles_by_line_id(line_ids, line_id) + max_relief = np.maximum(poles[:, 2].max() - poles[:, 2].min(), max_relief) + + return (max_relief // z_cell_size + 2) * z_cell_size + + +def normalize_vertically(poles: np.ndarray, relief: float) -> np.ndarray: """ - Given a set of pole locations, normalize the vertical component to the minimum - and maximum elevations of the survey electrodes, rounded to the nearest cell thickness. + Given a set of pole locations, normalize the vertical component to the maximum relief across all lines. This ensures that the drape mesh has uniform vertical discretization across all survey lines. :param poles: Array of pole locations to normalize. - :param survey: PotentialElectrode object containing the survey electrode locations. - :param z_cell_size: Cell size in the vertical direction for rounding the minimum and maximum + :param relief: Maximum relief across all survey lines, rounded to the nearest cell thickness. :return: Array of pole locations with normalized vertical component. """ - min_elev = np.min(np.r_[survey.vertices[:, 2], survey.complement.vertices[:, 2]]) - max_elev = np.max(np.r_[survey.vertices[:, 2], survey.complement.vertices[:, 2]]) - - delta = ((max_elev - min_elev) // z_cell_size + 2) * z_cell_size min_poles_z = poles[:, 2].min() poles[:, 2] -= min_poles_z - poles[:, 2] *= delta / poles[:, 2].max() + poles[:, 2] *= relief / poles[:, 2].max() - # Shift back vertically and round to the nearest cell size to align - poles[:, 2] += (min_poles_z // z_cell_size - 1) * z_cell_size + # Shift back vertically + poles[:, 2] += min_poles_z return poles @@ -153,177 +163,3 @@ def get_poles_by_line_id(line_ids: IntegerData, uid: int) -> np.ndarray: ], ] ) - - -# class BaseBatch2DDriver(LineSweepDriver): -# """Base class for batch 2D DC and IP forward and inversion drivers.""" -# -# _params_class: type[BaseForwardOptions | BaseInversionOptions] -# _params_2d_class: type[BaseForwardOptions | BaseInversionOptions] -# -# _model_list: list[str] = [] -# -# def __init__(self, params): -# super().__init__(params) -# if params.file_control.files_only: -# sys.exit("Files written") -# -# def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID | float]: -# """ -# Transfer models from the input parameters to the output drape mesh. -# -# :param mesh: Destination DrapeModel object. -# """ -# models = {"starting_model": self.batch2d_params.models.starting_model} -# -# for model in self._model_list: -# models[model] = getattr(self.batch2d_params, model, None) -# -# if not self.batch2d_params.forward_only: -# for model in ["reference_model", "lower_bound", "upper_bound"]: -# models[model] = getattr(self.batch2d_params.models, model) -# -# if self.batch2d_params.models.gradient_rotation is not None: -# group_properties = {} -# for prop in self.batch2d_params.models.gradient_rotation.properties: -# model = self.batch2d_params.mesh.get_data(prop)[0] -# group_properties[model.name] = model -# -# models.update(group_properties) -# -# if self.batch2d_params.mesh is not None: -# xyz_in = get_locations(self.workspace, self.batch2d_params.mesh) -# xyz_out = mesh.centroids -# -# for name, model in models.items(): -# if model is None: -# continue -# elif isinstance(model, Data): -# model_values = weighted_average( -# xyz_in, xyz_out, [model.values], n=1 -# )[0] -# else: -# model_values = model * np.ones(len(xyz_out)) -# -# model_object = mesh.add_data({name: {"values": model_values}}) -# models[name] = model_object -# -# if ( -# not self.batch2d_params.forward_only -# and self.batch2d_params.models.gradient_rotation is not None -# ): -# pg = PropertyGroup( -# mesh, -# properties=[models[prop] for prop in group_properties], -# property_group_type=self.batch2d_params.models.gradient_rotation.property_group_type, -# ) -# models["gradient_rotation"] = pg -# del models["azimuth"] -# del models["dip"] -# -# return models -# -# def write_files(self, lookup): -# """Write ui.geoh5 and ui.json files for sweep trials.""" -# -# kwargs_2d = {} -# with fetch_active_workspace(self.workspace, mode="r+"): -# for uid, trial in lookup.items(): -# if trial["status"] != "pending": -# continue -# -# filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" -# -# if filepath.exists(): -# warnings.warn( -# f"File {filepath} already exists but 'status' marked as 'pending'. " -# "Over-writing file." -# ) -# filepath.unlink() -# -# with Workspace.create(filepath) as iter_workspace: -# cell_mask: np.ndarray = ( -# self.batch2d_params.line_selection.line_object.values -# == trial["line_id"] -# ) -# -# if not np.any(cell_mask): -# continue -# -# receiver_entity = extract_dcip_survey( -# iter_workspace, self.inversion_data.entity, cell_mask -# ) -# current_entity = receiver_entity.current_electrodes -# receiver_locs = np.vstack( -# [receiver_entity.vertices, current_entity.vertices] -# ) -# -# mesh = get_drape_model( -# iter_workspace, -# "Models", -# receiver_locs, -# [ -# self.batch2d_params.drape_model.u_cell_size, -# self.batch2d_params.drape_model.v_cell_size, -# ], -# self.batch2d_params.drape_model.depth_core, -# [self.batch2d_params.drape_model.horizontal_padding] * 2 -# + [self.batch2d_params.drape_model.vertical_padding, 1], -# self.batch2d_params.drape_model.expansion_factor, -# )[0] -# -# model_parameters = self.transfer_models(mesh) -# -# for key in self._params_2d_class.model_fields: -# param = getattr(self.batch2d_params, key, None) -# if key not in ["title", "inversion_type"]: -# kwargs_2d[key] = param -# -# self.batch2d_params.active_cells.topography_object.copy( -# parent=iter_workspace, copy_children=True -# ) -# -# kwargs_2d.update( -# dict( -# **{ -# "geoh5": iter_workspace, -# "mesh": mesh, -# "data_object": receiver_entity, -# "line_selection": LineSelectionOptions( -# line_object=receiver_entity.get_data( -# self.batch2d_params.line_selection.line_object.name -# )[0], -# line_id=trial["line_id"], -# ), -# "out_group": None, -# }, -# **model_parameters, -# ) -# ) -# -# params = self._params_2d_class(**kwargs_2d) -# params.write_ui_json(Path(self.working_directory) / f"{uid}.ui.json") -# -# lookup[uid]["status"] = "written" -# -# _ = self.update_lookup(lookup) # pylint: disable=no-member -# -# @property -# def inversion_data(self) -> InversionData: -# """Inversion data""" -# if getattr(self, "_inversion_data", None) is None: -# with fetch_active_workspace(self.workspace, mode="r+"): -# self._inversion_data = InversionData( -# self.workspace, self.batch2d_params -# ) -# -# return self._inversion_data -# -# @property -# def inversion_topography(self): -# """Inversion topography""" -# if getattr(self, "_inversion_topography", None) is None: -# self._inversion_topography = InversionTopography( -# self.workspace, self.batch2d_params -# ) -# return self._inversion_topography From 777166334231fa94e5ef86d053a10d05d8edfcd5 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 4 Mar 2026 23:06:14 -0800 Subject: [PATCH 04/24] Full run of test --- simpeg_drivers/electricals/base_2d.py | 4 ++-- simpeg_drivers/electricals/options.py | 12 ---------- tests/run_tests/driver_dc_b2d_test.py | 33 ++++----------------------- 3 files changed, 6 insertions(+), 43 deletions(-) diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 3f774ad5..13569700 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -50,8 +50,8 @@ def inversion_mesh(self) -> InversionMesh: """Inversion mesh""" if getattr(self, "_inversion_mesh", None) is None: with fetch_active_workspace(self.workspace, mode="r+"): - entity = self.params.mesh - if entity is None: + entity = None + if self.params.mesh is None: entity = create_mesh_by_line_id( self.workspace, self.params.line_selection.line_object, diff --git a/simpeg_drivers/electricals/options.py b/simpeg_drivers/electricals/options.py index 7a528abf..7cca3219 100644 --- a/simpeg_drivers/electricals/options.py +++ b/simpeg_drivers/electricals/options.py @@ -23,15 +23,3 @@ class IPModelOptions(ConductivityModelOptions): """ lower_bound: float | FloatData | None = 0 - - -class FileControlOptions(BaseModel): - """ - File control parameters for pseudo 3D simulations. - - :param files_only: Boolean to only write files. - :param cleanup: Boolean to cleanup files. - """ - - files_only: bool = False - cleanup: bool = True diff --git a/tests/run_tests/driver_dc_b2d_test.py b/tests/run_tests/driver_dc_b2d_test.py index 3a3822d2..940bab6e 100644 --- a/tests/run_tests/driver_dc_b2d_test.py +++ b/tests/run_tests/driver_dc_b2d_test.py @@ -11,10 +11,8 @@ from __future__ import annotations -import json from pathlib import Path -from geoh5py.groups import SimPEGGroup from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( @@ -25,9 +23,6 @@ DC2DForwardOptions, DC2DInversionOptions, ) -from simpeg_drivers.electricals.options import ( - FileControlOptions, -) from simpeg_drivers.options import ( DrapeModelOptions, LineSelectionOptions, @@ -36,7 +31,6 @@ SyntheticsComponents, ) from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, ModelOptions, SurveyOptions, SyntheticsComponentsOptions, @@ -47,7 +41,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 1.0878862593748388, "phi_d": 2050, "phi_m": 37.8} +target_run = {"data_norm": 1.101767837151429, "phi_d": 2210, "phi_m": 21.4} def test_dc_p3d_fwr_run( @@ -95,7 +89,7 @@ def test_dc_p3d_run( with Workspace(workpath) as geoh5: components = SyntheticsComponents(geoh5) - fwr_group = geoh5.get_entity("Direct Current (DC) 2D Batch Forward")[0] + fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] survey = fwr_group.get_entity("survey")[0] potential = survey.get_data("Iteration_0_potential")[0] @@ -103,14 +97,6 @@ def test_dc_p3d_run( params = DC2DInversionOptions.build( geoh5=geoh5, mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), topography_object=components.topography, data_object=potential.parent, potential_channel=potential, @@ -129,24 +115,13 @@ def test_dc_p3d_run( percentile=100, upper_bound=10, cooling_rate=1, - file_control=FileControlOptions(cleanup=False), ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - - basepath = workpath.parent - with open(basepath / "lookup.json", encoding="utf8") as f: - lookup = json.load(f) - middle_line_id = next(k for k, v in lookup.items() if v["line_id"] == 101) - - with Workspace(basepath / f"{middle_line_id}.ui.geoh5", mode="r") as workspace: - middle_inversion_group = next( - k for k in workspace.groups if isinstance(k, SimPEGGroup) - ) + driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) output = get_inversion_output( - basepath / f"{middle_line_id}.ui.geoh5", middle_inversion_group.uid + driver.params.geoh5.h5file, driver.params.out_group.uid ) if geoh5.open(): output["data"] = potential.values From 719eb275b2fa416eda8d1a762cc98041217133dc Mon Sep 17 00:00:00 2001 From: dominiquef Date: Sat, 7 Mar 2026 17:42:10 -0800 Subject: [PATCH 05/24] Continue update --- .../uijson/direct_current_2d_forward.ui.json | 9 +- .../direct_current_2d_inversion.ui.json | 13 +- .../induced_polarization_2d_forward.ui.json | 9 +- .../induced_polarization_2d_inversion.ui.json | 13 +- simpeg_drivers/components/data.py | 34 +- simpeg_drivers/electricals/base_2d.py | 15 +- .../direct_current/two_dimensions/options.py | 4 +- simpeg_drivers/electricals/driver.py | 252 --------------- .../two_dimensions/options.py | 4 +- simpeg_drivers/line_sweep/__init__.py | 9 - simpeg_drivers/line_sweep/driver.py | 303 ------------------ simpeg_drivers/options.py | 4 +- simpeg_drivers/utils/nested.py | 48 +-- simpeg_drivers/utils/surveys.py | 26 +- simpeg_drivers/utils/utils.py | 2 +- .../driver_2d_rotated_gradients_test.py | 201 ------------ ...=> driver_dc_2d_rotated_gradients_test.py} | 66 ++-- tests/run_tests/driver_dc_2d_test.py | 129 +++++--- tests/run_tests/driver_dc_b2d_test.py | 143 --------- tests/run_tests/driver_ip_b2d_test.py | 180 ----------- 20 files changed, 220 insertions(+), 1244 deletions(-) delete mode 100644 simpeg_drivers/electricals/driver.py delete mode 100644 simpeg_drivers/line_sweep/__init__.py delete mode 100644 simpeg_drivers/line_sweep/driver.py delete mode 100644 tests/run_tests/driver_2d_rotated_gradients_test.py rename tests/run_tests/{driver_dc_b2d_rotated_gradients_test.py => driver_dc_2d_rotated_gradients_test.py} (73%) delete mode 100644 tests/run_tests/driver_dc_b2d_test.py delete mode 100644 tests/run_tests/driver_ip_b2d_test.py diff --git a/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json b/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json index cc08960c..708e2548 100644 --- a/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json +++ b/simpeg_drivers-assets/uijson/direct_current_2d_forward.ui.json @@ -19,12 +19,17 @@ }, "line_object": { "association": "Cell", - "dataType": "Referenced", + "dataType": [ + "Integer", + "Referenced" + ], "group": "Survey", "main": true, "label": "Line ID", "parent": "data_object", "value": "", + "optional": true, + "enabled": false, "tooltip": "Selects the data representing the different lines in the survey." }, "line_id": { @@ -35,6 +40,8 @@ "value": 1, "optional": true, "enabled": false, + "dependency": "line_object", + "dependencyType": "enabled", "tooltip": "Selects the line of data to be processed." }, "potential_channel_bool": true, diff --git a/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json b/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json index 3a471d24..63263429 100644 --- a/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json +++ b/simpeg_drivers-assets/uijson/direct_current_2d_inversion.ui.json @@ -19,22 +19,29 @@ }, "line_object": { "association": "Cell", - "dataType": "Referenced", - "group": "Data", + "dataType": [ + "Integer", + "Referenced" + ], + "group": "Survey", "main": true, "label": "Line ID", "parent": "data_object", "value": "", + "optional": true, + "enabled": false, "tooltip": "Selects the data representing the different lines in the survey." }, "line_id": { - "group": "Data", + "group": "Survey", "main": true, "min": 1, "label": "Line number", "value": 1, "optional": true, "enabled": false, + "dependency": "line_object", + "dependencyType": "enabled", "tooltip": "Selects the line of data to be processed." }, "potential_channel": { diff --git a/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json index 16e9f8c5..27854c2c 100644 --- a/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json +++ b/simpeg_drivers-assets/uijson/induced_polarization_2d_forward.ui.json @@ -19,12 +19,17 @@ }, "line_object": { "association": "Cell", - "dataType": "Referenced", + "dataType": [ + "Integer", + "Referenced" + ], "group": "Survey", "main": true, "label": "Line ID", "parent": "data_object", "value": "", + "optional": true, + "enabled": false, "tooltip": "Selects the data representing the different lines in the survey." }, "line_id": { @@ -35,6 +40,8 @@ "value": 1, "optional": true, "enabled": false, + "dependency": "line_object", + "dependencyType": "enabled", "tooltip": "Selects the line of data to be processed." }, "chargeability_channel_bool": true, diff --git a/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json b/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json index 60437986..ba6947d2 100644 --- a/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json +++ b/simpeg_drivers-assets/uijson/induced_polarization_2d_inversion.ui.json @@ -19,22 +19,29 @@ }, "line_object": { "association": "Cell", - "dataType": "Referenced", - "group": "Data", + "dataType": [ + "Integer", + "Referenced" + ], + "group": "Survey", "main": true, "label": "Line ID", "parent": "data_object", "value": "", + "optional": true, + "enabled": false, "tooltip": "Selects the data representing the different lines in the survey." }, "line_id": { - "group": "Data", + "group": "Survey", "main": true, "min": 1, "label": "Line number", "value": 1, "optional": true, "enabled": false, + "dependency": "line_object", + "dependencyType": "enabled", "tooltip": "Selects the line of data to be processed." }, "chargeability_channel": { diff --git a/simpeg_drivers/components/data.py b/simpeg_drivers/components/data.py index 6a342c5d..75f5e3c9 100644 --- a/simpeg_drivers/components/data.py +++ b/simpeg_drivers/components/data.py @@ -17,10 +17,10 @@ import numpy as np from geoh5py.objects import LargeLoopGroundTEMReceivers, PotentialElectrode -from scipy.sparse import csgraph, csr_matrix from scipy.spatial import cKDTree from simpeg.electromagnetics.static.utils.static_utils import geometric_factor +from simpeg_drivers.utils.surveys import get_parts_from_electrodes from simpeg_drivers.utils.utils import drape_2_tensor from .factories import ( @@ -140,19 +140,7 @@ def parts(self): Return parts indices from the entity. """ if isinstance(self.entity, PotentialElectrode): - edge_array = csr_matrix( - ( - np.ones(self.entity.n_cells * 2), - ( - np.kron(self.entity.cells[:, 0], [1, 1]), - self.entity.cells.flatten(), - ), - ), - shape=(self.entity.n_vertices, self.entity.n_vertices), - ) - - connections = csgraph.connected_components(edge_array)[1] - return connections[self.entity.cells[:, 0]] + return get_parts_from_electrodes(self.entity) if isinstance(self.entity, LargeLoopGroundTEMReceivers): return self.entity.tx_id_property.values @@ -360,7 +348,7 @@ def create_survey(self): survey.cells = self.entity.cells if "2d" in self.params.inversion_type: - survey.line_ids = self.entity.get_data("line_ids")[0].values[ + survey.line_ids = self.params.line_selection.line_object.values[ survey_factory.sorting ] @@ -398,11 +386,17 @@ def update_params(self, data_dict, uncert_dict): setattr(self.params, f"{comp}_uncertainty", uncert_dict[comp]) if getattr(self.params, "line_selection", None) is not None: - new_line = self.params.line_selection.line_object.copy( - parent=self.entity, - values=self.params.line_selection.line_object.values[self.mask], - ) - self.params.line_selection.line_object = new_line + if self.params.line_selection.line_object is None: + parts = get_parts_from_electrodes(self.entity) + _, u_part = np.unique(parts, return_inverse=True) + line_ids = self.entity.add_data({"Line IDs": {"values": u_part + 1}}) + else: + line_ids = self.params.line_selection.line_object.copy( + parent=self.entity, + values=self.params.line_selection.line_object.values[self.mask], + ) + + self.params.line_selection.line_object = line_ids @property def survey(self): diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 13569700..8906e725 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -12,28 +12,17 @@ from __future__ import annotations import numpy as np -from geoapps_utils.utils.locations import get_locations -from geoapps_utils.utils.numerical import weighted_average -from geoh5py.data import Data, IntegerData -from geoh5py.groups import PropertyGroup -from geoh5py.objects import DrapeModel, PotentialElectrode +from geoh5py.data import IntegerData +from geoh5py.objects import DrapeModel from geoh5py.shared.merging.drape_model import DrapeModelMerger from geoh5py.ui_json.ui_json import fetch_active_workspace from geoh5py.workspace import Workspace -from simpeg_drivers.components.data import InversionData from simpeg_drivers.components.meshes import InversionMesh -from simpeg_drivers.components.topography import InversionTopography -from simpeg_drivers.components.windows import InversionWindow from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.line_sweep.driver import LineSweepDriver from simpeg_drivers.options import ( - BaseForwardOptions, - BaseInversionOptions, DrapeModelOptions, - LineSelectionOptions, ) -from simpeg_drivers.utils.surveys import extract_dcip_survey from simpeg_drivers.utils.utils import get_drape_model diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py index 89292cb8..468efad5 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py @@ -47,7 +47,7 @@ class DC2DForwardOptions(BaseForwardOptions): data_object: PotentialElectrode potential_channel_bool: bool = True - line_selection: LineSelectionOptions + line_selection: LineSelectionOptions = LineSelectionOptions() mesh: DrapeModel | None = None drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions @@ -75,7 +75,7 @@ class DC2DInversionOptions(BaseInversionOptions): data_object: PotentialElectrode potential_channel: FloatData potential_uncertainty: float | FloatData | None = None - line_selection: LineSelectionOptions + line_selection: LineSelectionOptions = LineSelectionOptions() mesh: DrapeModel | None = None drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions diff --git a/simpeg_drivers/electricals/driver.py b/simpeg_drivers/electricals/driver.py deleted file mode 100644 index 0571d901..00000000 --- a/simpeg_drivers/electricals/driver.py +++ /dev/null @@ -1,252 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -import sys -import uuid -import warnings -from pathlib import Path - -import numpy as np -from geoapps_utils.utils.locations import get_locations -from geoapps_utils.utils.numerical import weighted_average -from geoh5py.data import Data -from geoh5py.groups import PropertyGroup -from geoh5py.objects import DrapeModel -from geoh5py.ui_json.ui_json import fetch_active_workspace -from geoh5py.workspace import Workspace - -from simpeg_drivers.components.data import InversionData -from simpeg_drivers.components.meshes import InversionMesh -from simpeg_drivers.components.topography import InversionTopography -from simpeg_drivers.components.windows import InversionWindow -from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.line_sweep.driver import LineSweepDriver -from simpeg_drivers.options import ( - BaseForwardOptions, - BaseInversionOptions, - LineSelectionOptions, -) -from simpeg_drivers.utils.surveys import extract_dcip_survey -from simpeg_drivers.utils.utils import get_drape_model - - -class Base2DDriver(InversionDriver): - """Base class for 2D DC and IP forward and inversion drivers.""" - - @property - def inversion_mesh(self) -> InversionMesh: - """Inversion mesh""" - if getattr(self, "_inversion_mesh", None) is None: - with fetch_active_workspace(self.workspace, mode="r+"): - if self.params.mesh is None: - self.params.mesh = self.create_drape_mesh() - - self._inversion_mesh = InversionMesh(self.workspace, self.params) - return self._inversion_mesh - - def create_drape_mesh(self) -> DrapeModel: - """Create a drape mesh for the inversion.""" - current_entity = self.params.data_object.current_electrodes - receiver_locs = np.vstack( - [self.params.data_object.vertices, current_entity.vertices] - ) - with fetch_active_workspace(self.workspace): - mesh = get_drape_model( - self.workspace, - "Models", - receiver_locs, - [ - self.params.drape_model.u_cell_size, - self.params.drape_model.v_cell_size, - ], - self.params.drape_model.depth_core, - [self.params.drape_model.horizontal_padding] * 2 - + [self.params.drape_model.vertical_padding, 1], - self.params.drape_model.expansion_factor, - )[0] - - return mesh - - -class BaseBatch2DDriver(LineSweepDriver): - """Base class for batch 2D DC and IP forward and inversion drivers.""" - - _params_class: type[BaseForwardOptions | BaseInversionOptions] - _params_2d_class: type[BaseForwardOptions | BaseInversionOptions] - - _model_list: list[str] = [] - - def __init__(self, params): - super().__init__(params) - if params.file_control.files_only: - sys.exit("Files written") - - def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID | float]: - """ - Transfer models from the input parameters to the output drape mesh. - - :param mesh: Destination DrapeModel object. - """ - models = {"starting_model": self.batch2d_params.models.starting_model} - - for model in self._model_list: - models[model] = getattr(self.batch2d_params, model, None) - - if not self.batch2d_params.forward_only: - for model in ["reference_model", "lower_bound", "upper_bound"]: - models[model] = getattr(self.batch2d_params.models, model) - - if self.batch2d_params.models.gradient_rotation is not None: - group_properties = {} - for prop in self.batch2d_params.models.gradient_rotation.properties: - model = self.batch2d_params.mesh.get_data(prop)[0] - group_properties[model.name] = model - - models.update(group_properties) - - if self.batch2d_params.mesh is not None: - xyz_in = get_locations(self.workspace, self.batch2d_params.mesh) - xyz_out = mesh.centroids - - for name, model in models.items(): - if model is None: - continue - elif isinstance(model, Data): - model_values = weighted_average( - xyz_in, xyz_out, [model.values], n=1 - )[0] - else: - model_values = model * np.ones(len(xyz_out)) - - model_object = mesh.add_data({name: {"values": model_values}}) - models[name] = model_object - - if ( - not self.batch2d_params.forward_only - and self.batch2d_params.models.gradient_rotation is not None - ): - pg = PropertyGroup( - mesh, - properties=[models[prop] for prop in group_properties], - property_group_type=self.batch2d_params.models.gradient_rotation.property_group_type, - ) - models["gradient_rotation"] = pg - del models["azimuth"] - del models["dip"] - - return models - - def write_files(self, lookup): - """Write ui.geoh5 and ui.json files for sweep trials.""" - - kwargs_2d = {} - with fetch_active_workspace(self.workspace, mode="r+"): - for uid, trial in lookup.items(): - if trial["status"] != "pending": - continue - - filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" - - if filepath.exists(): - warnings.warn( - f"File {filepath} already exists but 'status' marked as 'pending'. " - "Over-writing file." - ) - filepath.unlink() - - with Workspace.create(filepath) as iter_workspace: - cell_mask: np.ndarray = ( - self.batch2d_params.line_selection.line_object.values - == trial["line_id"] - ) - - if not np.any(cell_mask): - continue - - receiver_entity = extract_dcip_survey( - iter_workspace, self.inversion_data.entity, cell_mask - ) - current_entity = receiver_entity.current_electrodes - receiver_locs = np.vstack( - [receiver_entity.vertices, current_entity.vertices] - ) - - mesh = get_drape_model( - iter_workspace, - "Models", - receiver_locs, - [ - self.batch2d_params.drape_model.u_cell_size, - self.batch2d_params.drape_model.v_cell_size, - ], - self.batch2d_params.drape_model.depth_core, - [self.batch2d_params.drape_model.horizontal_padding] * 2 - + [self.batch2d_params.drape_model.vertical_padding, 1], - self.batch2d_params.drape_model.expansion_factor, - )[0] - - model_parameters = self.transfer_models(mesh) - - for key in self._params_2d_class.model_fields: - param = getattr(self.batch2d_params, key, None) - if key not in ["title", "inversion_type"]: - kwargs_2d[key] = param - - self.batch2d_params.active_cells.topography_object.copy( - parent=iter_workspace, copy_children=True - ) - - kwargs_2d.update( - dict( - **{ - "geoh5": iter_workspace, - "mesh": mesh, - "data_object": receiver_entity, - "line_selection": LineSelectionOptions( - line_object=receiver_entity.get_data( - self.batch2d_params.line_selection.line_object.name - )[0], - line_id=trial["line_id"], - ), - "out_group": None, - }, - **model_parameters, - ) - ) - - params = self._params_2d_class(**kwargs_2d) - params.write_ui_json(Path(self.working_directory) / f"{uid}.ui.json") - - lookup[uid]["status"] = "written" - - _ = self.update_lookup(lookup) # pylint: disable=no-member - - @property - def inversion_data(self) -> InversionData: - """Inversion data""" - if getattr(self, "_inversion_data", None) is None: - with fetch_active_workspace(self.workspace, mode="r+"): - self._inversion_data = InversionData( - self.workspace, self.batch2d_params - ) - - return self._inversion_data - - @property - def inversion_topography(self): - """Inversion topography""" - if getattr(self, "_inversion_topography", None) is None: - self._inversion_topography = InversionTopography( - self.workspace, self.batch2d_params - ) - return self._inversion_topography diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py index 4322efdc..29a33611 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py @@ -48,7 +48,7 @@ class IP2DForwardOptions(BaseForwardOptions): data_object: PotentialElectrode chargeability_channel_bool: bool = True - line_selection: LineSelectionOptions + line_selection: LineSelectionOptions = LineSelectionOptions() mesh: DrapeModel | None = None drape_model: DrapeModelOptions = DrapeModelOptions() models: IPModelOptions @@ -76,7 +76,7 @@ class IP2DInversionOptions(BaseInversionOptions): data_object: PotentialElectrode chargeability_channel: FloatData chargeability_uncertainty: float | FloatData | None = None - line_selection: LineSelectionOptions + line_selection: LineSelectionOptions = LineSelectionOptions() mesh: DrapeModel | None = None drape_model: DrapeModelOptions = DrapeModelOptions() models: IPModelOptions diff --git a/simpeg_drivers/line_sweep/__init__.py b/simpeg_drivers/line_sweep/__init__.py deleted file mode 100644 index df32b204..00000000 --- a/simpeg_drivers/line_sweep/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/line_sweep/driver.py b/simpeg_drivers/line_sweep/driver.py deleted file mode 100644 index c5351ce9..00000000 --- a/simpeg_drivers/line_sweep/driver.py +++ /dev/null @@ -1,303 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -import json -import re -from pathlib import Path - -import numpy as np -from geoapps_utils.param_sweeps.driver import SweepDriver, SweepParams -from geoapps_utils.param_sweeps.generate import generate -from geoapps_utils.run import load_ui_json_as_dict -from geoapps_utils.utils.importing import GeoAppsError -from geoh5py.data import FilenameData -from geoh5py.groups import ContainerGroup, SimPEGGroup -from geoh5py.objects import DrapeModel, PotentialElectrode, Surface -from geoh5py.shared.utils import fetch_active_workspace -from geoh5py.ui_json import InputFile -from geoh5py.workspace import Workspace - -from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.options import BaseInversionOptions -from simpeg_drivers.utils.utils import active_from_xyz, drape_to_octree - - -class LineSweepDriver(SweepDriver, InversionDriver): - """Line Sweep driver for batch 2D forward and inversion drivers.""" - - _params_class = SweepParams - - def __init__(self, params): - self.batch2d_params = params - self.cleanup = params.file_control.cleanup - - params = self.setup_params() - params.inversion_type = self.batch2d_params.inversion_type - super().__init__(params) - - @property - def out_group(self): - """The SimPEGGroup""" - return self._out_group - - @out_group.setter - def out_group(self, value: SimPEGGroup): - if not isinstance(value, SimPEGGroup): - raise TypeError("Output group must be a SimPEGGroup.") - - self.batch2d_params.out_group = value - self.batch2d_params.update_out_group_options() - self._out_group = value - - def validate_out_group(self, out_group: SimPEGGroup | None) -> SimPEGGroup: - """ - Validate or create a SimPEGGroup to store results. - - :param out_group: Output group from selection. - """ - if isinstance(out_group, SimPEGGroup): - return out_group - - with fetch_active_workspace(self.workspace, mode="r+"): - out_group = self.workspace.get_entity(self.batch2d_params.title)[0] - if out_group is None: - out_group = SimPEGGroup.create( - self.workspace, name=self.batch2d_params.title - ) - - return out_group - - def run(self): - """ - Run the line sweep driver. - - TODO: Add parallelization on GEOPY-2490 - """ - with fetch_active_workspace(self.workspace, mode="r+"): - if not isinstance(self.out_group, SimPEGGroup): - raise GeoAppsError( - f"Output group should be a valid SimPEGGroup, received: {type(self.out_group)}." - ) - - lookup = self.get_lookup() - self.write_files(lookup) - - for name, trial in lookup.items(): - file_path = Path(self.working_directory) / f"{name}.ui.json" - if trial["status"] == "complete": - continue - - trial["status"] = "processing" - self.update_lookup(lookup) - params_dict = load_ui_json_as_dict(file_path) - driver = self.driver_class_from_name( - params_dict["inversion_type"], - forward_only=params_dict["forward_only"], - ) - driver.start(file_path) - - trial["status"] = "complete" - self.update_lookup(lookup) - - self.collect_results() - - if self.cleanup: - self.file_cleanup() - - def setup_params(self): - h5_file_path = Path(self.workspace.h5file).resolve() - ui_json_path = h5_file_path.parent / ( - re.sub(r"\.ui$", "", h5_file_path.stem) + ".ui.json" - ) - if not (ui_json_path).is_file(): - with fetch_active_workspace(self.workspace): - self.batch2d_params.write_ui_json( - path=h5_file_path.parent / ui_json_path.name - ) - generate( - ui_json_path, - parameters=["line_id"], - update_values={"conda_environment": "simpeg_drivers"}, - ) - ifile = InputFile.read_ui_json( - h5_file_path.parent - / (re.sub(r"\.ui$", "", h5_file_path.stem) + "_sweep.ui.json") - ) - with fetch_active_workspace(self.workspace): - lines = self.batch2d_params.line_selection.line_object.values - ifile.data["line_id_start"] = int(lines.min()) - ifile.data["line_id_end"] = int(lines.max()) - ifile.data["line_id_n"] = len(np.unique(lines)) - sweep_params = SweepParams.build(ifile) - sweep_params.geoh5 = self.workspace - return sweep_params - - def file_cleanup(self): - """Remove files associated with the parameter sweep.""" - path = Path(self.workspace.h5file).parent - with open(path / "lookup.json", encoding="utf8") as f: - files = list(json.load(f)) - - files = [f"{f}.ui.json" for f in files] + [f"{f}.ui.geoh5" for f in files] - files += ["lookup.json"] - files += [f.name for f in path.glob("*_sweep.ui.json")] - - for file in files: - (path / file).unlink(missing_ok=True) - - @staticmethod - def line_files(path: str | Path): - with open(Path(path) / "lookup.json", encoding="utf8") as file: - line_files = {v["line_id"]: k for k, v in json.load(file).items()} - return line_files - - def collect_results(self): - path = Path(self.workspace.h5file).parent - files = LineSweepDriver.line_files(str(path)) - line_ids = self.batch2d_params.line_selection.line_object.values - data = {} - drape_models = [] - - out_lines = [] - log_lines = [] - for line in np.unique(line_ids): - with Workspace(f"{path / files[line]}.ui.geoh5") as ws: - out_group = next( - group for group in ws.groups if isinstance(group, SimPEGGroup) - ) - run_group = ContainerGroup.create( - self.workspace, name=f"Line {line}", parent=self.out_group - ) - local_simpeg_group = out_group.copy( - parent=run_group, copy_children=True, copy_relatives=True - ) - survey = next( - child - for child in local_simpeg_group.children - if isinstance(child, PotentialElectrode) - ) - line_data = survey.get_entity( - self.batch2d_params.line_selection.line_object.name - ) - - if not line_data: - raise GeoAppsError(f"Line {line} not found in {survey.name}") - - line_indices = line_ids == line - data = self.collect_line_data(survey, line_indices, data) - mesh = next( - child - for child in local_simpeg_group.children - if isinstance(child, DrapeModel) - ) - filedata = [ - k for k in out_group.children if isinstance(k, FilenameData) - ] - for fdat in filedata: - if ".out" in fdat.name: - out_lines += [f"Line {line} from file {files[line]}\n"] - out_lines += fdat.file_bytes.decode(encoding="utf8").split( - sep="\n" - ) - out_lines += ["\n"] - - if ".log" in fdat.name: - log_lines += fdat.file_bytes.decode(encoding="utf8").split( - sep="\n" - ) - log_lines += ["\n"] - - fdat.copy(parent=out_group) - - drape_models.append(mesh) - - # Write new log files to disk - with open(ws.h5file.parent / "SimPEG.out", "w", encoding="utf8") as f: - f.write("".join(out_lines)) - - with open(ws.h5file.parent / "SimPEG.log", "w", encoding="utf8") as f: - f.write("".join(log_lines)) - - self.batch2d_params.data_object.add_data(data) - - if self.batch2d_params.mesh is None: - return - - # interpolate drape model children common to all drape models into octree - active = active_from_xyz( - self.batch2d_params.mesh, - self.inversion_topography.locations, - triangulation=getattr( - self.batch2d_params.active_cells.topography_object, "cells", None - ), - ) - common_children = set.intersection( - *[{c.name for c in d.children} for d in drape_models] - ) - children = {n: [n] * len(drape_models) for n in common_children} - octree_model = drape_to_octree( - self.batch2d_params.mesh, drape_models, children, active, method="nearest" - ) - - # interpolate last iterations for each drape model into octree - iter_children = [ - [c.name for c in m.children if "iteration" in c.name.lower()] - for m in drape_models - ] - if any(iter_children): - iter_numbers = [ - [int(re.findall(r"\d+", n)[0]) for n in k] for k in iter_children - ] - last_iterations = [np.where(k == np.max(k))[0][0] for k in iter_numbers] - label = re.sub(r"\d+", "final", iter_children[0][0]) - children = { - label: [c[last_iterations[i]] for i, c in enumerate(iter_children)] - } - octree_model = drape_to_octree( - self.batch2d_params.mesh, - drape_models, - children, - active, - method="nearest", - ) - - octree_model.copy(parent=self.out_group) - - def collect_line_data(self, survey, line_indices, data): - """ - Fill chunks of values from one line - """ - for name in survey.get_data_list(): - if "Iteration" not in name: - continue - - child = survey.get_entity(name)[0] - if name not in data: - data[name] = {"values": np.zeros_like(line_indices) * np.nan} - - data[name]["values"][line_indices] = child.values - - if isinstance(self.batch2d_params, BaseInversionOptions): - label = re.sub(r"\d+", "final", name) - - if label not in data: - data[label] = {"values": np.zeros_like(line_indices) * np.nan} - - data[label]["values"][line_indices] = child.values - - return data - - @property - def workspace(self): - """Application workspace.""" - return self.batch2d_params.geoh5 diff --git a/simpeg_drivers/options.py b/simpeg_drivers/options.py index ec6028a7..3d53cde0 100644 --- a/simpeg_drivers/options.py +++ b/simpeg_drivers/options.py @@ -512,12 +512,12 @@ class LineSelectionOptions(BaseModel): arbitrary_types_allowed=True, ) line_id: int | None = None - line_object: ReferencedData + line_object: IntegerData | ReferencedData | None = None @field_validator("line_object", mode="before") @classmethod def validate_cell_association(cls, value): - if value.association is not DataAssociationEnum.CELL: + if value and value.association is not DataAssociationEnum.CELL: raise ValueError("Line identifier must be associated with cells.") return value diff --git a/simpeg_drivers/utils/nested.py b/simpeg_drivers/utils/nested.py index cf67a2a4..4e4823e3 100644 --- a/simpeg_drivers/utils/nested.py +++ b/simpeg_drivers/utils/nested.py @@ -64,10 +64,10 @@ def create_nested_mesh( if base_mesh.dim == 1: return base_mesh - return create_nested_2dmesh(survey, base_mesh) + return create_nested_2d_tensor(survey, base_mesh) -def create_nested_2dmesh( +def create_nested_2d_tensor( survey: BaseSurvey, base_mesh: TensorMesh, ) -> TensorMesh: @@ -275,7 +275,7 @@ def create_simulation( if isinstance(simulation, BaseEM1DSimulation): local_mesh = simulation.layers_mesh - actives = np.ones(simulation.layers_mesh.n_cells, dtype=bool) + local_actives = np.ones(simulation.layers_mesh.n_cells, dtype=bool) model_slice = np.arange( indices, simulation.mesh.n_cells, simulation.mesh.shape_cells[0] )[::-1] @@ -303,26 +303,28 @@ def create_simulation( 3 if getattr(simulation, "model_type", None) == "vector" else 1 ), ) - actives = mapping.local_active + local_actives = mapping.local_active # For DCIP-2D else: - actives_2d = simulation.active_cells.reshape( - simulation.mesh.shape_cells, order="F" - ) - local_active = actives_2d[simulation.mesh.parts == tile_id, :] - actives = local_active.flatten(order="F") - # Create a projection from the global active cells to the local active cells + active_mesh_part = np.isin( + simulation.mesh.parts, np.unique(local_survey.line_ids) + ) n_actives = simulation.active_cells.sum() activate_ind = np.zeros(simulation.mesh.n_cells, dtype=int) activate_ind[np.where(simulation.active_cells)[0]] = np.arange(n_actives) activate_ind = activate_ind.reshape(simulation.mesh.shape_cells, order="F") - local_active_ind = activate_ind[ - simulation.mesh.parts == tile_id, : - ].flatten(order="F")[actives] - mapping = maps.Projection(n_actives, local_active_ind) - n_actives = int(actives.sum()) + actives_2d = simulation.active_cells.reshape( + simulation.mesh.shape_cells, order="F" + ) + local_actives = actives_2d[active_mesh_part, :].flatten(order="F") + local_active_inds = activate_ind[active_mesh_part, :].flatten(order="F")[ + local_actives + ] + mapping = maps.Projection(n_actives, local_active_inds) + + n_actives = int(local_actives.sum()) if getattr(simulation, "_chiMap", None) is not None: if simulation.model_type == "vector": kwargs["chiMap"] = maps.IdentityMap(nP=n_actives * 3) @@ -330,22 +332,24 @@ def create_simulation( else: kwargs["chiMap"] = maps.IdentityMap(nP=n_actives) - kwargs["active_cells"] = actives + kwargs["active_cells"] = local_actives if getattr(simulation, "_rhoMap", None) is not None: kwargs["rhoMap"] = maps.IdentityMap(nP=n_actives) - kwargs["active_cells"] = actives + kwargs["active_cells"] = local_actives if getattr(simulation, "_sigmaMap", None) is not None: kwargs["sigmaMap"] = maps.ExpMap(local_mesh) * maps.InjectActiveCells( - local_mesh, actives, value_inactive=np.log(1e-8) + local_mesh, local_actives, value_inactive=np.log(1e-8) ) if getattr(simulation, "_etaMap", None) is not None: - kwargs["etaMap"] = maps.InjectActiveCells(local_mesh, actives, value_inactive=0) + kwargs["etaMap"] = maps.InjectActiveCells( + local_mesh, local_actives, value_inactive=0 + ) proj = maps.InjectActiveCells( local_mesh, - actives, + local_actives, value_inactive=1e-8, ) kwargs["sigma"] = proj * mapping * simulation.sigma[simulation.active_cells] @@ -439,6 +443,10 @@ def create_survey( slice_inds = slice_from_ordering(survey, indices, channel=channel) new_survey.ordering = survey.ordering[slice_inds, :] + + if hasattr(survey, "line_ids"): + new_survey.line_ids = survey.line_ids[slice_inds] + if hasattr(survey, "dobs") and survey.dobs is not None: # Return the subset of data that belongs to the tile diff --git a/simpeg_drivers/utils/surveys.py b/simpeg_drivers/utils/surveys.py index 78da3ff8..d01c840a 100644 --- a/simpeg_drivers/utils/surveys.py +++ b/simpeg_drivers/utils/surveys.py @@ -16,6 +16,7 @@ from geoapps_utils.utils.numerical import traveling_salesman from geoh5py import Workspace from geoh5py.objects import PotentialElectrode +from scipy.sparse import csgraph, csr_matrix from scipy.spatial import cKDTree from simpeg.survey import BaseSurvey @@ -82,7 +83,7 @@ def compute_alongline_distance(points: np.ndarray, ordered: bool = True): return distances -def extract_dcip_survey( +def copy_potentials_from_mask( workspace: Workspace, survey: PotentialElectrode, cell_mask: np.ndarray ): """ @@ -145,6 +146,29 @@ def get_unique_locations(survey: BaseSurvey) -> np.ndarray: return np.unique(locations, axis=0) +def get_parts_from_electrodes(survey: PotentialElectrode) -> np.ndarray: + """ + Get part numbers from a survey containing PotentialElectrode objects. + + :param survey: PotentialElectrode survey object. + + :return: Array of part numbers corresponding to each cell in the survey. + """ + edge_array = csr_matrix( + ( + np.ones(survey.n_cells * 2), + ( + np.kron(survey.cells[:, 0], [1, 1]), + survey.cells.flatten(), + ), + ), + shape=(survey.n_vertices, survey.n_vertices), + ) + + connections = csgraph.connected_components(edge_array)[1] + return connections[survey.cells[:, 0]] + + def compute_em_projections(locations, simulation): """ Pre-compute projections for the receivers for efficiency. diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 2520d97f..407bc13c 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -310,7 +310,7 @@ def drape_2_tensor(drape_model: DrapeModel, return_sorting: bool = False) -> tup # Skip indices for ghost points count = -1 - part = 0 + part = 1 parts = [] cell_widths = [] section = [] diff --git a/tests/run_tests/driver_2d_rotated_gradients_test.py b/tests/run_tests/driver_2d_rotated_gradients_test.py deleted file mode 100644 index 7c3df653..00000000 --- a/tests/run_tests/driver_2d_rotated_gradients_test.py +++ /dev/null @@ -1,201 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -from __future__ import annotations - -from pathlib import Path - -import numpy as np -from geoapps_utils.modelling.plates import PlateModel -from geoh5py.groups.property_group import PropertyGroup -from geoh5py.workspace import Workspace - -from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( - DC2DForwardDriver, - DC2DInversionDriver, -) -from simpeg_drivers.electricals.direct_current.two_dimensions.options import ( - DC2DForwardOptions, - DC2DInversionOptions, -) -from simpeg_drivers.options import ( - DrapeModelOptions, - LineSelectionOptions, -) -from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents -from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, - ModelOptions, - SurveyOptions, - SyntheticsComponentsOptions, -) -from tests.utils.targets import check_target, get_inversion_output, get_workspace - - -# To test the full run and validate the inversion. -# Move this file out of the test directory and run. - -target_run = {"data_norm": 0.644727102942656, "phi_d": 182, "phi_m": 66.3} - - -def test_dc2d_rotated_grad_fwr_run( - tmp_path: Path, - n_electrodes=10, - n_lines=3, - refinement=(4, 6), -): - filepath = Path(tmp_path) / "inversion_test.ui.geoh5" - with get_workspace(filepath) as geoh5: - # Run the forward - components = SyntheticsComponents( - geoh5=geoh5, - options=SyntheticsComponentsOptions( - method="direct current 2d", - survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), - model=ModelOptions( - background=0.01, - anomaly=10.0, - plate=PlateModel( - strike_length=1000.0, - dip_length=150.0, - width=20.0, - origin=(50.0, 0.0, -30), - direction=90, - dip=45, - ), - ), - ), - ) - - line_selection = LineSelectionOptions( - line_object=components.survey.get_data("line_ids")[0], - line_id=101, - ) - params = DC2DForwardOptions.build( - geoh5=geoh5, - data_object=components.survey, - line_selection=line_selection, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), - starting_model=components.model, - topography_object=components.topography, - ) - fwr_driver = DC2DForwardDriver(params) - fwr_driver.run() - - -def test_dc2d_rotated_grad_run( - tmp_path: Path, - max_iterations=1, - pytest=True, -): - workpath = tmp_path / "inversion_test.ui.geoh5" - if pytest: - workpath = ( - tmp_path.parent - / "test_dc2d_rotated_grad_fwr_run0" - / "inversion_test.ui.geoh5" - ) - - with Workspace(workpath) as geoh5: - potential = geoh5.get_entity("Iteration_0_potential")[0] - components = SyntheticsComponents(geoh5) - - orig_potential = potential.values.copy() - mesh = geoh5.get_entity("Models")[0] - - # Create property group with orientation - i = np.ones(mesh.n_cells) - j = np.zeros(mesh.n_cells) - k = np.ones(mesh.n_cells) - - data_list = mesh.add_data( - { - "i": {"values": i}, - "j": {"values": j}, - "k": {"values": k}, - } - ) - pg = PropertyGroup(mesh, properties=data_list, property_group_type="3D vector") - - # Run the inverse - params = DC2DInversionOptions.build( - geoh5=geoh5, - mesh=mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), - topography_object=components.topography, - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, - ), - data_object=potential.parent, - gradient_rotation=pg, - potential_channel=potential, - potential_uncertainty=1e-3, - model_type="Resistivity (Ohm-m)", - starting_model=100.0, - reference_model=100.0, - s_norm=1.0, - x_norm=0.0, - z_norm=2.0, - max_global_iterations=max_iterations, - initial_beta=None, - initial_beta_ratio=1e0, - percentile=100, - lower_bound=0.1, - cooling_rate=1, - starting_chi_factor=1.0, - chi_factor=0.1, - ) - params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - - driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - - with Workspace(driver.params.geoh5.h5file) as run_ws: - output = get_inversion_output( - driver.params.geoh5.h5file, driver.params.out_group.uid - ) - output["data"] = orig_potential - - if pytest: - check_target(output, target_run) - nan_ind = np.isnan(run_ws.get_entity("Iteration_0_model")[0].values) - inactive_ind = run_ws.get_entity("active_cells")[0].values == 0 - assert np.all(nan_ind == inactive_ind) - - -if __name__ == "__main__": - # Full run - test_dc2d_rotated_grad_fwr_run( - Path("./"), - n_electrodes=20, - n_lines=3, - refinement=(4, 4), - ) - - test_dc2d_rotated_grad_run( - Path("./"), - max_iterations=20, - pytest=False, - ) diff --git a/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py b/tests/run_tests/driver_dc_2d_rotated_gradients_test.py similarity index 73% rename from tests/run_tests/driver_dc_b2d_rotated_gradients_test.py rename to tests/run_tests/driver_dc_2d_rotated_gradients_test.py index 761101a8..e7081b6d 100644 --- a/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py +++ b/tests/run_tests/driver_dc_2d_rotated_gradients_test.py @@ -19,26 +19,21 @@ from geoh5py.groups import PropertyGroup, SimPEGGroup from geoh5py.workspace import Workspace -from simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.driver import ( - DCBatch2DForwardDriver, - DCBatch2DInversionDriver, +from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( + DC2DForwardDriver, + DC2DInversionDriver, ) -from simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.options import ( - DCBatch2DForwardOptions, - DCBatch2DInversionOptions, -) -from simpeg_drivers.electricals.options import ( - FileControlOptions, +from simpeg_drivers.electricals.direct_current.two_dimensions.options import ( + DC2DForwardOptions, + DC2DInversionOptions, ) from simpeg_drivers.options import ( DrapeModelOptions, - LineSelectionOptions, ) from simpeg_drivers.utils.synthetics.driver import ( SyntheticsComponents, ) from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, ModelOptions, SurveyOptions, SyntheticsComponentsOptions, @@ -52,13 +47,17 @@ target_run = {"data_norm": 1.1067294238524659, "phi_d": 55.6, "phi_m": 7.08} -def test_dc_rotated_p3d_fwr_run( - tmp_path: Path, n_electrodes=10, n_lines=3, refinement=(4, 6) -): +def test_dc_rotated_2d_fwr_run(tmp_path: Path, n_electrodes=10, n_lines=3): opts = SyntheticsComponentsOptions( - method="direct current pseudo 3d", + method="direct current 2d", survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), + mesh=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=50.0, + expansion_factor=1.1, + vertical_padding=100.0, + ), model=ModelOptions( background=0.01, anomaly=10.0, @@ -74,29 +73,18 @@ def test_dc_rotated_p3d_fwr_run( ) with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: components = SyntheticsComponents(geoh5, options=opts) - params = DCBatch2DForwardOptions.build( + params = DC2DForwardOptions.build( geoh5=geoh5, mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), topography_object=components.topography, data_object=components.survey, starting_model=components.model, - line_selection=LineSelectionOptions( - line_object=components.survey.get_data("line_ids")[0] - ), ) - fwr_driver = DCBatch2DForwardDriver(params) + fwr_driver = DC2DForwardDriver(params) fwr_driver.run() -def test_dc_rotated_gradient_p3d_run( +def test_dc_rotated_gradient_2d_run( tmp_path: Path, max_iterations=1, pytest=True, @@ -104,13 +92,13 @@ def test_dc_rotated_gradient_p3d_run( workpath = tmp_path / "inversion_test.ui.geoh5" if pytest: workpath = ( - tmp_path.parent / "test_dc_rotated_p3d_fwr_run0" / "inversion_test.ui.geoh5" + tmp_path.parent / "test_dc_rotated_2d_fwr_run0" / "inversion_test.ui.geoh5" ) with Workspace(workpath) as geoh5: components = SyntheticsComponents(geoh5) - fwr_group = geoh5.get_entity("Direct Current (DC) 2D Batch Forward")[0] + fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] survey = fwr_group.get_entity("survey")[0] potential = survey.get_data("Iteration_0_potential")[0] # Create property group with orientation @@ -130,7 +118,7 @@ def test_dc_rotated_gradient_p3d_run( ) # Run the inverse - params = DCBatch2DInversionOptions.build( + params = DC2DInversionOptions.build( geoh5=geoh5, mesh=components.mesh, drape_model=DrapeModelOptions( @@ -146,9 +134,6 @@ def test_dc_rotated_gradient_p3d_run( gradient_rotation=pg, potential_channel=potential, potential_uncertainty=1e-3, - line_selection=LineSelectionOptions( - line_object=potential.parent.get_entity("line_ids")[0] - ), starting_model=1e-2, reference_model=1e-2, s_norm=0.0, @@ -160,11 +145,10 @@ def test_dc_rotated_gradient_p3d_run( percentile=100, upper_bound=10, cooling_rate=1, - file_control=FileControlOptions(cleanup=False), ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - DCBatch2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) + DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) basepath = workpath.parent with open(basepath / "lookup.json", encoding="utf8") as f: @@ -187,10 +171,8 @@ def test_dc_rotated_gradient_p3d_run( if __name__ == "__main__": # Full run - test_dc_rotated_p3d_fwr_run( - Path("./"), n_electrodes=20, n_lines=3, refinement=(4, 4) - ) - test_dc_rotated_gradient_p3d_run( + test_dc_rotated_2d_fwr_run(Path("./"), n_electrodes=20, n_lines=3) + test_dc_rotated_gradient_2d_run( Path("./"), max_iterations=20, pytest=False, diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index 4d3d80bc..d60f34bf 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -13,7 +13,6 @@ from pathlib import Path -import numpy as np from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( @@ -32,7 +31,6 @@ SyntheticsComponents, ) from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, ModelOptions, SurveyOptions, SyntheticsComponentsOptions, @@ -43,87 +41,129 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.6072705410748107, "phi_d": 681, "phi_m": 95.6} +target_run = {"data_norm": 1.101767837151429, "phi_d": 2210, "phi_m": 21.4} -def test_dc_2d_fwr_run( +def test_dc_p3d_fwr_run( tmp_path: Path, n_electrodes=10, n_lines=3, - refinement=(4, 6), ): # Run the forward opts = SyntheticsComponentsOptions( method="direct current 2d", survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), - model=ModelOptions(background=0.01, anomaly=10), + mesh=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=50.0, + expansion_factor=1.1, + vertical_padding=100.0, + ), + model=ModelOptions(background=0.01, anomaly=10.0), ) with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: - components = SyntheticsComponents(geoh5, options=opts) - - line_selection = LineSelectionOptions( - line_object=components.survey.get_data("line_ids")[0], - line_id=101, - ) + components = SyntheticsComponents(geoh5=geoh5, options=opts) params = DC2DForwardOptions.build( geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, data_object=components.survey, - line_selection=line_selection, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), starting_model=components.model, - topography_object=components.topography, + line_selection=LineSelectionOptions(), ) fwr_driver = DC2DForwardDriver(params) fwr_driver.run() -def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): +def test_dc_p3d_run( + tmp_path: Path, + max_iterations=1, + pytest=True, +): workpath = tmp_path / "inversion_test.ui.geoh5" if pytest: - workpath = tmp_path.parent / "test_dc_2d_fwr_run0" / "inversion_test.ui.geoh5" + workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: - potential = geoh5.get_entity("Iteration_0_potential")[0] components = SyntheticsComponents(geoh5) + fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] + survey = fwr_group.get_entity("survey")[0] + potential = survey.get_data("Iteration_0_potential")[0] # Run the inverse params = DC2DInversionOptions.build( geoh5=geoh5, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), + mesh=components.mesh, topography_object=components.topography, + data_object=potential.parent, + potential_channel=potential, + potential_uncertainty=1e-3, line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, + line_object=potential.parent.get_entity("Line IDs")[0] ), + starting_model=1e-2, + reference_model=1e-2, + s_norm=0.0, + x_norm=1.0, + z_norm=1.0, + max_global_iterations=max_iterations, + initial_beta=None, + initial_beta_ratio=10.0, + percentile=100, + upper_bound=10, + cooling_rate=1, + ) + params.write_ui_json(path=tmp_path / "Inv_run.ui.json") + + driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) + + output = get_inversion_output( + driver.params.geoh5.h5file, driver.params.out_group.uid + ) + if geoh5.open(): + output["data"] = potential.values + if pytest: + check_target(output, target_run) + + +def test_dc_single_run( + tmp_path: Path, + max_iterations=1, + pytest=True, +): + workpath = tmp_path / "inversion_test.ui.geoh5" + if pytest: + workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" + + with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) + fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] + survey = fwr_group.get_entity("survey")[0] + potential = survey.get_data("Iteration_0_potential")[0] + + # Run the inverse + params = DC2DInversionOptions.build( + geoh5=geoh5, + title="Direct Current Single 2D Inversion", + mesh=components.mesh, + topography_object=components.topography, data_object=potential.parent, potential_channel=potential, potential_uncertainty=1e-3, - model_type="Resistivity (Ohm-m)", - starting_model=100.0, - reference_model=100.0, + line_selection=LineSelectionOptions( + line_object=potential.parent.get_entity("Line IDs")[0], line_id=2 + ), + starting_model=1e-2, + reference_model=1e-2, s_norm=0.0, x_norm=1.0, z_norm=1.0, max_global_iterations=max_iterations, initial_beta=None, - initial_beta_ratio=1e0, + initial_beta_ratio=10.0, percentile=100, - lower_bound=0.1, + upper_bound=10, cooling_rate=1, ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") @@ -134,20 +174,19 @@ def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): driver.params.geoh5.h5file, driver.params.out_group.uid ) if geoh5.open(): - output["data"] = potential.values[np.isfinite(potential.values)] + output["data"] = potential.values if pytest: check_target(output, target_run) if __name__ == "__main__": # Full run - test_dc_2d_fwr_run( + test_dc_p3d_fwr_run( Path("./"), n_electrodes=20, n_lines=3, - refinement=(4, 4), ) - test_dc_2d_run( + test_dc_p3d_run( Path("./"), max_iterations=20, pytest=False, diff --git a/tests/run_tests/driver_dc_b2d_test.py b/tests/run_tests/driver_dc_b2d_test.py deleted file mode 100644 index 940bab6e..00000000 --- a/tests/run_tests/driver_dc_b2d_test.py +++ /dev/null @@ -1,143 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from pathlib import Path - -from geoh5py.workspace import Workspace - -from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( - DC2DForwardDriver, - DC2DInversionDriver, -) -from simpeg_drivers.electricals.direct_current.two_dimensions.options import ( - DC2DForwardOptions, - DC2DInversionOptions, -) -from simpeg_drivers.options import ( - DrapeModelOptions, - LineSelectionOptions, -) -from simpeg_drivers.utils.synthetics.driver import ( - SyntheticsComponents, -) -from simpeg_drivers.utils.synthetics.options import ( - ModelOptions, - SurveyOptions, - SyntheticsComponentsOptions, -) -from tests.utils.targets import check_target, get_inversion_output, get_workspace - - -# To test the full run and validate the inversion. -# Move this file out of the test directory and run. - -target_run = {"data_norm": 1.101767837151429, "phi_d": 2210, "phi_m": 21.4} - - -def test_dc_p3d_fwr_run( - tmp_path: Path, - n_electrodes=10, - n_lines=3, -): - # Run the forward - opts = SyntheticsComponentsOptions( - method="direct current 2d", - survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=50.0, - expansion_factor=1.1, - vertical_padding=100.0, - ), - model=ModelOptions(background=0.01, anomaly=10.0), - ) - with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: - components = SyntheticsComponents(geoh5=geoh5, options=opts) - params = DC2DForwardOptions.build( - geoh5=geoh5, - mesh=components.mesh, - topography_object=components.topography, - data_object=components.survey, - starting_model=components.model, - line_selection=LineSelectionOptions( - line_object=components.survey.get_data("line_ids")[0] - ), - ) - fwr_driver = DC2DForwardDriver(params) - fwr_driver.run() - - -def test_dc_p3d_run( - tmp_path: Path, - max_iterations=1, - pytest=True, -): - workpath = tmp_path / "inversion_test.ui.geoh5" - if pytest: - workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" - - with Workspace(workpath) as geoh5: - components = SyntheticsComponents(geoh5) - fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] - survey = fwr_group.get_entity("survey")[0] - potential = survey.get_data("Iteration_0_potential")[0] - - # Run the inverse - params = DC2DInversionOptions.build( - geoh5=geoh5, - mesh=components.mesh, - topography_object=components.topography, - data_object=potential.parent, - potential_channel=potential, - potential_uncertainty=1e-3, - line_selection=LineSelectionOptions( - line_object=potential.parent.get_entity("line_ids")[0] - ), - starting_model=1e-2, - reference_model=1e-2, - s_norm=0.0, - x_norm=1.0, - z_norm=1.0, - max_global_iterations=max_iterations, - initial_beta=None, - initial_beta_ratio=10.0, - percentile=100, - upper_bound=10, - cooling_rate=1, - ) - params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - - driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - - output = get_inversion_output( - driver.params.geoh5.h5file, driver.params.out_group.uid - ) - if geoh5.open(): - output["data"] = potential.values - if pytest: - check_target(output, target_run) - - -if __name__ == "__main__": - # Full run - test_dc_p3d_fwr_run( - Path("./"), - n_electrodes=20, - n_lines=3, - ) - test_dc_p3d_run( - Path("./"), - max_iterations=20, - pytest=False, - ) diff --git a/tests/run_tests/driver_ip_b2d_test.py b/tests/run_tests/driver_ip_b2d_test.py deleted file mode 100644 index 612610a2..00000000 --- a/tests/run_tests/driver_ip_b2d_test.py +++ /dev/null @@ -1,180 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - -from __future__ import annotations - -import json -from pathlib import Path - -from geoh5py.groups import SimPEGGroup -from geoh5py.workspace import Workspace - -from simpeg_drivers.electricals.induced_polarization.pseudo_three_dimensions.driver import ( - IPBatch2DForwardDriver, - IPBatch2DInversionDriver, -) -from simpeg_drivers.electricals.induced_polarization.pseudo_three_dimensions.options import ( - IPBatch2DForwardOptions, - IPBatch2DInversionOptions, -) -from simpeg_drivers.electricals.options import ( - FileControlOptions, -) -from simpeg_drivers.options import ( - ActiveCellsOptions, - DrapeModelOptions, - LineSelectionOptions, -) -from simpeg_drivers.utils.synthetics.driver import ( - SyntheticsComponents, -) -from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, - ModelOptions, - SurveyOptions, - SyntheticsComponentsOptions, -) -from tests.utils.targets import check_target, get_inversion_output, get_workspace - - -# To test the full run and validate the inversion. -# Move this file out of the test directory and run. - -target_run = {"data_norm": 0.09310413606088193, "phi_d": 217000, "phi_m": 5.09e-08} - - -def test_ip_p3d_fwr_run( - tmp_path: Path, - n_electrodes=10, - n_lines=3, - refinement=(4, 6), -): - # Run the forward - opts = SyntheticsComponentsOptions( - method="induced polarization pseudo 3d", - survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), - model=ModelOptions(background=1e-6, anomaly=1e-1), - ) - with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: - components = SyntheticsComponents(geoh5, options=opts) - - params = IPBatch2DForwardOptions.build( - geoh5=geoh5, - mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=100.0, - vertical_padding=100.0, - ), - active_cells=ActiveCellsOptions( - topography_object=components.topography, - ), - data_object=components.survey, - conductivity_model=1e-2, - starting_model=components.model, - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0] - ), - ) - - fwr_driver = IPBatch2DForwardDriver(params) - fwr_driver.run() - - -def test_ip_p3d_run( - tmp_path, - max_iterations=1, - pytest=True, -): - workpath = tmp_path / "inversion_test.ui.geoh5" - if pytest: - workpath = tmp_path.parent / "test_ip_p3d_fwr_run0" / "inversion_test.ui.geoh5" - - with Workspace(workpath) as geoh5: - components = SyntheticsComponents(geoh5) - fwr_group = geoh5.get_entity("Induced Polarization (IP) 2D Batch Forward")[0] - survey = fwr_group.get_entity("survey")[0] - chargeability = survey.get_data("Iteration_0_chargeability")[0] - - # Run the inverse - params = IPBatch2DInversionOptions.build( - geoh5=geoh5, - mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), - topography_object=components.topography, - data_object=chargeability.parent, - chargeability_channel=chargeability, - chargeability_uncertainty=2e-4, - line_selection=LineSelectionOptions( - line_object=chargeability.parent.get_entity("line_ids")[0], - ), - conductivity_model=1e-2, - starting_model=1e-6, - reference_model=1e-6, - s_norm=0.0, - x_norm=0.0, - z_norm=0.0, - length_scale_x=1.0, - length_scale_z=1.0, - max_global_iterations=max_iterations, - initial_beta=None, - initial_beta_ratio=1e0, - percentile=100, - upper_bound=0.1, - cooling_rate=1, - file_control=FileControlOptions(cleanup=False), - ) - params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - - IPBatch2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - - basepath = workpath.parent - with open(basepath / "lookup.json", encoding="utf8") as f: - lookup = json.load(f) - middle_line_id = next(k for k, v in lookup.items() if v["line_id"] == 101) - - with Workspace(basepath / f"{middle_line_id}.ui.geoh5", mode="r") as workspace: - middle_inversion_group = next( - k for k in workspace.groups if isinstance(k, SimPEGGroup) - ) - - output = get_inversion_output( - basepath / f"{middle_line_id}.ui.geoh5", middle_inversion_group.uid - ) - if geoh5.open(): - output["data"] = chargeability.values - if pytest: - check_target(output, target_run) - - -if __name__ == "__main__": - # Full run - test_ip_p3d_fwr_run( - Path("./"), - n_electrodes=20, - n_lines=3, - refinement=(4, 4), - ) - test_ip_p3d_run( - Path("./"), - max_iterations=20, - pytest=False, - ) From 1d66a57f82e0b4d96da354bffa20abf38294725e Mon Sep 17 00:00:00 2001 From: dominiquef Date: Sun, 8 Mar 2026 10:13:27 -0700 Subject: [PATCH 06/24] Remove ghost cells from model save. Add back app res --- simpeg_drivers/components/data.py | 10 ++++++++++ .../components/factories/directives_factory.py | 2 +- simpeg_drivers/components/models.py | 2 ++ 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/simpeg_drivers/components/data.py b/simpeg_drivers/components/data.py index 75f5e3c9..345b21b4 100644 --- a/simpeg_drivers/components/data.py +++ b/simpeg_drivers/components/data.py @@ -344,6 +344,16 @@ def create_survey(self): ) survey.cells = self.entity.cells + observed = self.entity.get_data("Observed_potential") + if observed: + self.entity.add_data( + { + "Observed_apparent_resistivity": { + "values": survey.apparent_resistivity * observed[0].values + } + } + ) + if "induced polarization" in self.params.inversion_type: survey.cells = self.entity.cells diff --git a/simpeg_drivers/components/factories/directives_factory.py b/simpeg_drivers/components/factories/directives_factory.py index c23c243d..69c7c84c 100644 --- a/simpeg_drivers/components/factories/directives_factory.py +++ b/simpeg_drivers/components/factories/directives_factory.py @@ -436,7 +436,7 @@ def assemble_keyword_arguments( if self.params.models.model_type == ModelTypeEnum.resistivity: kwargs["transforms"].append(lambda x: 1 / x) - if "1d" in self.factory_type: + if "1d" in self.factory_type or "2d" in self.factory_type: ghosts = ( np.squeeze(np.asarray(inversion_object.permutation.sum(axis=0))) == 0 ) diff --git a/simpeg_drivers/components/models.py b/simpeg_drivers/components/models.py index efb0dc22..d3b8a9fb 100644 --- a/simpeg_drivers/components/models.py +++ b/simpeg_drivers/components/models.py @@ -643,6 +643,8 @@ def _get_value(self, model: float | NumericData) -> np.ndarray: elif isinstance(model, int | float): nc = self.driver.inversion_mesh.mesh.n_cells model *= np.ones(nc) + elif isinstance(model, np.ndarray): + model = (self.driver.inversion_mesh.permutation @ model).astype(model.dtype) return model From b2cfc62119a900386a7cb3c78231ed64e90c50a0 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Mon, 9 Mar 2026 10:51:50 -0700 Subject: [PATCH 07/24] Avoid zero division --- simpeg_drivers/electricals/base_2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 8906e725..0b7c6b7b 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -128,7 +128,7 @@ def normalize_vertically(poles: np.ndarray, relief: float) -> np.ndarray: """ min_poles_z = poles[:, 2].min() poles[:, 2] -= min_poles_z - poles[:, 2] *= relief / poles[:, 2].max() + poles[:, 2] *= relief / np.maximum(poles[:, 2].max(), 1e-3) # Shift back vertically poles[:, 2] += min_poles_z From d1e5495ed2e6a341229e21e4dba382c7d94d1db9 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Mon, 9 Mar 2026 15:59:30 -0700 Subject: [PATCH 08/24] Update 2d tests --- .../driver_dc_2d_rotated_gradients_test.py | 53 ++++++------- tests/run_tests/driver_dc_2d_test.py | 78 +++++++++++++------ tests/run_tests/driver_ip_2d_test.py | 51 +++++++----- 3 files changed, 107 insertions(+), 75 deletions(-) diff --git a/tests/run_tests/driver_dc_2d_rotated_gradients_test.py b/tests/run_tests/driver_dc_2d_rotated_gradients_test.py index e7081b6d..3333ebde 100644 --- a/tests/run_tests/driver_dc_2d_rotated_gradients_test.py +++ b/tests/run_tests/driver_dc_2d_rotated_gradients_test.py @@ -11,12 +11,11 @@ from __future__ import annotations -import json from pathlib import Path import numpy as np from geoapps_utils.modelling.plates import PlateModel -from geoh5py.groups import PropertyGroup, SimPEGGroup +from geoh5py.groups import PropertyGroup from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( @@ -44,7 +43,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 1.1067294238524659, "phi_d": 55.6, "phi_m": 7.08} +target_run = {"data_norm": 10.305373769233688, "phi_d": 187000, "phi_m": 410} def test_dc_rotated_2d_fwr_run(tmp_path: Path, n_electrodes=10, n_lines=3): @@ -56,16 +55,17 @@ def test_dc_rotated_2d_fwr_run(tmp_path: Path, n_electrodes=10, n_lines=3): v_cell_size=5.0, depth_core=50.0, expansion_factor=1.1, - vertical_padding=100.0, + vertical_padding=200.0, + horizontal_padding=200.0, ), model=ModelOptions( - background=0.01, - anomaly=10.0, + background=0.001, + anomaly=1.0, plate=PlateModel( strike_length=1000.0, - dip_length=150.0, + dip_length=50.0, width=20.0, - origin=(0.0, 0.0, -50), + origin=(0.0, 0.0, 0.0), direction=90, dip=45, ), @@ -101,6 +101,13 @@ def test_dc_rotated_gradient_2d_run( fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] survey = fwr_group.get_entity("survey")[0] potential = survey.get_data("Iteration_0_potential")[0] + uncertainties = survey.add_data( + { + "Uncertainties": { + "values": np.abs(potential.values) * 0.05 + 1e-4, + } + } + ) # Create property group with orientation dip = np.ones(components.mesh.n_cells) * 45 azimuth = np.ones(components.mesh.n_cells) * 90 @@ -121,47 +128,31 @@ def test_dc_rotated_gradient_2d_run( params = DC2DInversionOptions.build( geoh5=geoh5, mesh=components.mesh, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), topography_object=components.topography, data_object=potential.parent, gradient_rotation=pg, potential_channel=potential, - potential_uncertainty=1e-3, - starting_model=1e-2, - reference_model=1e-2, + potential_uncertainty=uncertainties, + starting_model=1e-3, + reference_model=1e-3, s_norm=0.0, x_norm=0.0, z_norm=0.0, + length_scale_z=0.1, max_global_iterations=max_iterations, initial_beta=None, initial_beta_ratio=10.0, percentile=100, upper_bound=10, cooling_rate=1, + sens_wts_threshold=1.0, ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - - basepath = workpath.parent - with open(basepath / "lookup.json", encoding="utf8") as f: - lookup = json.load(f) - middle_line_id = next(k for k, v in lookup.items() if v["line_id"] == 101) - - with Workspace(basepath / f"{middle_line_id}.ui.geoh5", mode="r") as workspace: - middle_inversion_group = next( - k for k in workspace.groups if isinstance(k, SimPEGGroup) - ) + driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) output = get_inversion_output( - basepath / f"{middle_line_id}.ui.geoh5", middle_inversion_group.uid + driver.params.geoh5.h5file, driver.params.out_group.uid ) if geoh5.open(): output["data"] = potential.values diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index d60f34bf..782c16b4 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -13,6 +13,8 @@ from pathlib import Path +import numpy as np +from geoapps_utils.modelling.plates import PlateModel from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( @@ -41,10 +43,10 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 1.101767837151429, "phi_d": 2210, "phi_m": 21.4} +target_run = {"data_norm": 11.14351536256954, "phi_d": 6360, "phi_m": 245} -def test_dc_p3d_fwr_run( +def test_dc_2d_fwr_run( tmp_path: Path, n_electrodes=10, n_lines=3, @@ -58,9 +60,21 @@ def test_dc_p3d_fwr_run( v_cell_size=5.0, depth_core=50.0, expansion_factor=1.1, - vertical_padding=100.0, + vertical_padding=200.0, + horizontal_padding=200.0, + ), + model=ModelOptions( + background=0.001, + anomaly=1.0, + plate=PlateModel( + strike_length=1000.0, + dip_length=20.0, + width=20.0, + origin=(0.0, 0.0, 0.0), + direction=90, + dip=90, + ), ), - model=ModelOptions(background=0.01, anomaly=10.0), ) with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: components = SyntheticsComponents(geoh5=geoh5, options=opts) @@ -76,21 +90,27 @@ def test_dc_p3d_fwr_run( fwr_driver.run() -def test_dc_p3d_run( +def test_dc_2d_run( tmp_path: Path, max_iterations=1, pytest=True, ): workpath = tmp_path / "inversion_test.ui.geoh5" if pytest: - workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" + workpath = tmp_path.parent / "test_dc_2d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: components = SyntheticsComponents(geoh5) fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] survey = fwr_group.get_entity("survey")[0] potential = survey.get_data("Iteration_0_potential")[0] - + uncertainties = survey.add_data( + { + "Uncertainties": { + "values": np.abs(potential.values) * 0.05 + 1e-4, + } + } + ) # Run the inverse params = DC2DInversionOptions.build( geoh5=geoh5, @@ -98,12 +118,12 @@ def test_dc_p3d_run( topography_object=components.topography, data_object=potential.parent, potential_channel=potential, - potential_uncertainty=1e-3, + potential_uncertainty=uncertainties, line_selection=LineSelectionOptions( line_object=potential.parent.get_entity("Line IDs")[0] ), - starting_model=1e-2, - reference_model=1e-2, + starting_model=1e-3, + reference_model=1e-3, s_norm=0.0, x_norm=1.0, z_norm=1.0, @@ -134,14 +154,20 @@ def test_dc_single_run( ): workpath = tmp_path / "inversion_test.ui.geoh5" if pytest: - workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" + workpath = tmp_path.parent / "test_dc_2d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: components = SyntheticsComponents(geoh5) fwr_group = geoh5.get_entity("Direct Current 2D Forward")[0] survey = fwr_group.get_entity("survey")[0] potential = survey.get_data("Iteration_0_potential")[0] - + uncertainties = survey.add_data( + { + "Uncertainties": { + "values": np.abs(potential.values) * 0.05 + 1e-4, + } + } + ) # Run the inverse params = DC2DInversionOptions.build( geoh5=geoh5, @@ -150,12 +176,12 @@ def test_dc_single_run( topography_object=components.topography, data_object=potential.parent, potential_channel=potential, - potential_uncertainty=1e-3, + potential_uncertainty=uncertainties, line_selection=LineSelectionOptions( line_object=potential.parent.get_entity("Line IDs")[0], line_id=2 ), - starting_model=1e-2, - reference_model=1e-2, + starting_model=1e-3, + reference_model=1e-3, s_norm=0.0, x_norm=1.0, z_norm=1.0, @@ -168,25 +194,27 @@ def test_dc_single_run( ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") - driver = DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) + DC2DInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) - output = get_inversion_output( - driver.params.geoh5.h5file, driver.params.out_group.uid - ) - if geoh5.open(): - output["data"] = potential.values - if pytest: - check_target(output, target_run) + with Workspace(workpath) as geoh5: + inv_group = geoh5.get_entity("Direct Current Single 2D Inversion")[0] + mesh = inv_group.get_entity("mesh")[0] + model = mesh.get_entity("Iteration_1_model")[0] + + # Check that model values for lines 1 and 3 are close to the starting model (1e-3) and that line 2 has been updated. + np.testing.assert_almost_equal(np.nanmin(model.values[:2369]), 1e-3, decimal=3) + np.testing.assert_almost_equal(np.nanmin(model.values[-2368:]), 1e-3, decimal=3) + assert np.nanmax(model.values[2368:-2368]) > 1e-3 if __name__ == "__main__": # Full run - test_dc_p3d_fwr_run( + test_dc_2d_fwr_run( Path("./"), n_electrodes=20, n_lines=3, ) - test_dc_p3d_run( + test_dc_2d_run( Path("./"), max_iterations=20, pytest=False, diff --git a/tests/run_tests/driver_ip_2d_test.py b/tests/run_tests/driver_ip_2d_test.py index 38f3ef9c..d3189e7b 100644 --- a/tests/run_tests/driver_ip_2d_test.py +++ b/tests/run_tests/driver_ip_2d_test.py @@ -13,6 +13,7 @@ from pathlib import Path import numpy as np +from geoapps_utils.modelling.plates import PlateModel from geoh5py.workspace import Workspace from simpeg_drivers.electricals.induced_polarization.two_dimensions import ( @@ -23,12 +24,11 @@ IP2DForwardDriver, IP2DInversionDriver, ) -from simpeg_drivers.options import LineSelectionOptions from simpeg_drivers.utils.synthetics.driver import ( SyntheticsComponents, ) from simpeg_drivers.utils.synthetics.options import ( - MeshOptions, + DrapeModelOptions, ModelOptions, SurveyOptions, SyntheticsComponentsOptions, @@ -39,21 +39,38 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.09193270736303862, "phi_d": 211000, "phi_m": 3.11e-08} +target_run = {"data_norm": 0.1267394640365658, "phi_d": 16300, "phi_m": 6.05e-05} def test_ip_2d_fwr_run( tmp_path: Path, n_electrodes=10, n_lines=3, - refinement=(4, 6), ): # Run the forward opts = SyntheticsComponentsOptions( method="induced polarization 2d", survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), - mesh=MeshOptions(refinement=refinement), - model=ModelOptions(background=1e-6, anomaly=1e-1), + mesh=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=50.0, + expansion_factor=1.1, + vertical_padding=200.0, + horizontal_padding=200.0, + ), + model=ModelOptions( + background=1e-6, + anomaly=1e-1, + plate=PlateModel( + strike_length=1000.0, + dip_length=50.0, + width=20.0, + origin=(0.0, 0.0, 0.0), + direction=90, + dip=45, + ), + ), ) with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: components = SyntheticsComponents(geoh5, options=opts) @@ -65,10 +82,6 @@ def test_ip_2d_fwr_run( starting_model=components.model, conductivity_model=1e2, model_type="Resistivity (Ohm-m)", - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, - ), ) fwr_driver = IP2DForwardDriver(params) @@ -87,7 +100,13 @@ def test_ip_2d_run( with Workspace(workpath) as geoh5: components = SyntheticsComponents(geoh5) chargeability = geoh5.get_entity("Iteration_0_chargeability")[0] - + uncertainties = chargeability.parent.add_data( + { + "Uncertainties": { + "values": np.abs(chargeability.values) * 0.05 + 1e-4, + } + } + ) # Run the inverse params = IP2DInversionOptions.build( geoh5=geoh5, @@ -95,11 +114,7 @@ def test_ip_2d_run( topography_object=components.topography, data_object=chargeability.parent, chargeability_channel=chargeability, - chargeability_uncertainty=2e-4, - line_selection=LineSelectionOptions( - line_object=chargeability.parent.get_entity("line_ids")[0], - line_id=101, - ), + chargeability_uncertainty=uncertainties, starting_model=1e-6, reference_model=1e-6, conductivity_model=1e-2, @@ -107,8 +122,7 @@ def test_ip_2d_run( x_norm=0.0, z_norm=0.0, max_global_iterations=max_iterations, - initial_beta=None, - initial_beta_ratio=1e0, + initial_beta_ratio=1e-2, percentile=100, upper_bound=0.1, cooling_rate=1, @@ -132,7 +146,6 @@ def test_ip_2d_run( Path("./"), n_electrodes=20, n_lines=3, - refinement=(4, 4), ) test_ip_2d_run( Path("./"), From cb406253a09f11c07a369fe02f6b3aa683ce2aa9 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Mon, 9 Mar 2026 16:35:07 -0700 Subject: [PATCH 09/24] Move utilities --- simpeg_drivers/components/data.py | 3 +- simpeg_drivers/electricals/base_2d.py | 109 +----------------- simpeg_drivers/utils/surveys.py | 154 ++++++++++++++++++-------- simpeg_drivers/utils/utils.py | 30 +++-- 4 files changed, 135 insertions(+), 161 deletions(-) diff --git a/simpeg_drivers/components/data.py b/simpeg_drivers/components/data.py index 345b21b4..2966d9ae 100644 --- a/simpeg_drivers/components/data.py +++ b/simpeg_drivers/components/data.py @@ -398,8 +398,7 @@ def update_params(self, data_dict, uncert_dict): if getattr(self.params, "line_selection", None) is not None: if self.params.line_selection.line_object is None: parts = get_parts_from_electrodes(self.entity) - _, u_part = np.unique(parts, return_inverse=True) - line_ids = self.entity.add_data({"Line IDs": {"values": u_part + 1}}) + line_ids = self.entity.add_data({"Line IDs": {"values": parts + 1}}) else: line_ids = self.params.line_selection.line_object.copy( parent=self.entity, diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 0b7c6b7b..1126fced 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -11,19 +11,11 @@ from __future__ import annotations -import numpy as np -from geoh5py.data import IntegerData -from geoh5py.objects import DrapeModel -from geoh5py.shared.merging.drape_model import DrapeModelMerger from geoh5py.ui_json.ui_json import fetch_active_workspace -from geoh5py.workspace import Workspace from simpeg_drivers.components.meshes import InversionMesh from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.options import ( - DrapeModelOptions, -) -from simpeg_drivers.utils.utils import get_drape_model +from simpeg_drivers.utils.surveys import create_mesh_by_line_id class Base2DDriver(InversionDriver): @@ -53,102 +45,3 @@ def inversion_mesh(self) -> InversionMesh: ) return self._inversion_mesh - - -def create_mesh_by_line_id( - workspace: Workspace, - line_ids: IntegerData, - drape_options: DrapeModelOptions, - **object_kwargs, -) -> DrapeModel: - """ - Create a drape mesh for the dc resistivity survey lines. - - :param workspace: Workspace to create the drape mesh in. - :param line_ids: IntegerData object containing the line IDs for each vertex. - :param drape_options: DrapeModelOptions containing the parameters for the drape mesh - :param object_kwargs: Additional keyword arguments to pass to the DrapeModelMerger.create_object method. - - :return: A DrapeModel object containing the merged drape mesh for all survey lines. - """ - drape_models = [] - temp_work = Workspace() - - relief = get_max_line_relief(line_ids, drape_options.v_cell_size) - - for line_id in np.unique(line_ids.values): - poles = get_poles_by_line_id(line_ids, line_id) - poles = np.unique(poles, axis=0) - poles = normalize_vertically(poles, relief) - - drape_model = get_drape_model( - temp_work, - poles, - [ - drape_options.u_cell_size, - drape_options.v_cell_size, - ], - drape_options.depth_core, - [drape_options.horizontal_padding] * 2 - + [drape_options.vertical_padding, 1], - drape_options.expansion_factor, - ) - drape_models.append(drape_model) - - entity = DrapeModelMerger.create_object(workspace, drape_models, **object_kwargs) - - return entity - - -def get_max_line_relief(line_ids, z_cell_size) -> float: - """ - Get the maximum relief across all survey lines, rounded to the nearest cell thickness. - - :param line_ids: IntegerData object containing the line IDs for each vertex. - :param z_cell_size: Cell size in the vertical direction for the drape mesh. - """ - max_relief = 0 - for line_id in np.unique(line_ids.values): - poles = get_poles_by_line_id(line_ids, line_id) - max_relief = np.maximum(poles[:, 2].max() - poles[:, 2].min(), max_relief) - - return (max_relief // z_cell_size + 2) * z_cell_size - - -def normalize_vertically(poles: np.ndarray, relief: float) -> np.ndarray: - """ - Given a set of pole locations, normalize the vertical component to the maximum relief across all lines. - - This ensures that the drape mesh has uniform vertical discretization across all survey lines. - - :param poles: Array of pole locations to normalize. - :param relief: Maximum relief across all survey lines, rounded to the nearest cell thickness. - - :return: Array of pole locations with normalized vertical component. - """ - min_poles_z = poles[:, 2].min() - poles[:, 2] -= min_poles_z - poles[:, 2] *= relief / np.maximum(poles[:, 2].max(), 1e-3) - - # Shift back vertically - poles[:, 2] += min_poles_z - - return poles - - -def get_poles_by_line_id(line_ids: IntegerData, uid: int) -> np.ndarray: - """Get the vertices associated with a given line ID.""" - mn_mask = line_ids.values == uid - - unique_tx = np.unique(line_ids.parent.ab_cell_id.values[mn_mask]) - - ab_mask = np.isin(line_ids.parent.complement.ab_cell_id.values, unique_tx) - - return np.vstack( - [ - line_ids.parent.vertices[line_ids.parent.cells[mn_mask].flatten()], - line_ids.parent.current_electrodes.vertices[ - line_ids.parent.current_electrodes.cells[ab_mask].flatten() - ], - ] - ) diff --git a/simpeg_drivers/utils/surveys.py b/simpeg_drivers/utils/surveys.py index d01c840a..13e55603 100644 --- a/simpeg_drivers/utils/surveys.py +++ b/simpeg_drivers/utils/surveys.py @@ -13,13 +13,19 @@ import numpy as np from discretize import TreeMesh -from geoapps_utils.utils.numerical import traveling_salesman from geoh5py import Workspace -from geoh5py.objects import PotentialElectrode +from geoh5py.data import IntegerData +from geoh5py.objects import DrapeModel, PotentialElectrode +from geoh5py.shared.merging.drape_model import DrapeModelMerger from scipy.sparse import csgraph, csr_matrix from scipy.spatial import cKDTree from simpeg.survey import BaseSurvey +from simpeg_drivers.options import ( + DrapeModelOptions, +) +from simpeg_drivers.utils.utils import get_drape_model + def station_spacing( locations: np.ndarray, @@ -63,47 +69,6 @@ def counter_clockwise_sort(segments: np.ndarray, vertices: np.ndarray) -> np.nda return segments -def compute_alongline_distance(points: np.ndarray, ordered: bool = True): - """ - Convert from cartesian (x, y, values) points to (distance, values) locations. - - :param: points: Cartesian coordinates of points lying either roughly within a - plane or a line. - """ - if not ordered: - order = traveling_salesman(points) - points = points[order, :] - - distances = np.cumsum( - np.r_[0, np.linalg.norm(np.diff(points[:, :2], axis=0), axis=1)] - ) - if points.shape[1] > 2: - distances = np.c_[distances, points[:, 2:]] - - return distances - - -def copy_potentials_from_mask( - workspace: Workspace, survey: PotentialElectrode, cell_mask: np.ndarray -): - """ - Returns a survey containing data from a single line. - - :param workspace: geoh5py workspace containing a valid DCIP survey. - :param survey: PotentialElectrode object. - :param cell_mask: Boolean array of M-N pairs to include in the new survey. - """ - - if not np.any(cell_mask): - raise ValueError("No cells found in the mask.") - - active_poles = np.zeros(survey.n_vertices, dtype=bool) - active_poles[survey.cells[cell_mask, :].ravel()] = True - potentials = survey.copy(parent=workspace, mask=active_poles, cell_mask=cell_mask) - - return potentials - - def get_intersecting_cells(locations: np.ndarray, mesh: TreeMesh) -> np.ndarray: """ Find cells that intersect with a set of segments. @@ -166,7 +131,9 @@ def get_parts_from_electrodes(survey: PotentialElectrode) -> np.ndarray: ) connections = csgraph.connected_components(edge_array)[1] - return connections[survey.cells[:, 0]] + parts = connections[survey.cells[:, 0]] + _, u_part = np.unique(parts, return_inverse=True) + return u_part def compute_em_projections(locations, simulation): @@ -206,3 +173,102 @@ def compute_dc_projections(locations, cells, simulation): proj_mn -= projection[cells[indices, 1], :] receiver.spatialP = proj_mn # pylint: disable=protected-access + + +def create_mesh_by_line_id( + workspace: Workspace, + line_ids: IntegerData, + drape_options: DrapeModelOptions, + **object_kwargs, +) -> DrapeModel: + """ + Create a drape mesh for the dc resistivity survey lines. + + :param workspace: Workspace to create the drape mesh in. + :param line_ids: IntegerData object containing the line IDs for each vertex. + :param drape_options: DrapeModelOptions containing the parameters for the drape mesh + :param object_kwargs: Additional keyword arguments to pass to the DrapeModelMerger.create_object method. + + :return: A DrapeModel object containing the merged drape mesh for all survey lines. + """ + drape_models = [] + temp_work = Workspace() + + relief = get_max_line_relief(line_ids, drape_options.v_cell_size) + + for line_id in np.unique(line_ids.values): + poles = get_poles_by_line_id(line_ids, line_id) + poles = np.unique(poles, axis=0) + poles = normalize_vertically(poles, relief) + + drape_model = get_drape_model( + temp_work, + poles, + [ + drape_options.u_cell_size, + drape_options.v_cell_size, + ], + drape_options.depth_core, + [drape_options.horizontal_padding] * 2 + + [drape_options.vertical_padding, 1], + drape_options.expansion_factor, + ) + drape_models.append(drape_model) + + entity = DrapeModelMerger.create_object(workspace, drape_models, **object_kwargs) + + return entity + + +def get_max_line_relief(line_ids: IntegerData, z_cell_size: float) -> float: + """ + Get the maximum relief across all survey lines, rounded to the nearest cell thickness. + + :param line_ids: IntegerData object containing the line IDs for each vertex. + :param z_cell_size: Cell size in the vertical direction for the drape mesh. + """ + max_relief = 0 + for line_id in np.unique(line_ids.values): + poles = get_poles_by_line_id(line_ids, line_id) + max_relief = np.maximum(poles[:, 2].max() - poles[:, 2].min(), max_relief) + + return (max_relief // z_cell_size + 2) * z_cell_size + + +def get_poles_by_line_id(line_ids: IntegerData, uid: int) -> np.ndarray: + """Get the vertices associated with a given line ID.""" + mn_mask = line_ids.values == uid + + unique_tx = np.unique(line_ids.parent.ab_cell_id.values[mn_mask]) + + ab_mask = np.isin(line_ids.parent.complement.ab_cell_id.values, unique_tx) + + return np.vstack( + [ + line_ids.parent.vertices[line_ids.parent.cells[mn_mask].flatten()], + line_ids.parent.current_electrodes.vertices[ + line_ids.parent.current_electrodes.cells[ab_mask].flatten() + ], + ] + ) + + +def normalize_vertically(poles: np.ndarray, relief: float) -> np.ndarray: + """ + Given a set of pole locations, normalize the vertical component to the maximum relief across all lines. + + This ensures that the drape mesh has uniform vertical discretization across all survey lines. + + :param poles: Array of pole locations to normalize. + :param relief: Maximum relief across all survey lines, rounded to the nearest cell thickness. + + :return: Array of pole locations with normalized vertical component. + """ + min_poles_z = poles[:, 2].min() + poles[:, 2] -= min_poles_z + poles[:, 2] *= relief / np.maximum(poles[:, 2].max(), 1e-3) + + # Shift back vertically + poles[:, 2] += min_poles_z + + return poles diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 407bc13c..3d7fc4dd 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -22,7 +22,7 @@ from geoapps_utils.utils.numerical import running_mean, traveling_salesman from geoh5py import Workspace from geoh5py.data import NumericData -from geoh5py.groups import Group, SimPEGGroup +from geoh5py.groups import SimPEGGroup from geoh5py.objects import DrapeModel, Octree from geoh5py.objects.surveys.direct_current import PotentialElectrode from geoh5py.objects.surveys.electromagnetics.airborne_app_con import ( @@ -32,14 +32,10 @@ from geoh5py.shared import INTEGER_NDV from geoh5py.ui_json import InputFile from grid_apps.utils import octree_2_treemesh -from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator, interp1d -from scipy.sparse import csr_matrix, diags -from scipy.spatial import ConvexHull, Delaunay, cKDTree +from scipy.interpolate import interp1d +from scipy.spatial import ConvexHull, cKDTree from simpeg_drivers import DRIVER_MAP -from simpeg_drivers.utils.surveys import ( - compute_alongline_distance, -) if TYPE_CHECKING: @@ -629,3 +625,23 @@ def simpeg_group_to_driver(group: SimPEGGroup, workspace: Workspace) -> Inversio params = inversion_driver._params_class.build(ifile) # pylint: disable=protected-access return inversion_driver(params) + + +def compute_alongline_distance(points: np.ndarray, ordered: bool = True): + """ + Convert from cartesian (x, y, values) points to (distance, values) locations. + + :param: points: Cartesian coordinates of points lying either roughly within a + plane or a line. + """ + if not ordered: + order = traveling_salesman(points) + points = points[order, :] + + distances = np.cumsum( + np.r_[0, np.linalg.norm(np.diff(points[:, :2], axis=0), axis=1)] + ) + if points.shape[1] > 2: + distances = np.c_[distances, points[:, 2:]] + + return distances From 3ff5a6a2942fd83b6a744039716ba66324f7dcd6 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Mon, 9 Mar 2026 16:44:35 -0700 Subject: [PATCH 10/24] Add unitest for new utils --- tests/utils_surveys_test.py | 53 +++++++++++++++++++++++++++++++++++-- 1 file changed, 51 insertions(+), 2 deletions(-) diff --git a/tests/utils_surveys_test.py b/tests/utils_surveys_test.py index 3546c940..13b31652 100644 --- a/tests/utils_surveys_test.py +++ b/tests/utils_surveys_test.py @@ -12,9 +12,15 @@ import numpy as np from geoh5py import Workspace -from geoh5py.objects import Points +from geoh5py.objects import CurrentElectrode, Points, PotentialElectrode -from simpeg_drivers.utils.surveys import counter_clockwise_sort, station_spacing +from simpeg_drivers.options import DrapeModelOptions +from simpeg_drivers.utils.surveys import ( + counter_clockwise_sort, + create_mesh_by_line_id, + get_parts_from_electrodes, + station_spacing, +) def create_test_survey( @@ -83,3 +89,46 @@ def test_counterclockwise_sort(): ccw_sorted = counter_clockwise_sort(segments, vertices) np.testing.assert_equal(ccw_sorted[0, :], [0, 5]) + + +def get_dc_survey(workspace): + """ + Create a DC survey with 4 current electrodes and 10 potential electrodes. + """ + vertices = np.random.randn(4, 3) + currents = CurrentElectrode.create(workspace, vertices=vertices, parts=[0, 0, 1, 1]) + currents.ab_cell_id = np.repeat([1, 2], 2) + + rx_vertices = np.random.randn(12, 3) + mn_pairs = np.c_[np.arange(11), np.arange(1, 12)] # Remove connection + mn_pairs = np.delete(mn_pairs, 5, axis=0) + potentials = PotentialElectrode.create( + workspace, + vertices=rx_vertices, + cells=mn_pairs, + ) + potentials.ab_cell_id = np.repeat([1, 2], 5) + + potentials.current_electrodes = currents + return potentials + + +def test_parts_from_electrodes(): + workspace = Workspace() + survey = get_dc_survey(workspace) + + line_ids = get_parts_from_electrodes(survey) + assert len(np.unique(line_ids)) == 2 + assert np.all(line_ids[:5] == 0) + assert np.all(line_ids[5:] == 1) + + +def test_drape_from_line_id(tmp_path): + + with Workspace.create(tmp_path / f"{__name__}.geoh5") as ws: + survey = get_dc_survey(ws) + drape = create_mesh_by_line_id( + ws, survey.ab_cell_id, DrapeModelOptions(), name="test_drape" + ) + + assert drape.name == "test_drape" From 61ea2d45e1944ef3c5736134afeb24c07da18325 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Mar 2026 23:48:53 +0000 Subject: [PATCH 11/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/uijson_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/uijson_test.py b/tests/uijson_test.py index db58145a..3931b837 100644 --- a/tests/uijson_test.py +++ b/tests/uijson_test.py @@ -21,10 +21,10 @@ from geoh5py.ui_json.annotations import Deprecated from packaging.version import Version from pydantic import AliasChoices, Field +from simpeg_drivers.line_sweep.driver import LineSweepDriver import simpeg_drivers from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.line_sweep.driver import LineSweepDriver from simpeg_drivers.options import Deprecations, IRLSOptions from simpeg_drivers.potential_fields.gravity.options import GravityInversionOptions from simpeg_drivers.potential_fields.gravity.uijson import GravityInversionUIJson From 77995ec89ae960a0cdd703a5a86bf8046109351a Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 09:44:31 -0700 Subject: [PATCH 12/24] Add Deprecations for old batch 2d --- simpeg_drivers/__init__.py | 4 +- simpeg_drivers/driver.py | 5 +- simpeg_drivers/electricals/base_2d.py | 53 +++++++++++++++++++ .../direct_current/two_dimensions/driver.py | 14 ++++- .../direct_current/two_dimensions/options.py | 18 ++----- .../two_dimensions/driver.py | 14 ++++- .../two_dimensions/options.py | 28 +++------- tests/uijson_test.py | 21 +++++--- 8 files changed, 105 insertions(+), 52 deletions(-) diff --git a/simpeg_drivers/__init__.py b/simpeg_drivers/__init__.py index f2c8c63b..2bc66099 100644 --- a/simpeg_drivers/__init__.py +++ b/simpeg_drivers/__init__.py @@ -69,7 +69,7 @@ def assets_path() -> Path: }, ), "direct current pseudo 3d": ( - "simpeg_drivers.electricals.direct_current.pseudo_three_dimensions.driver", + "simpeg_drivers.electricals.direct_current.two_dimensions.driver", { "forward": "DCBatch2DForwardDriver", "inversion": "DCBatch2DInversionDriver", @@ -115,7 +115,7 @@ def assets_path() -> Path: }, ), "induced polarization pseudo 3d": ( - "simpeg_drivers.electricals.induced_polarization.pseudo_three_dimensions.driver", + "simpeg_drivers.electricals.induced_polarization.two_dimensions.driver", { "forward": "IPBatch2DForwardDriver", "inversion": "IPBatch2DInversionDriver", diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index cf333289..728ac298 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -506,19 +506,18 @@ def params(self) -> BaseForwardOptions | BaseInversionOptions: @params.setter def params( self, - val: BaseForwardOptions | BaseInversionOptions | SweepParams, + val: BaseForwardOptions | BaseInversionOptions, ): if not isinstance( val, ( BaseForwardOptions, BaseInversionOptions, - SweepParams, BaseJointOptions, ), ): raise TypeError( - "Parameters must be of type 'BaseInversionOptions', 'BaseForwardOptions' or 'SweepParams'." + "Parameters must be of type 'BaseInversionOptions', 'BaseForwardOptions' or 'BaseJointOptions'." ) self._params = val diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 1126fced..e7dee40d 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -11,13 +11,25 @@ from __future__ import annotations +from logging import getLogger + +from geoh5py.objects import DrapeModel, Octree, PotentialElectrode from geoh5py.ui_json.ui_json import fetch_active_workspace +from pydantic import field_validator, model_validator from simpeg_drivers.components.meshes import InversionMesh from simpeg_drivers.driver import InversionDriver +from simpeg_drivers.options import ( + CoreOptions, + DrapeModelOptions, + LineSelectionOptions, +) from simpeg_drivers.utils.surveys import create_mesh_by_line_id +logger = getLogger(__name__) + + class Base2DDriver(InversionDriver): """ Base class for 2D DC and IP forward and inversion drivers. @@ -45,3 +57,44 @@ def inversion_mesh(self) -> InversionMesh: ) return self._inversion_mesh + + +class DeprecatedBatch2DDriver(Base2DDriver): + """Direct Current 2D forward driver.""" + + def __init__(self, *args, **kwargs): + logger.warning( + "The Batch2D classes will be deprecated in version 0.5.0. " + "Please use the non-batch classes instead. Results may be affected.", + ) + + super().__init__(*args, **kwargs) + + +class Base2DOptions(CoreOptions): + """ + Base options for the Direct Current 2D forward and inverse driver. + + + :param data_object: Potential electrode object. + :param line_selection: Line selection parameters. + :param mesh: Optional mesh object if providing a heterogeneous model. + :param drape_model: Drape model parameters. + """ + + data_object: PotentialElectrode + line_selection: LineSelectionOptions = LineSelectionOptions() + mesh: DrapeModel | Octree | None = None + drape_model: DrapeModelOptions = DrapeModelOptions() + + @field_validator("mesh", mode="before") + @classmethod + def mesh_cannot_be_octree(cls, value: Octree | DrapeModel): + if isinstance(value, Octree): + logger.warning( + "DC 2D forward and inversion no longer support Octree meshes as input. " + "The program will attempt to transfer models onto the newly created DrapeModel. " + "Results may be affected." + ) + return None + return value diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py index 5642d81a..fc26fca6 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py @@ -11,11 +11,23 @@ from __future__ import annotations -from simpeg_drivers.electricals.base_2d import Base2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver, DeprecatedBatch2DDriver from .options import DC2DForwardOptions, DC2DInversionOptions +class DCBatch2DForwardDriver(DeprecatedBatch2DDriver): + """Deprecated - Direct Current Batch Direct Current 2D forward driver.""" + + _params_class = DC2DForwardOptions + + +class DCBatch2DInversionDriver(DeprecatedBatch2DDriver): + """Deprecated - Direct Current Batch 2D inversion driver.""" + + _params_class = DC2DInversionOptions + + class DC2DForwardDriver(Base2DDriver): """Direct Current 2D forward driver.""" diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py index 468efad5..d60b882a 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/options.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/options.py @@ -15,25 +15,21 @@ from typing import ClassVar from geoh5py.data import FloatData -from geoh5py.objects import DrapeModel, PotentialElectrode from simpeg_drivers import assets_path +from simpeg_drivers.electricals.base_2d import Base2DOptions from simpeg_drivers.options import ( BaseForwardOptions, BaseInversionOptions, ConductivityModelOptions, - DrapeModelOptions, - LineSelectionOptions, ) -class DC2DForwardOptions(BaseForwardOptions): +class DC2DForwardOptions(BaseForwardOptions, Base2DOptions): """ Direct Current 2D forward options. :param potential_channel_bool: Potential channel boolean. - :param line_selection: Line selection parameters. - :param drape_model: Drape model parameters. """ name: ClassVar[str] = "Direct Current 2D Forward" @@ -45,15 +41,11 @@ class DC2DForwardOptions(BaseForwardOptions): physical_property: str = "conductivity" inversion_type: str = "direct current 2d" - data_object: PotentialElectrode potential_channel_bool: bool = True - line_selection: LineSelectionOptions = LineSelectionOptions() - mesh: DrapeModel | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions -class DC2DInversionOptions(BaseInversionOptions): +class DC2DInversionOptions(BaseInversionOptions, Base2DOptions): """ Direct Current 2D inversion options. @@ -72,10 +64,6 @@ class DC2DInversionOptions(BaseInversionOptions): physical_property: str = "conductivity" inversion_type: str = "direct current 2d" - data_object: PotentialElectrode potential_channel: FloatData potential_uncertainty: float | FloatData | None = None - line_selection: LineSelectionOptions = LineSelectionOptions() - mesh: DrapeModel | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() models: ConductivityModelOptions diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py index b1322e92..96139fc6 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py @@ -11,7 +11,7 @@ from __future__ import annotations -from simpeg_drivers.electricals.base_2d import Base2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver, DeprecatedBatch2DDriver from .options import ( IP2DForwardOptions, @@ -19,6 +19,18 @@ ) +class IPBatch2DForwardDriver(DeprecatedBatch2DDriver): + """Deprecated - Direct Current Batch Direct Current 2D forward driver.""" + + _params_class = IP2DForwardOptions + + +class IPBatch2DInversionDriver(DeprecatedBatch2DDriver): + """Deprecated - Direct Current Batch 2D inversion driver.""" + + _params_class = IP2DInversionOptions + + class IP2DForwardDriver(Base2DDriver): """Induced Polarization 2D forward driver.""" diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py index 29a33611..e3b9f5dd 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py @@ -15,26 +15,19 @@ from typing import ClassVar from geoh5py.data import FloatData -from geoh5py.objects import DrapeModel, PotentialElectrode from simpeg_drivers import assets_path +from simpeg_drivers.electricals.base_2d import Base2DOptions from simpeg_drivers.electricals.options import IPModelOptions -from simpeg_drivers.options import ( - BaseForwardOptions, - BaseInversionOptions, - DrapeModelOptions, - LineSelectionOptions, -) +from simpeg_drivers.options import BaseForwardOptions, BaseInversionOptions -class IP2DForwardOptions(BaseForwardOptions): +class IP2DForwardOptions(BaseForwardOptions, Base2DOptions): """ Induced Polarization 2D forward options. :param chargeability_channel_bool: Chargeability channel boolean. - :param mesh: Optional mesh object if providing a heterogeneous model. - :param drape_model: Drape model parameters. - :param line_selection: Line selection parameters. + :param models: Set of models options for IP """ name: ClassVar[str] = "Induced Polarization 2D Forward" @@ -46,22 +39,17 @@ class IP2DForwardOptions(BaseForwardOptions): physical_property: str = "chargeability" inversion_type: str = "induced polarization 2d" - data_object: PotentialElectrode chargeability_channel_bool: bool = True - line_selection: LineSelectionOptions = LineSelectionOptions() - mesh: DrapeModel | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() models: IPModelOptions -class IP2DInversionOptions(BaseInversionOptions): +class IP2DInversionOptions(BaseInversionOptions, Base2DOptions): """ Induced Polarization 2D inversion options. :param chargeability_channel: Chargeability data channel. :param chargeability_uncertainty: Chargeability data uncertainty channel. - :param line_selection: Line selection parameters. - :param drape_model: Drape model parameters. + :param models: Set of models options for IP """ name: ClassVar[str] = "Induced Polarization 2D Inversion" @@ -73,10 +61,6 @@ class IP2DInversionOptions(BaseInversionOptions): physical_property: str = "chargeability" inversion_type: str = "induced polarization 2d" - data_object: PotentialElectrode chargeability_channel: FloatData chargeability_uncertainty: float | FloatData | None = None - line_selection: LineSelectionOptions = LineSelectionOptions() - mesh: DrapeModel | None = None - drape_model: DrapeModelOptions = DrapeModelOptions() models: IPModelOptions diff --git a/tests/uijson_test.py b/tests/uijson_test.py index db58145a..95f1f1b2 100644 --- a/tests/uijson_test.py +++ b/tests/uijson_test.py @@ -24,13 +24,13 @@ import simpeg_drivers from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.line_sweep.driver import LineSweepDriver from simpeg_drivers.options import Deprecations, IRLSOptions from simpeg_drivers.potential_fields.gravity.options import GravityInversionOptions from simpeg_drivers.potential_fields.gravity.uijson import GravityInversionUIJson from simpeg_drivers.uijson import SimPEGDriversUIJson from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents from simpeg_drivers.utils.synthetics.options import ( + DrapeModelOptions, MeshOptions, ModelOptions, SurveyOptions, @@ -307,7 +307,7 @@ def test_gravity_uijson(tmp_path): } -def test_legacy_uijson(tmp_path: Path): +def test_legacy_uijson(tmp_path: Path, caplog): """ Loop over all uijson files in the legacy directory and check that the read and run still works. @@ -326,6 +326,9 @@ def test_legacy_uijson(tmp_path: Path): if inversion_type not in CHANNEL_NAME: continue + if "pseudo" in inversion_type: + pass + forward = ifile.data.get("forward_only", None) work_path = version_path / ( @@ -339,7 +342,7 @@ def test_legacy_uijson(tmp_path: Path): n_stations=10, n_lines=3, ), - mesh=MeshOptions(), + mesh=DrapeModelOptions() if "2d" in inversion_type else MeshOptions(), model=ModelOptions( background=1.0, anomaly=2.0, @@ -392,10 +395,12 @@ def test_legacy_uijson(tmp_path: Path): driver = InversionDriver.from_input_file(ifile.data) - if hasattr(driver.params, "cooling_factor"): - assert driver.params.cooling_factor == 4.0 + with caplog.at_level(logging.WARNING): + params = driver._params_class.build(ifile) # pylint: disable=protected-access + driver = driver(params) - if isinstance(driver, LineSweepDriver): - continue + if "pseudo" in inversion_type: + assert "no longer support Octree meshes" in caplog.text + assert "The Batch2D classes will be deprecated" in caplog.text - assert driver.inversion + assert driver.models From a032fba405f7140a480aa74c0f77b6f86e2283b6 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 11:05:34 -0700 Subject: [PATCH 13/24] Rework logic for dip dir orientations --- simpeg_drivers/components/models.py | 6 ++-- simpeg_drivers/options.py | 23 ++++-------- simpeg_drivers/utils/regularization.py | 48 +++++++++++--------------- tests/utils_regularization_test.py | 24 ++++++++++--- 4 files changed, 49 insertions(+), 52 deletions(-) diff --git a/simpeg_drivers/components/models.py b/simpeg_drivers/components/models.py index d3b8a9fb..8e0a7771 100644 --- a/simpeg_drivers/components/models.py +++ b/simpeg_drivers/components/models.py @@ -458,14 +458,14 @@ def gradient_dip(self) -> np.ndarray | None: if self._gradient_dip.model is None: return None - return self._gradient_dip.model.copy() + return np.deg2rad(self._gradient_dip.model) @property def gradient_direction(self) -> np.ndarray | None: if self._gradient_direction.model is None: return None - return self._gradient_direction.model.copy() + return np.deg2rad(self._gradient_direction.model) def remove_air(self, active_cells: np.ndarray): """Use active cells vector to remove air cells from model""" @@ -643,8 +643,6 @@ def _get_value(self, model: float | NumericData) -> np.ndarray: elif isinstance(model, int | float): nc = self.driver.inversion_mesh.mesh.n_cells model *= np.ones(nc) - elif isinstance(model, np.ndarray): - model = (self.driver.inversion_mesh.permutation @ model).astype(model.dtype) return model diff --git a/simpeg_drivers/options.py b/simpeg_drivers/options.py index 4d9035e0..556c11a2 100644 --- a/simpeg_drivers/options.py +++ b/simpeg_drivers/options.py @@ -18,7 +18,7 @@ import numpy as np from geoapps_utils.base import Options -from geoapps_utils.utils.numerical import weighted_average +from geoh5py import Workspace from geoh5py.data import ( BooleanData, DataAssociationEnum, @@ -282,19 +282,19 @@ class ModelOptions(BaseModel): y_norm: float | FloatData | None = 2.0 z_norm: float | FloatData = 2.0 - _gradient_orientations: np.ndarray | None = None + _gradient_orientations: list[FloatData] | None = None @property - def gradient_direction(self) -> np.ndarray: + def gradient_direction(self) -> FloatData | None: if self.gradient_orientations is None: return None - return self.gradient_orientations[:, 0] + return self.gradient_orientations[0] @property - def gradient_dip(self) -> np.ndarray: + def gradient_dip(self) -> FloatData | None: if self.gradient_orientations is None: return None - return self.gradient_orientations[:, 1] + return self.gradient_orientations[1] @property def gradient_orientations(self) -> tuple(float, float): @@ -307,16 +307,7 @@ def gradient_orientations(self) -> tuple(float, float): if self._gradient_orientations is None and self.gradient_rotation is not None: orientations = direction_and_dip(self.gradient_rotation) - - angles = np.deg2rad(orientations) - # Deal with aircells here - orientations = weighted_average( - self.gradient_rotation.parent.centroids, - self.gradient_rotation.parent.centroids, - [angles[:, 0], angles[:, 1]], - ) - - self._gradient_orientations = np.vstack(orientations).T + self._gradient_orientations = orientations return self._gradient_orientations diff --git a/simpeg_drivers/utils/regularization.py b/simpeg_drivers/utils/regularization.py index 96802c6b..27b97dee 100644 --- a/simpeg_drivers/utils/regularization.py +++ b/simpeg_drivers/utils/regularization.py @@ -16,6 +16,7 @@ x_rotation_matrix, z_rotation_matrix, ) +from geoh5py.data import Data from geoh5py.groups import PropertyGroup from geoh5py.groups.property_group_type import GroupTypeEnum from simpeg.regularization import SparseSmoothness @@ -436,35 +437,15 @@ def rotated_gradient( return unit_grad -def ensure_dip_direction_convention( - orientations: np.ndarray, group_type: str -) -> np.ndarray: - """ - Ensure orientations array has dip and direction convention. - - :param orientations: Array of orientations. Either n * 2 if Strike & dip - or Dip direction & dip group_type, or n * 3 if 3D Vector group_type defining the normal of the dipping plane. - :param group_type as specified in geoh5py.GroupTypeEnum. - """ - - if group_type == GroupTypeEnum.VECTOR: - orientations = np.rad2deg(cartesian_normal_to_direction_and_dip(orientations)) - - if group_type in [GroupTypeEnum.STRIKEDIP]: - orientations[:, 0] = 90.0 + orientations[:, 0] - - return orientations - - -def direction_and_dip(property_group: PropertyGroup) -> list[np.ndarray]: +def direction_and_dip(property_group: PropertyGroup) -> list[Data]: """Conversion of orientation group to direction and dip.""" group_type = property_group.property_group_type - if group_type not in [ - GroupTypeEnum.VECTOR, - GroupTypeEnum.STRIKEDIP, - GroupTypeEnum.DIPDIR, - ]: + + if group_type == GroupTypeEnum.DIPDIR: + return [property_group.parent.get_data(k)[0] for k in property_group.properties] + + if group_type not in [GroupTypeEnum.VECTOR, GroupTypeEnum.STRIKEDIP]: raise ValueError( "Property group does not contain orientation data. " "Type must be one of '3D vector', 'Strike & dip', or " @@ -475,7 +456,20 @@ def direction_and_dip(property_group: PropertyGroup) -> list[np.ndarray]: [property_group.parent.get_data(k)[0].values for k in property_group.properties] ).T - return ensure_dip_direction_convention(orientations, group_type) + if group_type == GroupTypeEnum.VECTOR: + orientations = np.rad2deg(cartesian_normal_to_direction_and_dip(orientations)) + + if group_type in [GroupTypeEnum.STRIKEDIP]: + orientations[:, 0] = 90.0 + orientations[:, 0] + + dir_dip = property_group.parent.add_data( + { + "azimuth": {"values": orientations[:, 0]}, + "dip": {"values": orientations[:, 1]}, + } + ) + + return dir_dip def set_rotated_operators( diff --git a/tests/utils_regularization_test.py b/tests/utils_regularization_test.py index a4bdffc9..9cd56854 100644 --- a/tests/utils_regularization_test.py +++ b/tests/utils_regularization_test.py @@ -10,12 +10,14 @@ import numpy as np from discretize import TreeMesh +from geoh5py import Workspace +from geoh5py.objects import Points from simpeg_drivers.utils.regularization import ( cell_neighbors, cell_neighbors_along_axis, cell_neighbors_lists, - ensure_dip_direction_convention, + direction_and_dip, ) @@ -91,7 +93,8 @@ def test_collect_all_neighbors(): assert [15, 15, 15] not in neighbor_centers -def test_ensure_dip_direction_convention(): +def test_ensure_dip_direction_convention(tmp_path): + # Rotate the vertical unit vector 37 degrees to the west and then 16 # degrees counter-clockwise about the z-axis. Should result in a dip # direction of 270 - 16 = 254 and a dip of 37. @@ -128,6 +131,17 @@ def test_ensure_dip_direction_convention(): arbitrary_vector.tolist(), ] ) - dir_dip = ensure_dip_direction_convention(orientations, group_type="3D vector") - assert np.allclose(dir_dip[:, 0], [90, 0, 270, 180] * 2 + [254]) - assert np.allclose(dir_dip[:, 1], [45] * 4 + [30] * 4 + [37]) + with Workspace.create(tmp_path / f"{__name__}.geoh5") as ws: + pts = Points.create(ws, vertices=np.random.randn(orientations.shape[0], 3)) + data_list = pts.add_data( + { + f"{ax}_vec": {"values": ori} + for ax, ori in zip("xyz", orientations.T, strict=True) + } + ) + prop_group = pts.create_property_group( + name="vector", property_group_type="3D vector", properties=data_list + ) + dir_dip = direction_and_dip(prop_group) + assert np.allclose(dir_dip[0].values, [90, 0, 270, 180] * 2 + [254]) + assert np.allclose(dir_dip[1].values, [45] * 4 + [30] * 4 + [37]) From bbf95e556c90b92f6e354a31ddb86b4f815fcc9f Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 11:47:42 -0700 Subject: [PATCH 14/24] Simplify logic --- simpeg_drivers/__init__.py | 8 +++---- .../factories/simulation_factory.py | 7 ++---- simpeg_drivers/electricals/base_2d.py | 24 +++++++++---------- .../direct_current/two_dimensions/driver.py | 14 +---------- .../two_dimensions/driver.py | 14 +---------- 5 files changed, 20 insertions(+), 47 deletions(-) diff --git a/simpeg_drivers/__init__.py b/simpeg_drivers/__init__.py index 2bc66099..3cc0ed27 100644 --- a/simpeg_drivers/__init__.py +++ b/simpeg_drivers/__init__.py @@ -71,8 +71,8 @@ def assets_path() -> Path: "direct current pseudo 3d": ( "simpeg_drivers.electricals.direct_current.two_dimensions.driver", { - "forward": "DCBatch2DForwardDriver", - "inversion": "DCBatch2DInversionDriver", + "forward": "DC2DForwardDriver", + "inversion": "DC2DInversionDriver", }, ), "fdem": ( @@ -117,8 +117,8 @@ def assets_path() -> Path: "induced polarization pseudo 3d": ( "simpeg_drivers.electricals.induced_polarization.two_dimensions.driver", { - "forward": "IPBatch2DForwardDriver", - "inversion": "IPBatch2DInversionDriver", + "forward": "IP2DForwardDriver", + "inversion": "IP2DInversionDriver", }, ), "joint cross gradient": ( diff --git a/simpeg_drivers/components/factories/simulation_factory.py b/simpeg_drivers/components/factories/simulation_factory.py index ec45a293..8c47269f 100644 --- a/simpeg_drivers/components/factories/simulation_factory.py +++ b/simpeg_drivers/components/factories/simulation_factory.py @@ -67,7 +67,7 @@ def concrete_object(self): return simulation.Simulation3DIntegral - if self.factory_type in ["direct current 3d", "direct current pseudo 3d"]: + if self.factory_type == "direct current 3d": from simpeg.electromagnetics.static.resistivity import simulation return simulation.Simulation3DNodal @@ -77,10 +77,7 @@ def concrete_object(self): return simulation_2d.Simulation2DNodal - if self.factory_type in [ - "induced polarization 3d", - "induced polarization pseudo 3d", - ]: + if self.factory_type == "induced polarization 3d": from simpeg.electromagnetics.static.induced_polarization import simulation return simulation.Simulation3DNodal diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index e7dee40d..1deccd2d 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -51,6 +51,7 @@ def inversion_mesh(self) -> InversionMesh: self.params.drape_model, parent=self.out_group, ) + self.params.mesh = entity self._inversion_mesh = InversionMesh( self.workspace, self.params, entity=entity @@ -59,18 +60,6 @@ def inversion_mesh(self) -> InversionMesh: return self._inversion_mesh -class DeprecatedBatch2DDriver(Base2DDriver): - """Direct Current 2D forward driver.""" - - def __init__(self, *args, **kwargs): - logger.warning( - "The Batch2D classes will be deprecated in version 0.5.0. " - "Please use the non-batch classes instead. Results may be affected.", - ) - - super().__init__(*args, **kwargs) - - class Base2DOptions(CoreOptions): """ Base options for the Direct Current 2D forward and inverse driver. @@ -98,3 +87,14 @@ def mesh_cannot_be_octree(cls, value: Octree | DrapeModel): ) return None return value + + @field_validator("inversion_type", mode="before") + @classmethod + def deprecated_pseudo(cls, value: str): + if "pseudo 3d" in value: + logger.warning( + "The Batch2D classes will be deprecated in version 0.5.0. " + "Please use the non-batch classes instead. Results may be affected.", + ) + value = value.replace("pseudo 3d", "2d") + return value diff --git a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py index fc26fca6..5642d81a 100644 --- a/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/direct_current/two_dimensions/driver.py @@ -11,23 +11,11 @@ from __future__ import annotations -from simpeg_drivers.electricals.base_2d import Base2DDriver, DeprecatedBatch2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver from .options import DC2DForwardOptions, DC2DInversionOptions -class DCBatch2DForwardDriver(DeprecatedBatch2DDriver): - """Deprecated - Direct Current Batch Direct Current 2D forward driver.""" - - _params_class = DC2DForwardOptions - - -class DCBatch2DInversionDriver(DeprecatedBatch2DDriver): - """Deprecated - Direct Current Batch 2D inversion driver.""" - - _params_class = DC2DInversionOptions - - class DC2DForwardDriver(Base2DDriver): """Direct Current 2D forward driver.""" diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py index 96139fc6..b1322e92 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/driver.py @@ -11,7 +11,7 @@ from __future__ import annotations -from simpeg_drivers.electricals.base_2d import Base2DDriver, DeprecatedBatch2DDriver +from simpeg_drivers.electricals.base_2d import Base2DDriver from .options import ( IP2DForwardOptions, @@ -19,18 +19,6 @@ ) -class IPBatch2DForwardDriver(DeprecatedBatch2DDriver): - """Deprecated - Direct Current Batch Direct Current 2D forward driver.""" - - _params_class = IP2DForwardOptions - - -class IPBatch2DInversionDriver(DeprecatedBatch2DDriver): - """Deprecated - Direct Current Batch 2D inversion driver.""" - - _params_class = IP2DInversionOptions - - class IP2DForwardDriver(Base2DDriver): """Induced Polarization 2D forward driver.""" From 5dbb4b918c87a5aa058237263b02c67ba8f3d68d Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 12:55:32 -0700 Subject: [PATCH 15/24] Fix issue with non-monotonic line ids --- simpeg_drivers/components/topography.py | 4 +++- simpeg_drivers/electricals/base_2d.py | 12 +++++++----- simpeg_drivers/utils/nested.py | 17 +++++++++++++---- 3 files changed, 23 insertions(+), 10 deletions(-) diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index 40f20a39..4e8e8125 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -91,8 +91,10 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: "magnetotellurics", "direct current 3d", "direct current 2d", + "direct current pseudo 3d", "induced polarization 3d", "induced polarization 2d", + "induced polarization pseudo 3d", "apparent conductivity", ] or isinstance(data.entity, LargeLoopGroundEMSurvey) @@ -105,7 +107,7 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: ) else: - if any(k in self.params.inversion_type for k in ["2d", "p3d"]): + if any(k in self.params.inversion_type for k in ["2d", "pseudo"]): vertices = self.locations cells = getattr( self.params.active_cells.topography_object, "cells", None diff --git a/simpeg_drivers/electricals/base_2d.py b/simpeg_drivers/electricals/base_2d.py index 1deccd2d..fb8fc7f6 100644 --- a/simpeg_drivers/electricals/base_2d.py +++ b/simpeg_drivers/electricals/base_2d.py @@ -88,13 +88,15 @@ def mesh_cannot_be_octree(cls, value: Octree | DrapeModel): return None return value - @field_validator("inversion_type", mode="before") + @model_validator(mode="before") @classmethod - def deprecated_pseudo(cls, value: str): - if "pseudo 3d" in value: + def deprecated_pseudo(cls, data: dict): + if "pseudo 3d" in data.get("inversion_type", ""): logger.warning( "The Batch2D classes will be deprecated in version 0.5.0. " "Please use the non-batch classes instead. Results may be affected.", ) - value = value.replace("pseudo 3d", "2d") - return value + data["inversion_type"] = data["inversion_type"].replace("pseudo 3d", "2d") + data["line_selection"]["line_id"] = None + + return data diff --git a/simpeg_drivers/utils/nested.py b/simpeg_drivers/utils/nested.py index 4e4823e3..94c3a20d 100644 --- a/simpeg_drivers/utils/nested.py +++ b/simpeg_drivers/utils/nested.py @@ -304,12 +304,21 @@ def create_simulation( ), ) local_actives = mapping.local_active - # For DCIP-2D + # For DCIP-2D, create a projection from the global active cells to + # the local active cells else: - # Create a projection from the global active cells to the local active cells - active_mesh_part = np.isin( - simulation.mesh.parts, np.unique(local_survey.line_ids) + # Map the line_ids to the mesh parts (assumes sequential numbering) + line_number = ( + np.where( + np.isin( + np.unique(simulation.survey.line_ids), + np.unique(local_survey.line_ids), + ) + )[0] + + 1 ) + + active_mesh_part = np.isin(simulation.mesh.parts, line_number) n_actives = simulation.active_cells.sum() activate_ind = np.zeros(simulation.mesh.n_cells, dtype=int) activate_ind[np.where(simulation.active_cells)[0]] = np.arange(n_actives) From 51de9966285ad8e46d09aefabccf7dce68b91617 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 12:56:44 -0700 Subject: [PATCH 16/24] Remove old import --- tests/uijson_test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/uijson_test.py b/tests/uijson_test.py index 78f6a32c..95f1f1b2 100644 --- a/tests/uijson_test.py +++ b/tests/uijson_test.py @@ -21,7 +21,6 @@ from geoh5py.ui_json.annotations import Deprecated from packaging.version import Version from pydantic import AliasChoices, Field -from simpeg_drivers.line_sweep.driver import LineSweepDriver import simpeg_drivers from simpeg_drivers.driver import InversionDriver From fc548729e6db5c017ad8cea13bcfe80e1f9f4a63 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 14:27:00 -0700 Subject: [PATCH 17/24] Fix loose handling of synthetic entities --- simpeg_drivers/utils/synthetics/driver.py | 44 +++++++------------ .../driver_joint_cross_gradient_test.py | 37 +++++++++------- .../driver_joint_pgi_homogeneous_test.py | 24 ++++------ 3 files changed, 43 insertions(+), 62 deletions(-) diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index ec921484..50fc87cf 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -43,67 +43,53 @@ def __init__( self._model: FloatData | None = None @property - def topography(self): - entity = self.geoh5.get_entity("topography")[0] - assert isinstance(entity, Surface | type(None)) - if entity is None: - assert self.options is not None + def topography(self) -> Surface: + if self._topography is None: entity = get_topography_surface( geoh5=self.geoh5, options=self.options, ) - self._topography = entity + self._topography = entity return self._topography @property - def survey(self): - entity = self.geoh5.get_entity(self.options.survey.name)[0] - assert isinstance(entity, ObjectBase | type(None)) - if entity is None: - assert self.options is not None + def survey(self) -> ObjectBase: + if self._survey is None: entity = get_survey( geoh5=self.geoh5, method=self.options.method, options=self.options.survey, ) - self._survey = entity + self._survey = entity return self._survey @property - def mesh(self): - entity = self.geoh5.get_entity("mesh")[0] - assert isinstance(entity, Octree | DrapeModel | type(None)) - if entity is None: - assert self.options is not None + def mesh(self) -> Octree | DrapeModel: + if self._mesh is None: entity = get_mesh( self.options.method, survey=self.survey, topography=self.topography, options=self.options.mesh, ) - self._mesh = entity + self._mesh = entity return self._mesh @property - def active(self): - entity = self.mesh.get_entity(self.options.active.name)[0] - assert isinstance(entity, BooleanData | type(None)) - if entity is None: + def active(self) -> FloatData: + if self._active is None: entity = get_active(self.mesh, self.topography) - self._active = entity + self._active = entity return self._active @property - def model(self): - entity = self.mesh.get_entity(self.options.model.name)[0] - assert isinstance(entity, FloatData | type(None)) - if entity is None: - assert self.options is not None + def model(self) -> FloatData: + if self._model is None: entity = get_model( method=self.options.method, mesh=self.mesh, active=self.active.values, options=self.options.model, ) - self._model = entity + self._model = entity return self._model diff --git a/tests/run_tests/driver_joint_cross_gradient_test.py b/tests/run_tests/driver_joint_cross_gradient_test.py index 6f339c4a..d306a53d 100644 --- a/tests/run_tests/driver_joint_cross_gradient_test.py +++ b/tests/run_tests/driver_joint_cross_gradient_test.py @@ -11,8 +11,7 @@ from pathlib import Path import numpy as np -from geoh5py.data import FloatData -from geoh5py.groups import SimPEGGroup +from geoh5py.objects import CurrentElectrode, Octree, Points from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.three_dimensions import ( @@ -170,25 +169,29 @@ def test_joint_cross_gradient_inv_run( topography = geoh5.get_entity("topography")[0] drivers = [] orig_data = [] - - for suffix in "ABC": - components = SyntheticsComponents( - geoh5=geoh5, - options=SyntheticsComponentsOptions( - method="joint", - survey=SurveyOptions(name=f"survey {suffix}"), - mesh=MeshOptions(name=f"mesh {suffix}"), - model=ModelOptions(name=f"model {suffix}"), - active=SyntheticsActiveCellsOptions(name=f"active {suffix}"), - ), + origin = None + for name in [ + "Gravity Forward", + "Magnetic Vector Forward", + "Direct Current 3D Forward", + ]: + group = geoh5.get_entity(name)[0] + mesh = next(child for child in group.children if isinstance(child, Octree)) + survey = next( + child + for child in group.children + if isinstance(child, Points) and not isinstance(child, CurrentElectrode) ) - mesh = components.mesh - survey = components.survey + if origin is None: + origin = mesh.origin + else: + mesh.origin = origin + data = next(k for k in survey.children if "Iteration_0" in k.name) orig_data.append(data.values) - if suffix == "A": + if name == "Gravity Forward": params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, @@ -205,7 +208,7 @@ def test_joint_cross_gradient_inv_run( chi_factor=0.8, ) drivers.append(GravityInversionDriver(params)) - elif suffix == "C": + elif name == "Direct Current 3D Forward": params = DC3DInversionOptions.build( geoh5=geoh5, mesh=mesh, diff --git a/tests/run_tests/driver_joint_pgi_homogeneous_test.py b/tests/run_tests/driver_joint_pgi_homogeneous_test.py index 8c1bdfe4..111eddb5 100644 --- a/tests/run_tests/driver_joint_pgi_homogeneous_test.py +++ b/tests/run_tests/driver_joint_pgi_homogeneous_test.py @@ -13,9 +13,8 @@ from pathlib import Path import numpy as np -from geoh5py.data import FloatData from geoh5py.groups.property_group import GroupTypeEnum, PropertyGroup -from geoh5py.groups.simpeg import SimPEGGroup +from geoh5py.objects import Octree, Points from geoh5py.workspace import Workspace from simpeg_drivers.joint.joint_petrophysics.driver import JointPetrophysicsDriver @@ -144,21 +143,14 @@ def test_homogeneous_run( petrophysics = None gradient_rotation = None - for suffix in "AB": - components = SyntheticsComponents( - geoh5=geoh5, - options=SyntheticsComponentsOptions( - method="joint", - survey=SurveyOptions(name=f"survey {suffix}"), - mesh=MeshOptions(name=f"mesh {suffix}"), - model=ModelOptions(name=f"model {suffix}"), - active=SyntheticsActiveCellsOptions(name=f"active {suffix}"), - ), + for name in ["Gravity Forward", "Magnetic Vector Forward"]: + group = geoh5.get_entity(name)[0] + mesh = next(child for child in group.children if isinstance(child, Octree)) + survey = next( + child for child in group.children if isinstance(child, Points) ) - mesh = components.mesh - survey = components.survey - if suffix == "A": + if name == "Gravity Forward": global_mesh = mesh.copy(parent=geoh5) model = global_mesh.get_entity("starting_model")[0] @@ -197,7 +189,7 @@ def test_homogeneous_run( ref_model = mesh.get_entity("starting_model")[0].copy(name="ref_model") ref_model.values = ref_model.values / 2.0 - if suffix == "A": + if name == "Gravity Forward": params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, From 607de34b260e629dae1058c25c04639aa8e6a8f4 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 15:07:25 -0700 Subject: [PATCH 18/24] Improve joint cross test --- .../driver_joint_cross_gradient_test.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/tests/run_tests/driver_joint_cross_gradient_test.py b/tests/run_tests/driver_joint_cross_gradient_test.py index d306a53d..2ad636a6 100644 --- a/tests/run_tests/driver_joint_cross_gradient_test.py +++ b/tests/run_tests/driver_joint_cross_gradient_test.py @@ -56,7 +56,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 53.295822303985325, "phi_d": 7940, "phi_m": 0.0255} +target_run = {"data_norm": 53.295822303985325, "phi_d": 646, "phi_m": 1.84} INDUCING_FIELD = (50000.0, 90.0, 0.0) @@ -199,7 +199,7 @@ def test_joint_cross_gradient_inv_run( topography_object=topography, data_object=survey, gz_channel=data, - gz_uncertainty=1e-2, + gz_uncertainty=5e-3, starting_model=0.0, reference_model=0.0, upper_bound=1.0, @@ -209,6 +209,13 @@ def test_joint_cross_gradient_inv_run( ) drivers.append(GravityInversionDriver(params)) elif name == "Direct Current 3D Forward": + uncertainties = survey.add_data( + { + "Uncertainties": { + "values": np.abs(data.values) * 0.05 + 1e-4, + } + } + ) params = DC3DInversionOptions.build( geoh5=geoh5, mesh=mesh, @@ -217,7 +224,7 @@ def test_joint_cross_gradient_inv_run( data_object=survey, potential_channel=data, model_type="Resistivity (Ohm-m)", - potential_uncertainty=5e-4, + potential_uncertainty=uncertainties, tile_spatial=1, starting_model=100.0, reference_model=100.0, @@ -255,7 +262,7 @@ def test_joint_cross_gradient_inv_run( group_c=drivers[2].out_group, group_c_multiplier=1.0, max_global_iterations=max_iterations, - initial_beta_ratio=1e1, + initial_beta_ratio=1e-1, cross_gradient_weight_a_b=1e0, cross_gradient_weight_c_a=1e0, cross_gradient_weight_c_b=1e0, @@ -280,13 +287,13 @@ def test_joint_cross_gradient_inv_run( # the scaling from its total misfit. np.testing.assert_allclose( driver.directives.scale_misfits.scalings, - [0.5011, 0.5, 0.5, 0.5, 1.0], + [0.52747, 0.5, 0.5, 0.5, 1.0], atol=1e-3, ) # Check that scaling * chi factor is reflected in data misfit multipliers np.testing.assert_allclose( driver.data_misfit.multipliers, - [0.4009, 0.4, 0.5, 0.5, 1.0], + [0.421978, 0.4, 0.5, 0.5, 1.0], atol=1e-3, ) From ce4d125defa9629ca703a481e5246cd1e9dffa26 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 15:25:24 -0700 Subject: [PATCH 19/24] Re-instate mechanics of Synthetic --- simpeg_drivers/utils/synthetics/driver.py | 59 ++++++++++++++--------- 1 file changed, 37 insertions(+), 22 deletions(-) diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index 50fc87cf..4aa72711 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -45,51 +45,66 @@ def __init__( @property def topography(self) -> Surface: if self._topography is None: - entity = get_topography_surface( - geoh5=self.geoh5, - options=self.options, - ) + entity = self.geoh5.get_entity("topography")[0] + + if entity is None: + entity = get_topography_surface( + geoh5=self.geoh5, + options=self.options, + ) self._topography = entity return self._topography @property def survey(self) -> ObjectBase: if self._survey is None: - entity = get_survey( - geoh5=self.geoh5, - method=self.options.method, - options=self.options.survey, - ) + entity = self.geoh5.get_entity(self.options.survey.name)[0] + + if entity is None: + entity = get_survey( + geoh5=self.geoh5, + method=self.options.method, + options=self.options.survey, + ) self._survey = entity return self._survey @property def mesh(self) -> Octree | DrapeModel: if self._mesh is None: - entity = get_mesh( - self.options.method, - survey=self.survey, - topography=self.topography, - options=self.options.mesh, - ) + entity = self.geoh5.get_entity("mesh")[0] + + if entity is None: + entity = get_mesh( + self.options.method, + survey=self.survey, + topography=self.topography, + options=self.options.mesh, + ) self._mesh = entity return self._mesh @property def active(self) -> FloatData: if self._active is None: - entity = get_active(self.mesh, self.topography) + entity = self.mesh.get_entity(self.options.active.name)[0] + + if entity is None: + entity = get_active(self.mesh, self.topography) self._active = entity + return self._active @property def model(self) -> FloatData: if self._model is None: - entity = get_model( - method=self.options.method, - mesh=self.mesh, - active=self.active.values, - options=self.options.model, - ) + entity = self.mesh.get_entity(self.options.model.name)[0] + if entity is None: + entity = get_model( + method=self.options.method, + mesh=self.mesh, + active=self.active.values, + options=self.options.model, + ) self._model = entity return self._model From 86d77051bafbaa48251044fc20528e24278bcf3a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 10 Mar 2026 22:26:43 +0000 Subject: [PATCH 20/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- simpeg_drivers/utils/synthetics/driver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py index 4aa72711..f2a7f03e 100644 --- a/simpeg_drivers/utils/synthetics/driver.py +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -92,7 +92,7 @@ def active(self) -> FloatData: if entity is None: entity = get_active(self.mesh, self.topography) self._active = entity - + return self._active @property From 836320837fda1447b1ac78037bde1b3d13b4c931 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 10 Mar 2026 17:50:01 -0700 Subject: [PATCH 21/24] Fix test --- tests/run_tests/driver_ground_tem_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/run_tests/driver_ground_tem_test.py b/tests/run_tests/driver_ground_tem_test.py index 1a396871..26b3b9e4 100644 --- a/tests/run_tests/driver_ground_tem_test.py +++ b/tests/run_tests/driver_ground_tem_test.py @@ -252,7 +252,7 @@ def test_ground_tem_run(tmp_path: Path, max_iterations=1, pytest=True): output = get_inversion_output( driver.params.geoh5.h5file, driver.params.out_group.uid ) - assert driver.inversion_data.entity.tx_id_property.name == "Transmitter ID" + assert driver.inversion_data.entity.tx_id_property.name == "tx_id" output["data"] = orig_dBzdt if pytest: check_target(output, target_run) From 88c242a6d4e44a57108d8eac51907ecd75f55ce3 Mon Sep 17 00:00:00 2001 From: domfournier Date: Wed, 11 Mar 2026 13:57:45 -0700 Subject: [PATCH 22/24] Update simpeg_drivers/utils/utils.py Co-authored-by: benk-mira <81254271+benk-mira@users.noreply.github.com> --- simpeg_drivers/utils/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 3d7fc4dd..e3967a18 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -627,7 +627,7 @@ def simpeg_group_to_driver(group: SimPEGGroup, workspace: Workspace) -> Inversio return inversion_driver(params) -def compute_alongline_distance(points: np.ndarray, ordered: bool = True): +def compute_alongline_distance(points: np.ndarray, ordered: bool = True) -> np.ndarray: """ Convert from cartesian (x, y, values) points to (distance, values) locations. From 03c8adf7f1f03fce9385aaf46b95a2ea19042fb8 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 11 Mar 2026 14:23:53 -0700 Subject: [PATCH 23/24] Flatten options for ip --- .../three_dimensions/options.py | 2 +- .../two_dimensions/options.py | 7 ++++-- simpeg_drivers/electricals/options.py | 25 ------------------- simpeg_drivers/options.py | 8 ++++++ 4 files changed, 14 insertions(+), 28 deletions(-) delete mode 100644 simpeg_drivers/electricals/options.py diff --git a/simpeg_drivers/electricals/induced_polarization/three_dimensions/options.py b/simpeg_drivers/electricals/induced_polarization/three_dimensions/options.py index 739b099a..780e256a 100644 --- a/simpeg_drivers/electricals/induced_polarization/three_dimensions/options.py +++ b/simpeg_drivers/electricals/induced_polarization/three_dimensions/options.py @@ -18,10 +18,10 @@ from geoh5py.objects.surveys.direct_current import PotentialElectrode from simpeg_drivers import assets_path -from simpeg_drivers.electricals.options import IPModelOptions from simpeg_drivers.options import ( BaseForwardOptions, BaseInversionOptions, + IPModelOptions, ) diff --git a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py index e3b9f5dd..b92fcfcd 100644 --- a/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py +++ b/simpeg_drivers/electricals/induced_polarization/two_dimensions/options.py @@ -18,8 +18,11 @@ from simpeg_drivers import assets_path from simpeg_drivers.electricals.base_2d import Base2DOptions -from simpeg_drivers.electricals.options import IPModelOptions -from simpeg_drivers.options import BaseForwardOptions, BaseInversionOptions +from simpeg_drivers.options import ( + BaseForwardOptions, + BaseInversionOptions, + IPModelOptions, +) class IP2DForwardOptions(BaseForwardOptions, Base2DOptions): diff --git a/simpeg_drivers/electricals/options.py b/simpeg_drivers/electricals/options.py deleted file mode 100644 index 7cca3219..00000000 --- a/simpeg_drivers/electricals/options.py +++ /dev/null @@ -1,25 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2023-2026 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from __future__ import annotations - -from geoh5py.data import FloatData -from pydantic import BaseModel - -from simpeg_drivers.options import ConductivityModelOptions - - -class IPModelOptions(ConductivityModelOptions): - """ - ModelOptions class with defaulted lower bound. - """ - - lower_bound: float | FloatData | None = 0 diff --git a/simpeg_drivers/options.py b/simpeg_drivers/options.py index 556c11a2..c25589c5 100644 --- a/simpeg_drivers/options.py +++ b/simpeg_drivers/options.py @@ -632,3 +632,11 @@ def component_uncertainty(self, component: str) -> np.ndarray | None: data *= np.ones_like(self.component_data(component)[None]) return {None: data} + + +class IPModelOptions(ConductivityModelOptions): + """ + ModelOptions class with defaulted lower bound. + """ + + lower_bound: float | FloatData | None = 0 From 9c0ab8dda9078a7e07ff056772354f271d7ea445 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 11 Mar 2026 14:25:33 -0700 Subject: [PATCH 24/24] Docstrings and return types in utils --- simpeg_drivers/utils/utils.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 3d7fc4dd..19462324 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -351,7 +351,7 @@ def cell_width_from_centers(centers: np.ndarray) -> np.ndarray: return np.r_[half_dx[0] * 2, (half_dx[:-1] + half_dx[1:]), half_dx[-1] * 2] -def floating_active(mesh: TensorMesh | TreeMesh, active: np.ndarray): +def floating_active(mesh: TensorMesh | TreeMesh, active: np.ndarray) -> bool: """ True if there are any active cells in the air @@ -482,6 +482,8 @@ def get_containing_cells( :param mesh: Computational mesh object :param data: Inversion data object + + :returns: Array of unique cell indices that contain data locations """ if isinstance(mesh, TreeMesh): if isinstance(data.entity, PotentialElectrode): @@ -537,13 +539,15 @@ def active_from_xyz( topo: np.ndarray, grid_reference="center", triangulation: np.ndarray | None = None, -): +) -> np.ndarray: """Returns an active cell index array below a surface :param mesh: Mesh object :param topo: Array of xyz locations :param grid_reference: Cell reference. Must be "center", "top", or "bottom" :param method: Interpolation method. Must be "linear", or "nearest" + + :return: Array of active cell indices below the surface defined by 'topo'. """ mesh_dim = 2 if isinstance(mesh, DrapeModel) else 3 @@ -606,6 +610,8 @@ def simpeg_group_to_driver(group: SimPEGGroup, workspace: Workspace) -> Inversio :param group: SimPEGGroup object. :param workspace: Workspace object. + + :returns: InversionDriver object. """ ui_json = deepcopy(group.options) @@ -631,8 +637,15 @@ def compute_alongline_distance(points: np.ndarray, ordered: bool = True): """ Convert from cartesian (x, y, values) points to (distance, values) locations. - :param: points: Cartesian coordinates of points lying either roughly within a + :param points: Cartesian coordinates of points lying either roughly within a plane or a line. + :param ordered: Flag to indicate whether points are already ordered along a line. + If False, then the points will be ordered using a traveling salesman algorithm + before computing distances. + + :returns: Array of shape (n_points, n_features) where the first column contains + distances along the line and the remaining columns contain the corresponding + values from the input 'points' array. """ if not ordered: order = traveling_salesman(points)