From 35b59d69504f1ddb358d6d5838843776043272ee Mon Sep 17 00:00:00 2001 From: dominiquef Date: Mon, 5 Feb 2024 16:36:11 -0800 Subject: [PATCH 01/17] Force Cell association to data on dc constants --- .../direct_current/pseudo_three_dimensions/constants.py | 6 +++--- .../electricals/direct_current/two_dimensions/constants.py | 6 +++--- .../pseudo_three_dimensions/constants.py | 6 +++--- .../induced_polarization/two_dimensions/constants.py | 6 +++--- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py index 1d25c5c4e..8b26b059d 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py @@ -141,7 +141,7 @@ "icon": "PotentialElectrode", "inversion_type": "direct current pseudo 3d", "line_object": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Referenced", "group": "Data", "main": True, @@ -159,7 +159,7 @@ }, "potential_channel_bool": True, "potential_channel": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, @@ -168,7 +168,7 @@ "value": None, }, "potential_uncertainty": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, diff --git a/geoapps/inversion/electricals/direct_current/two_dimensions/constants.py b/geoapps/inversion/electricals/direct_current/two_dimensions/constants.py index 59ed244fd..cc8db9774 100644 --- a/geoapps/inversion/electricals/direct_current/two_dimensions/constants.py +++ b/geoapps/inversion/electricals/direct_current/two_dimensions/constants.py @@ -135,7 +135,7 @@ "icon": "PotentialElectrode", "inversion_type": "direct current 2d", "line_object": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Referenced", "group": "Data", "main": True, @@ -166,7 +166,7 @@ }, "potential_channel_bool": True, "potential_channel": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, @@ -176,7 +176,7 @@ "value": None, }, "potential_uncertainty": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, diff --git a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py index 837f6dc13..ecbf23556 100644 --- a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py @@ -143,7 +143,7 @@ "icon": "PotentialElectrode", "inversion_type": "induced polarization pseudo 3d", "line_object": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Referenced", "group": "Data", "main": True, @@ -161,7 +161,7 @@ }, "chargeability_channel_bool": True, "chargeability_channel": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, @@ -170,7 +170,7 @@ "value": None, }, "chargeability_uncertainty": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, diff --git a/geoapps/inversion/electricals/induced_polarization/two_dimensions/constants.py b/geoapps/inversion/electricals/induced_polarization/two_dimensions/constants.py index 44aad2675..2ff8d634b 100644 --- a/geoapps/inversion/electricals/induced_polarization/two_dimensions/constants.py +++ b/geoapps/inversion/electricals/induced_polarization/two_dimensions/constants.py @@ -136,7 +136,7 @@ "icon": "PotentialElectrode", "inversion_type": "induced polarization 2d", "line_object": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Referenced", "group": "Data", "main": True, @@ -167,7 +167,7 @@ }, "chargeability_channel_bool": True, "chargeability_channel": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, @@ -176,7 +176,7 @@ "value": None, }, "chargeability_uncertainty": { - "association": ["Cell", "Vertex"], + "association": "Cell", "dataType": "Float", "group": "Data", "main": True, From c3b6556bd96c0ef66aabbeac032190270d11c566 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 13:19:22 -0800 Subject: [PATCH 02/17] Move transfer model to seperate method --- .../pseudo_three_dimensions/driver.py | 71 ++++++++++--------- 1 file changed, 36 insertions(+), 35 deletions(-) diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py index 806ad7520..5c6a59c18 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py @@ -8,11 +8,13 @@ from __future__ import annotations import sys +import uuid from copy import deepcopy from pathlib import Path import numpy as np from geoh5py.data import Data +from geoh5py.objects import DrapeModel from geoh5py.workspace import Workspace from geoapps.inversion.components.data import InversionData @@ -42,6 +44,33 @@ def __init__(self, params: DirectCurrentPseudo3DParams): # pylint: disable=W023 if params.files_only: sys.exit("Files written") + def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID]: + xyz_in = get_locations(self.workspace, self.pseudo3d_params.mesh) + models = {"starting_model": self.pseudo3d_params.starting_model} + if not self.pseudo3d_params.forward_only: + models.update( + { + "reference_model": self.pseudo3d_params.reference_model, + "lower_bound": self.pseudo3d_params.lower_bound, + "upper_bound": self.pseudo3d_params.upper_bound, + } + ) + + xyz_out = mesh.centroids + model_uids = deepcopy(models) + 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}}) + model_uids[name] = model_object.uid + + return model_uids + def write_files(self, lookup): """Write ui.geoh5 and ui.json files for sweep trials.""" @@ -54,22 +83,12 @@ def write_files(self, lookup): self._inversion_topography = InversionTopography( self.workspace, self.pseudo3d_params ) - xyz_in = get_locations(self.workspace, self.pseudo3d_params.mesh) - models = {"starting_model": self.pseudo3d_params.starting_model} - if not forward_only: - models.update( - { - "reference_model": self.pseudo3d_params.reference_model, - "lower_bound": self.pseudo3d_params.lower_bound, - "upper_bound": self.pseudo3d_params.upper_bound, - } - ) - - for uuid, trial in lookup.items(): + + for uid, trial in lookup.items(): if trial["status"] != "pending": continue - filepath = Path(self.working_directory) / f"{uuid}.ui.geoh5" + filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" with Workspace(filepath) as iter_workspace: receiver_entity = extract_dcip_survey( iter_workspace, @@ -96,23 +115,7 @@ def write_files(self, lookup): self.pseudo3d_params.expansion_factor, )[0] - iter_workspace.remove_entity(receiver_entity.current_electrodes) - iter_workspace.remove_entity(receiver_entity) - - xyz_out = mesh.centroids - model_uids = deepcopy(models) - 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}}) - model_uids[name] = model_object.uid + model_uids = self.transfer_models(mesh) for key in ifile.data: param = getattr(self.pseudo3d_params, key, None) @@ -122,15 +125,13 @@ def write_files(self, lookup): self.pseudo3d_params.topography_object.copy( parent=iter_workspace, copy_children=True ) - self.pseudo3d_params.data_object.copy( - parent=iter_workspace, copy_children=True - ) ifile.data.update( dict( **{ "geoh5": iter_workspace, "mesh": mesh, + "data_object": receiver_entity, "line_id": trial["line_id"], "out_group": None, }, @@ -138,9 +139,9 @@ def write_files(self, lookup): ) ) - ifile.name = f"{uuid}.ui.json" + ifile.name = f"{uid}.ui.json" ifile.path = self.working_directory # pylint: disable=E1101 ifile.write_ui_json() - lookup[uuid]["status"] = "written" + lookup[uid]["status"] = "written" _ = self.update_lookup(lookup) # pylint: disable=E1101 From eebc15956c0dd1395b7eb8ece3dc4faf27ab1c68 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 13:20:43 -0800 Subject: [PATCH 03/17] Remove global map, simplify logic for 2D --- .../components/factories/survey_factory.py | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/geoapps/inversion/components/factories/survey_factory.py b/geoapps/inversion/components/factories/survey_factory.py index 5f421eff8..8d91ec127 100644 --- a/geoapps/inversion/components/factories/survey_factory.py +++ b/geoapps/inversion/components/factories/survey_factory.py @@ -280,10 +280,6 @@ def _dcip_arguments(self, data=None, local_index=None): self.params.line_object.values, self.params.line_id, ) - self.local_index = np.arange(receiver_entity.n_cells) - data.global_map = [ - k for k in receiver_entity.children if k.name == "Global Map" - ][0].values source_ids, order = np.unique( receiver_entity.ab_cell_id.values[self.local_index], return_index=True @@ -291,12 +287,8 @@ def _dcip_arguments(self, data=None, local_index=None): currents = receiver_entity.current_electrodes if "2d" in self.params.inversion_type: - receiver_locations = receiver_entity.vertices - source_locations = currents.vertices - if local_index is not None: - receiver_locations = data.drape_locations(receiver_locations) - source_locations = data.drape_locations(source_locations) - + receiver_locations = data.drape_locations(receiver_entity.vertices) + source_locations = data.drape_locations(currents.vertices) else: receiver_locations = data.locations source_locations = currents.vertices @@ -332,11 +324,6 @@ def _dcip_arguments(self, data=None, local_index=None): self.local_index = np.hstack(self.local_index) - if "2d" in self.factory_type: - current_entity = receiver_entity.current_electrodes - self.params.geoh5.remove_entity(receiver_entity) - self.params.geoh5.remove_entity(current_entity) - return [sources] def _tdem_arguments(self, data=None, local_index=None, mesh=None): From 68afee72fc94a8541bc4ec73f4ff1c4479486a0c Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 13:26:19 -0800 Subject: [PATCH 04/17] Rely on geoh5py to extract survey from indices. Fix issue with cell ordering --- geoapps/utils/surveys.py | 122 +++++++++------------------------------ 1 file changed, 26 insertions(+), 96 deletions(-) diff --git a/geoapps/utils/surveys.py b/geoapps/utils/surveys.py index f64569d1f..cc812a141 100644 --- a/geoapps/utils/surveys.py +++ b/geoapps/utils/surveys.py @@ -17,8 +17,7 @@ import numpy as np from discretize import TensorMesh, TreeMesh -from geoh5py.data import FloatData -from geoh5py.objects import CurrentElectrode, PotentialElectrode +from geoh5py.objects import PotentialElectrode from geoh5py.objects.surveys.electromagnetics.base import LargeLoopGroundEMSurvey from geoh5py.workspace import Workspace from scipy.spatial import cKDTree @@ -56,7 +55,7 @@ def get_containing_cells( locations = data.drape_locations(get_unique_locations(data.survey)) xi = np.searchsorted(mesh.nodes_x, locations[:, 0]) - 1 yi = np.searchsorted(mesh.nodes_y, locations[:, -1]) - 1 - inds = xi + yi * mesh.shape_cells[0] + inds = xi * mesh.shape_cells[1] + yi else: raise TypeError("Mesh must be 'TreeMesh' or 'TensorMesh'") @@ -165,7 +164,7 @@ def compute_alongline_distance(points: np.ndarray, ordered: bool = True): def survey_lines(survey, start_loc: list[int | float], save: str | None = None): """ - Build an array of line ids for a survey laid out in a line biased grid. + Build an array of line ids for a dc survey laid out in a line biased grid. :param: survey: geoh5py.objects.surveys object with .vertices attribute or xyz array. :param: start_loc: Easting and Northing of a survey extremity from which the @@ -173,12 +172,10 @@ def survey_lines(survey, start_loc: list[int | float], save: str | None = None): :save: Name assigned to line id (ReferencedData) object if not None. """ - # extract xy locations and create linear indexing - try: - locs = survey.vertices[:, :2] - except AttributeError: - locs = survey[:, :2] + if survey.vertices is None: + raise ValueError("Survey object must have a vertices attribute.") + locs = survey.vertices[:, :2] nodes = np.arange(len(locs)).tolist() # find the id of the closest point to the starting location @@ -221,23 +218,32 @@ def survey_lines(survey, start_loc: list[int | float], save: str | None = None): loc = locs[next_id] lines += [line_id] # nodes run out before last id assigned - inds = np.argsort(inds) + sorted_id = np.array(lines)[np.argsort(inds)] + + # Convert to cell data + cell_ids = sorted_id[survey.cells] + + assert np.all(cell_ids[:, 0] == cell_ids[:, 1]) + + line_reference = {0: "Unknown"} + line_reference.update({k: str(k) for k in np.unique(lines)}) + if save is not None: survey.add_data( { save: { - "values": np.array(lines)[inds], - "association": "VERTEX", + "values": cell_ids[:, 0], + "association": "CELL", "entity_type": { "primitive_type": "REFERENCED", - "value_map": {k: str(k) for k in lines}, + "value_map": line_reference, }, } } ) - return np.array(lines)[inds] + return sorted_id def slice_and_map(obj: np.ndarray, slicer: np.ndarray | Callable): @@ -272,7 +278,6 @@ def extract_dcip_survey( survey: PotentialElectrode, lines: np.ndarray, line_id: int, - name: str = "Line", ): """ Returns a survey containing data from a single line. @@ -281,93 +286,18 @@ def extract_dcip_survey( :param: survey: PotentialElectrode object. :param: lines: Line indexer for survey. :param: line_id: Index of line to extract data from. - :param: name: Name prefix to assign to the new survey object (to be - suffixed with the line number). """ - current = survey.current_electrodes - - # Extract line locations and store map into full survey - survey_locs, survey_loc_map = slice_and_map(survey.vertices, lines == line_id) - - # Use line locations to slice cells and store map into full survey - func = lambda c: (c[0] in survey_loc_map) & (c[1] in survey_loc_map) - survey_cells, survey_cell_map = slice_and_map(survey.cells, func) - survey_cells = [[survey_loc_map[i] for i in c] for c in survey_cells] - - # Use line cells to slice ab_cell_ids - ab_cell_ids = survey.ab_cell_id.values[list(survey_cell_map)] - ab_cell_ids = np.array(ab_cell_ids, dtype=int) - 1 - - # Use line ab_cell_ids to slice current cells - current_cells, current_cell_map = slice_and_map( - current.cells, np.unique(ab_cell_ids) - ) - - # Use line current cells to slice current locs - current_locs, current_loc_map = slice_and_map( - current.vertices, np.unique(current_cells.ravel()) - ) - - # Remap global ids to local counterparts - ab_cell_ids = np.array([current_cell_map[i] for i in ab_cell_ids]) - current_cells = [[current_loc_map[i] for i in c] for c in current_cells] - - # Save objects - line_name = f"{name} {line_id}" - currents = CurrentElectrode.create( - workspace, - name=f"{line_name} (currents)", - vertices=current_locs, - cells=np.array(current_cells, dtype=np.int32), - allow_delete=True, - ) - currents.add_default_ab_cell_id() - - potentials = PotentialElectrode.create( - workspace, - name=line_name, - vertices=survey_locs, - allow_delete=True, - ) - potentials.cells = np.array(survey_cells, dtype=np.int32) - - # Add ab_cell_id as referenced data object - value_map = {k + 1: str(k + 1) for k in ab_cell_ids} - value_map.update({0: "Unknown"}) - potentials.add_data( - { - "A-B Cell ID": { - "values": np.array(ab_cell_ids + 1, dtype=np.int32), - "association": "CELL", - "entity_type": { - "primitive_type": "REFERENCED", - "value_map": value_map, - }, - }, - "Global Map": { - "values": np.array(list(survey_cell_map), dtype=np.int32), - "association": "CELL", - "entity_type": { - "primitive_type": "REFERENCED", - "value_map": {k: str(k) for k in survey_cell_map}, - }, - }, - } - ) - - # Attach current and potential objects and copy data slice into line survey - potentials.current_electrodes = currents - for c in survey.children: - if isinstance(c, FloatData) and "Pseudo" not in c.name: - potentials.add_data({c.name: {"values": c.values[list(survey_cell_map)]}}) + cell_mask = lines == line_id + 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 split_dcip_survey( +def plit_dcip_surveys( survey: PotentialElectrode, lines: np.ndarray, - name: str = "Line", workspace: Workspace | None = None, ): """ @@ -386,7 +316,7 @@ def split_dcip_survey( with ws.open(mode="r+") as ws: line_surveys = [] for line_id in np.unique(lines): - line_survey = extract_dcip_survey(ws, survey, lines, line_id, name) + line_survey = extract_dcip_survey(ws, survey, lines, line_id) line_surveys.append(line_survey) return line_surveys From fb4deab7a62f87609907065ee9dffdf55ef233be Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 15:22:23 -0800 Subject: [PATCH 05/17] Fix line length --- geoapps/inversion/components/topography.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/geoapps/inversion/components/topography.py b/geoapps/inversion/components/topography.py index a039457f5..f1f24c82e 100644 --- a/geoapps/inversion/components/topography.py +++ b/geoapps/inversion/components/topography.py @@ -102,7 +102,8 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: if floating_active(mesh.mesh, active_cells): warnings.warn( - "Active cell adjustment has created a patch of active cells in the air, likely due to a faulty survey location." + "Active cell adjustment has created a patch of active cells in the air, " + "likely due to a faulty survey location." ) return active_cells From cdc2880fe880c7dd8810b52c0c27d7d243741b55 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 15:24:02 -0800 Subject: [PATCH 06/17] Refactor line logic in collecting data --- geoapps/inversion/line_sweep/driver.py | 32 ++++++++++++++------------ 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/geoapps/inversion/line_sweep/driver.py b/geoapps/inversion/line_sweep/driver.py index e7cb3925d..f3b925c72 100644 --- a/geoapps/inversion/line_sweep/driver.py +++ b/geoapps/inversion/line_sweep/driver.py @@ -116,17 +116,19 @@ def line_files(path: str | Path): def collect_results(self): path = Path(self.workspace.h5file).parent files = LineSweepDriver.line_files(str(path)) - lines = np.unique(self.pseudo3d_params.line_object.values) - data_result = self.pseudo3d_params.data_object.copy( - parent=self.pseudo3d_params.out_group - ) - + line_ids = self.pseudo3d_params.line_object.values data = {} drape_models = [] - for line in lines: + for line in np.unique(line_ids): with Workspace(f"{path / files[line]}.ui.geoh5") as ws: survey = ws.get_entity("Data")[0] - data = self.collect_line_data(survey, data) + line_data = survey.get_entity(self.pseudo3d_params.line_object.name) + + if not line_data: + raise ValueError(f"Line {line} not found in {survey.name}") + + line_indices = line_ids == line + data = self.collect_line_data(survey, line_indices, data) mesh = ws.get_entity("Models")[0] filedata = [ @@ -143,7 +145,7 @@ def collect_results(self): mesh.name = "models" drape_models.append(mesh) - data_result.add_data(data) + self.pseudo3d_params.data_object.add_data(data) if self.pseudo3d_params.mesh is None: return @@ -186,16 +188,16 @@ def collect_results(self): octree_model.copy(parent=self.pseudo3d_params.out_group) - def collect_line_data(self, survey, data): + def collect_line_data(self, survey, line_indices, data): + """ + Fill chunks of values from one line + """ for child in survey.children: # initialize data values dictionary if "Iteration" in child.name and child.name not in data: - data[child.name] = {"values": np.zeros(survey.n_cells)} + data[child.name] = {"values": np.zeros_like(line_indices) * np.nan} - ind = None - for child in survey.children: # fill a chunk of values from one line + for child in survey.children: if "Iteration" in child.name: - if ind is None: - ind = ~np.isnan(child.values) - data[child.name]["values"][ind] = child.values[ind] + data[child.name]["values"][line_indices] = child.values return data From 500d805ff704a3812640ed8ce894aa66c42c417a Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 6 Feb 2024 15:32:54 -0800 Subject: [PATCH 07/17] Add dc survey generation. Update test --- geoapps/utils/surveys.py | 92 +----------------- geoapps/utils/testing.py | 109 +++++++++++---------- tests/run_tests/driver_dc_p3d_test.py | 8 +- tests/utils_test.py | 134 ++++---------------------- 4 files changed, 85 insertions(+), 258 deletions(-) diff --git a/geoapps/utils/surveys.py b/geoapps/utils/surveys.py index cc812a141..98bb8349f 100644 --- a/geoapps/utils/surveys.py +++ b/geoapps/utils/surveys.py @@ -22,8 +22,6 @@ from geoh5py.workspace import Workspace from scipy.spatial import cKDTree -from geoapps.utils.statistics import is_outlier - def get_containing_cells( mesh: TreeMesh | TensorMesh, data: InversionData @@ -162,90 +160,6 @@ def compute_alongline_distance(points: np.ndarray, ordered: bool = True): return distances -def survey_lines(survey, start_loc: list[int | float], save: str | None = None): - """ - Build an array of line ids for a dc survey laid out in a line biased grid. - - :param: survey: geoh5py.objects.surveys object with .vertices attribute or xyz array. - :param: start_loc: Easting and Northing of a survey extremity from which the - all other survey locations will be traversed and assigned line ids. - :save: Name assigned to line id (ReferencedData) object if not None. - - """ - if survey.vertices is None: - raise ValueError("Survey object must have a vertices attribute.") - - locs = survey.vertices[:, :2] - nodes = np.arange(len(locs)).tolist() - - # find the id of the closest point to the starting location - start_id = np.argmin(np.linalg.norm(locs - start_loc, axis=1)) - - # pop the starting location and index out of their respective lists - locs = locs.tolist() - loc = locs[start_id] - inds = [] - n = nodes.pop(start_id) - inds.append(n) - - # compute the tree of the remaining points and begin to traverse the tree - # in the direction of closest neighbors. Label points with same line id - # until an outlier is detected in the distance to the next closest point, - # then increase the line id. - tree = cKDTree(locs) - line_id = 1 # zero is reserved - lines = [] - distances = [] - while nodes: - lines.append(line_id) - dist, next_id = next_neighbor(tree, loc, nodes) - - outlier = False - if len(distances) > 1: - if np.allclose(dist, distances, atol=1e-6): - outlier = False - else: - outlier = is_outlier(distances, dist) - - if outlier: - line_id += 1 - distances = [] - else: - distances.append(dist) - - n = nodes.pop(nodes.index(next_id)) - inds.append(n) - loc = locs[next_id] - - lines += [line_id] # nodes run out before last id assigned - inds = np.argsort(inds) - sorted_id = np.array(lines)[np.argsort(inds)] - - # Convert to cell data - cell_ids = sorted_id[survey.cells] - - assert np.all(cell_ids[:, 0] == cell_ids[:, 1]) - - line_reference = {0: "Unknown"} - line_reference.update({k: str(k) for k in np.unique(lines)}) - - if save is not None: - survey.add_data( - { - save: { - "values": cell_ids[:, 0], - "association": "CELL", - "entity_type": { - "primitive_type": "REFERENCED", - "value_map": line_reference, - }, - } - } - ) - - return sorted_id - - def slice_and_map(obj: np.ndarray, slicer: np.ndarray | Callable): """ Slice an array and return both sliced array and global to local map. @@ -288,6 +202,10 @@ def extract_dcip_survey( :param: line_id: Index of line to extract data from. """ cell_mask = lines == line_id + + if not np.any(cell_mask): + raise ValueError(f"Line '{line_id}' not found in survey.") + 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) @@ -295,7 +213,7 @@ def extract_dcip_survey( return potentials -def plit_dcip_surveys( +def split_dcip_survey( survey: PotentialElectrode, lines: np.ndarray, workspace: Workspace | None = None, diff --git a/geoapps/utils/testing.py b/geoapps/utils/testing.py index 638b86263..553284cd8 100644 --- a/geoapps/utils/testing.py +++ b/geoapps/utils/testing.py @@ -36,7 +36,6 @@ from geoapps.driver_base.utils import active_from_xyz, treemesh_2_octree from geoapps.octree_creation.driver import OctreeDriver from geoapps.utils.models import get_drape_model -from geoapps.utils.surveys import survey_lines class Geoh5Tester: @@ -78,6 +77,63 @@ def make(self): return self.ws +def generate_dc_survey(workspace, x_loc, y_loc, z_loc=None): + """ + Utility function to generate a DC survey. + """ + # Create sources along line + if z_loc is None: + z_loc = np.zeros_like(x_loc) + + vertices = np.c_[x_loc.ravel(), y_loc.ravel(), z_loc.ravel()] + parts = np.kron(np.arange(x_loc.shape[0]), np.ones(x_loc.shape[1])).astype("int") + currents = CurrentElectrode.create(workspace, vertices=vertices, parts=parts) + currents.add_default_ab_cell_id() + n_dipoles = 9 + dipoles = [] + current_id = [] + + for val in currents.ab_cell_id.values: + cell_id = int(currents.ab_map[val]) - 1 + + for dipole in range(n_dipoles): + dipole_ids = currents.cells[cell_id, :] + 2 + dipole + + if ( + any(dipole_ids > (currents.n_vertices - 1)) + or len( + np.unique(parts[np.r_[currents.cells[cell_id, 0], dipole_ids[1]]]) + ) + > 1 + ): + continue + + dipoles += [dipole_ids] + current_id += [val] + + potentials = PotentialElectrode.create( + workspace, vertices=vertices, cells=np.vstack(dipoles).astype("uint32") + ) + line_id = potentials.vertices[potentials.cells[:, 0], 1] + line_id = (line_id - np.min(line_id) + 1).astype(np.int32) + line_reference = {0: "Unknown"} + line_reference.update({k: str(k) for k in np.unique(line_id)}) + potentials.add_data( + { + "line_ids": { + "values": line_id, + "type": "REFERENCED", + "value_map": line_reference, + } + } + ) + potentials.ab_cell_id = np.hstack(current_id).astype("int32") + potentials.current_electrodes = currents + currents.potential_electrodes = potentials + + return potentials + + def setup_inversion_workspace( work_dir, background=None, @@ -140,50 +196,7 @@ def topo_drape(x, y): ] if inversion_type in ["dcip", "dcip_2d"]: - ab_vertices = np.c_[ - X[:, :-2].flatten(), Y[:, :-2].flatten(), Z[:, :-2].flatten() - ] - mn_vertices = np.c_[ - X[:, 2:].flatten() + center[0], - Y[:, 2:].flatten() + center[1], - Z[:, 2:].flatten() + center[2], - ] - - parts = np.repeat(np.arange(n_lines), n_electrodes - 2).astype("int32") - currents = CurrentElectrode.create( - geoh5, name="survey (currents)", vertices=ab_vertices, parts=parts - ) - currents.add_default_ab_cell_id() - - N = 6 - dipoles = [] - current_id = [] - for val in currents.ab_cell_id.values: # For each source dipole - cell_id = int(currents.ab_map[val]) - 1 # Python 0 indexing - line = currents.parts[currents.cells[cell_id, 0]] - for m_n in range(N): - dipole_ids = (currents.cells[cell_id, :] + m_n).astype( - "uint32" - ) # Skip two poles - - # Shorten the array as we get to the end of the line - if any(dipole_ids > (len(mn_vertices) - 1)) or any( - currents.parts[dipole_ids] != line - ): - continue - - dipoles += [dipole_ids] # Save the receiver id - current_id += [val] # Save the source id - - survey = PotentialElectrode.create( - geoh5, - name="survey", - vertices=mn_vertices, - cells=np.vstack(dipoles).astype("uint32"), - ) - survey.current_electrodes = currents - survey.ab_cell_id = np.asarray(current_id).astype("int32") - currents.potential_electrodes = survey + survey = generate_dc_survey(geoh5, X, Y, Z) elif inversion_type == "magnetotellurics": survey = MTReceivers.create( @@ -383,13 +396,11 @@ def topo_drape(x, y): # Create a mesh if "2d" in inversion_type: - locs = np.unique(np.vstack([ab_vertices, mn_vertices]), axis=0) - lines = survey_lines(locs, [-100, -100]) - + lines = survey.get_entity("line_id")[0].values entity, mesh, _ = get_drape_model( # pylint: disable=W0632 geoh5, "Models", - locs[lines == 2], + survey.vertices[survey.cells[lines == 2, :], :], [cell_size[0], cell_size[2]], 100.0, [padding_distance] * 2 + [padding_distance] * 2, diff --git a/tests/run_tests/driver_dc_p3d_test.py b/tests/run_tests/driver_dc_p3d_test.py index f2f690f8e..7cba3fcc7 100644 --- a/tests/run_tests/driver_dc_p3d_test.py +++ b/tests/run_tests/driver_dc_p3d_test.py @@ -21,7 +21,6 @@ DirectCurrentPseudo3DParams, ) from geoapps.shared_utils.utils import get_inversion_output -from geoapps.utils.surveys import survey_lines from geoapps.utils.testing import check_target, setup_inversion_workspace # To test the full run and validate the inversion. @@ -50,8 +49,6 @@ def test_dc_p3d_fwr_run( drape_height=0.0, flatten=False, ) - - _ = survey_lines(survey, [-100, -100], save="line_ids") params = DirectCurrentPseudo3DParams( forward_only=True, geoh5=geoh5, @@ -87,7 +84,6 @@ def test_dc_p3d_run( out_group = geoh5.get_entity("Line 1")[0].parent mesh = out_group.get_entity("mesh")[0] # Finds the octree mesh topography = geoh5.get_entity("topography")[0] - _ = survey_lines(potential.parent, [-100, 100], save="line_IDs") # Run the inverse np.random.seed(0) @@ -103,7 +99,7 @@ def test_dc_p3d_run( data_object=potential.parent.uid, potential_channel=potential.uid, potential_uncertainty=1e-3, - line_object=geoh5.get_entity("line_IDs")[0].uid, + line_object=geoh5.get_entity("line_ids")[0].uid, starting_model=1e-2, reference_model=1e-2, s_norm=0.0, @@ -128,7 +124,7 @@ def test_dc_p3d_run( basepath = workpath.parent with open(basepath / "lookup.json", encoding="utf8") as f: lookup = json.load(f) - middle_line_id = [k for k, v in lookup.items() if v["line_id"] == 2][0] + middle_line_id = [k for k, v in lookup.items() if v["line_id"] == 101][0] with Workspace(basepath / f"{middle_line_id}.ui.geoh5", mode="r") as workspace: middle_inversion_group = [ diff --git a/tests/utils_test.py b/tests/utils_test.py index cc49c16de..30da1a766 100644 --- a/tests/utils_test.py +++ b/tests/utils_test.py @@ -18,7 +18,6 @@ from discretize import CylindricalMesh, TreeMesh from discretize.utils import mesh_builder_xyz from geoh5py.objects import Curve, Grid2D, Points -from geoh5py.objects.surveys.direct_current import CurrentElectrode, PotentialElectrode from geoh5py.workspace import Workspace from geoapps.driver_base.utils import active_from_xyz, running_mean, treemesh_2_octree @@ -53,9 +52,8 @@ find_unique_tops, new_neighbors, split_dcip_survey, - survey_lines, ) -from geoapps.utils.testing import Geoh5Tester +from geoapps.utils.testing import Geoh5Tester, generate_dc_survey from geoapps.utils.workspace import sorted_children_dict from . import PROJECT @@ -259,135 +257,39 @@ def test_new_neighbors(): assert neighbor_id[0] == 1 -def test_survey_lines(tmp_path: Path): - name = "TestCurrents" - n_data = 15 - path = tmp_path / r"testDC.geoh5" - workspace = Workspace.create(path) - - # Create sources along line - x_loc, y_loc = np.meshgrid(np.linspace(0, 10, n_data), np.arange(-1, 3)) - vertices = np.c_[x_loc.ravel(), y_loc.ravel(), np.zeros_like(x_loc).ravel()] - parts = np.kron(np.arange(4), np.ones(n_data)).astype("int") - currents = CurrentElectrode.create( - workspace, name=name, vertices=vertices, parts=parts - ) - currents.add_default_ab_cell_id() - potentials = PotentialElectrode.create( - workspace, name=name + "_rx", vertices=vertices - ) - n_dipoles = 9 - dipoles = [] - current_id = [] - for val in currents.ab_cell_id.values: - cell_id = int(currents.ab_map[val]) - 1 - - for dipole in range(n_dipoles): - dipole_ids = currents.cells[cell_id, :] + 2 + dipole - - if ( - any(dipole_ids > (potentials.n_vertices - 1)) - or len(np.unique(parts[dipole_ids])) > 1 - ): - continue - - dipoles += [dipole_ids] - current_id += [val] - - potentials.cells = np.vstack(dipoles).astype("uint32") - potentials.ab_cell_id = np.hstack(current_id).astype("int32") - potentials.current_electrodes = currents - currents.potential_electrodes = potentials - _ = survey_lines(potentials, [-1, -1], save="test_line_ids") - assert np.all( - np.unique(workspace.get_entity("test_line_ids")[0].values) == np.arange(1, 5) - ) - - def test_extract_dcip_survey(tmp_path: Path): - name = "TestCurrents" n_data = 12 path = tmp_path / r"testDC.geoh5" - workspace = Workspace.create(path) - # Create sources along line x_loc, y_loc = np.meshgrid(np.arange(n_data), np.arange(-1, 3)) - vertices = np.c_[x_loc.ravel(), y_loc.ravel(), np.zeros_like(x_loc).ravel()] - parts = np.kron(np.arange(4), np.ones(n_data)).astype("int") - currents = CurrentElectrode.create( - workspace, name=name, vertices=vertices, parts=parts - ) - currents.add_default_ab_cell_id() - potentials = PotentialElectrode.create( - workspace, name=name + "_rx", vertices=vertices - ) - n_dipoles = 9 - dipoles = [] - current_id = [] - for val in currents.ab_cell_id.values: - cell_id = int(currents.ab_map[val]) - 1 - for dipole in range(n_dipoles): - dipole_ids = currents.cells[cell_id, :] + 2 + dipole + with Workspace.create(path) as workspace: + potentials = generate_dc_survey(workspace, x_loc, y_loc) + + line_id = potentials.get_data("line_ids")[0].values - if ( - any(dipole_ids > (potentials.n_vertices - 1)) - or len(np.unique(parts[dipole_ids])) > 1 - ): - continue + with pytest.raises(ValueError, match="Line '3' not found in survey."): + extract_dcip_survey(workspace, potentials, line_id, 3) - dipoles += [dipole_ids] - current_id += [val] + surveys = extract_dcip_survey(workspace, potentials, line_id, 6) - potentials.cells = np.vstack(dipoles).astype("uint32") - potentials.ab_cell_id = np.hstack(current_id).astype("int32") - potentials.current_electrodes = currents - currents.potential_electrodes = potentials - extract_dcip_survey(workspace, potentials, parts, 3, "test_survey_line") - assert workspace.get_entity("test_survey_line 3")[0] is not None + line_field = surveys.get_entity("line_ids") + assert line_field + assert np.all(line_field[0].values == 6) def test_split_dcip_survey(tmp_path: Path): - name = "TestCurrents" n_data = 12 path = tmp_path / r"testDC.geoh5" - workspace = Workspace.create(path) - # Create sources along line x_loc, y_loc = np.meshgrid(np.arange(n_data), np.arange(-1, 3)) - vertices = np.c_[x_loc.ravel(), y_loc.ravel(), np.zeros_like(x_loc).ravel()] - parts = np.kron(np.arange(4), np.ones(n_data)).astype("int") - currents = CurrentElectrode.create( - workspace, name=name, vertices=vertices, parts=parts - ) - currents.add_default_ab_cell_id() - potentials = PotentialElectrode.create( - workspace, name=name + "_rx", vertices=vertices - ) - n_dipoles = 9 - dipoles = [] - current_id = [] - for val in currents.ab_cell_id.values: - cell_id = int(currents.ab_map[val]) - 1 - - for dipole in range(n_dipoles): - dipole_ids = currents.cells[cell_id, :] + 2 + dipole - - if ( - any(dipole_ids > (potentials.n_vertices - 1)) - or len(np.unique(parts[dipole_ids])) > 1 - ): - continue - - dipoles += [dipole_ids] - current_id += [val] - - potentials.cells = np.vstack(dipoles).astype("uint32") - potentials.ab_cell_id = np.hstack(current_id).astype("int32") - potentials.current_electrodes = currents - currents.potential_electrodes = potentials - split_dcip_survey(potentials, parts + 1, "DC Survey Line", workspace) - assert workspace.get_entity("DC Survey Line 4")[0] is not None + + with Workspace.create(path) as workspace: + potentials = generate_dc_survey(workspace, x_loc, y_loc) + line_id = potentials.get_data("line_ids")[0].values + + surveys = split_dcip_survey(potentials, line_id, workspace) + assert len(surveys) == len(np.unique(line_id)) def test_rectangular_block(): From bc842cf2df69a620c3c22b9a0794901a60d71882 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 07:55:35 -0800 Subject: [PATCH 08/17] Base class pseudo 3d --- .../pseudo_three_dimensions/driver.py | 128 +-------------- geoapps/inversion/electricals/driver.py | 146 ++++++++++++++++++ .../pseudo_three_dimensions/driver.py | 140 +---------------- 3 files changed, 153 insertions(+), 261 deletions(-) create mode 100644 geoapps/inversion/electricals/driver.py diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py index 5c6a59c18..5f8bac455 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/driver.py @@ -7,19 +7,6 @@ from __future__ import annotations -import sys -import uuid -from copy import deepcopy -from pathlib import Path - -import numpy as np -from geoh5py.data import Data -from geoh5py.objects import DrapeModel -from geoh5py.workspace import Workspace - -from geoapps.inversion.components.data import InversionData -from geoapps.inversion.components.topography import InversionTopography -from geoapps.inversion.components.windows import InversionWindow from geoapps.inversion.electricals.direct_current.pseudo_three_dimensions.constants import ( validations, ) @@ -29,119 +16,10 @@ from geoapps.inversion.electricals.direct_current.two_dimensions.params import ( DirectCurrent2DParams, ) -from geoapps.inversion.line_sweep.driver import LineSweepDriver -from geoapps.shared_utils.utils import get_locations, weighted_average -from geoapps.utils.models import get_drape_model -from geoapps.utils.surveys import extract_dcip_survey +from geoapps.inversion.electricals.driver import BasePseudo3DDriver -class DirectCurrentPseudo3DDriver(LineSweepDriver): +class DirectCurrentPseudo3DDriver(BasePseudo3DDriver): _params_class = DirectCurrentPseudo3DParams + _params_2d_class = DirectCurrent2DParams _validations = validations - - def __init__(self, params: DirectCurrentPseudo3DParams): # pylint: disable=W0235 - super().__init__(params) - if params.files_only: - sys.exit("Files written") - - def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID]: - xyz_in = get_locations(self.workspace, self.pseudo3d_params.mesh) - models = {"starting_model": self.pseudo3d_params.starting_model} - if not self.pseudo3d_params.forward_only: - models.update( - { - "reference_model": self.pseudo3d_params.reference_model, - "lower_bound": self.pseudo3d_params.lower_bound, - "upper_bound": self.pseudo3d_params.upper_bound, - } - ) - - xyz_out = mesh.centroids - model_uids = deepcopy(models) - 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}}) - model_uids[name] = model_object.uid - - return model_uids - - def write_files(self, lookup): - """Write ui.geoh5 and ui.json files for sweep trials.""" - - forward_only = self.pseudo3d_params.forward_only - ifile = DirectCurrent2DParams(forward_only=forward_only).input_file - - with self.workspace.open(mode="r+"): - self._window = InversionWindow(self.workspace, self.pseudo3d_params) - self._inversion_data = InversionData(self.workspace, self.pseudo3d_params) - self._inversion_topography = InversionTopography( - self.workspace, self.pseudo3d_params - ) - - for uid, trial in lookup.items(): - if trial["status"] != "pending": - continue - - filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" - with Workspace(filepath) as iter_workspace: - receiver_entity = extract_dcip_survey( - iter_workspace, - self.inversion_data.entity, - self.pseudo3d_params.line_object.values, - trial["line_id"], - ) - 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.pseudo3d_params.u_cell_size, - self.pseudo3d_params.v_cell_size, - ], - self.pseudo3d_params.depth_core, - [self.pseudo3d_params.horizontal_padding] * 2 - + [self.pseudo3d_params.vertical_padding, 1], - self.pseudo3d_params.expansion_factor, - )[0] - - model_uids = self.transfer_models(mesh) - - for key in ifile.data: - param = getattr(self.pseudo3d_params, key, None) - if key not in ["title", "inversion_type"]: - ifile.data[key] = param - - self.pseudo3d_params.topography_object.copy( - parent=iter_workspace, copy_children=True - ) - - ifile.data.update( - dict( - **{ - "geoh5": iter_workspace, - "mesh": mesh, - "data_object": receiver_entity, - "line_id": trial["line_id"], - "out_group": None, - }, - **model_uids, - ) - ) - - ifile.name = f"{uid}.ui.json" - ifile.path = self.working_directory # pylint: disable=E1101 - ifile.write_ui_json() - lookup[uid]["status"] = "written" - - _ = self.update_lookup(lookup) # pylint: disable=E1101 diff --git a/geoapps/inversion/electricals/driver.py b/geoapps/inversion/electricals/driver.py new file mode 100644 index 000000000..69c330023 --- /dev/null +++ b/geoapps/inversion/electricals/driver.py @@ -0,0 +1,146 @@ +# Copyright (c) 2024 Mira Geoscience Ltd. +# +# This file is part of geoapps. +# +# geoapps is distributed under the terms and conditions of the MIT License +# (see LICENSE file at the root of this source code package). +# +# This file is part of geoapps. +# +# geoapps 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 +from copy import deepcopy +from pathlib import Path + +import numpy as np +from geoh5py.data import Data +from geoh5py.objects import DrapeModel +from geoh5py.workspace import Workspace + +from geoapps.inversion.components.data import InversionData +from geoapps.inversion.components.topography import InversionTopography +from geoapps.inversion.components.windows import InversionWindow +from geoapps.inversion.line_sweep.driver import LineSweepDriver +from geoapps.inversion.params import BaseParams +from geoapps.shared_utils.utils import get_locations, weighted_average +from geoapps.utils.models import get_drape_model +from geoapps.utils.surveys import extract_dcip_survey + + +class BasePseudo3DDriver(LineSweepDriver): + _params_class: type(BaseParams) + _params_2d_class: type(BaseParams) + _validations: dict + _model_list: list[str] = [] + + def __init__(self, params): # pylint: disable=W0235 + super().__init__(params) + if params.files_only: + sys.exit("Files written") + + def transfer_models(self, mesh: DrapeModel) -> dict[str, uuid.UUID]: + xyz_in = get_locations(self.workspace, self.pseudo3d_params.mesh) + + models = {"starting_model": self.pseudo3d_params.starting_model} + + for model in self._model_list: + models[model] = getattr(self.pseudo3d_params, model) + + if not self.pseudo3d_params.forward_only: + for model in ["reference_model", "lower_bound", "upper_bound"]: + models[model] = getattr(self.pseudo3d_params, model) + + xyz_out = mesh.centroids + model_uids = deepcopy(models) + 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}}) + model_uids[name] = model_object.uid + + return model_uids + + def write_files(self, lookup): + """Write ui.geoh5 and ui.json files for sweep trials.""" + + forward_only = self.pseudo3d_params.forward_only + ifile = self._params_2d_class(forward_only=forward_only).input_file + + with self.workspace.open(mode="r+"): + self._window = InversionWindow(self.workspace, self.pseudo3d_params) + self._inversion_data = InversionData(self.workspace, self.pseudo3d_params) + self._inversion_topography = InversionTopography( + self.workspace, self.pseudo3d_params + ) + + for uid, trial in lookup.items(): + if trial["status"] != "pending": + continue + + filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" + with Workspace(filepath) as iter_workspace: + receiver_entity = extract_dcip_survey( + iter_workspace, + self.inversion_data.entity, + self.pseudo3d_params.line_object.values, + trial["line_id"], + ) + 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.pseudo3d_params.u_cell_size, + self.pseudo3d_params.v_cell_size, + ], + self.pseudo3d_params.depth_core, + [self.pseudo3d_params.horizontal_padding] * 2 + + [self.pseudo3d_params.vertical_padding, 1], + self.pseudo3d_params.expansion_factor, + )[0] + + model_uids = self.transfer_models(mesh) + + for key in ifile.data: + param = getattr(self.pseudo3d_params, key, None) + if key not in ["title", "inversion_type"]: + ifile.data[key] = param + + self.pseudo3d_params.topography_object.copy( + parent=iter_workspace, copy_children=True + ) + + ifile.data.update( + dict( + **{ + "geoh5": iter_workspace, + "mesh": mesh, + "data_object": receiver_entity, + "line_id": trial["line_id"], + "out_group": None, + }, + **model_uids, + ) + ) + + ifile.name = f"{uid}.ui.json" + ifile.path = self.working_directory # pylint: disable=E1101 + ifile.write_ui_json() + lookup[uid]["status"] = "written" + + _ = self.update_lookup(lookup) # pylint: disable=E1101 diff --git a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/driver.py b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/driver.py index b84942e77..d0c6806f8 100644 --- a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/driver.py +++ b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/driver.py @@ -7,17 +7,7 @@ from __future__ import annotations -import sys -from copy import deepcopy -from pathlib import Path - -import numpy as np -from geoh5py.data import Data -from geoh5py.workspace import Workspace - -from geoapps.inversion.components.data import InversionData -from geoapps.inversion.components.topography import InversionTopography -from geoapps.inversion.components.windows import InversionWindow +from geoapps.inversion.electricals.driver import BasePseudo3DDriver from geoapps.inversion.electricals.induced_polarization.pseudo_three_dimensions.constants import ( validations, ) @@ -27,132 +17,10 @@ from geoapps.inversion.electricals.induced_polarization.two_dimensions.params import ( InducedPolarization2DParams, ) -from geoapps.inversion.line_sweep.driver import LineSweepDriver -from geoapps.shared_utils.utils import get_locations, weighted_average -from geoapps.utils.models import get_drape_model -from geoapps.utils.surveys import extract_dcip_survey -class InducedPolarizationPseudo3DDriver(LineSweepDriver): +class InducedPolarizationPseudo3DDriver(BasePseudo3DDriver): _params_class = InducedPolarizationPseudo3DParams + _params_2d_class = InducedPolarization2DParams _validations = validations - - def __init__( - self, params: InducedPolarizationPseudo3DParams - ): # pylint: disable=W0235 - super().__init__(params) - if params.files_only: - sys.exit("Files written") - - def write_files(self, lookup: dict) -> None: - """ - Write ui.geoh5 and ui.json files for sweep trials. - - :param lookup: dictionary of trial hashes and trial - parameter values and status data. - """ - - forward_only = self.pseudo3d_params.forward_only - ifile = InducedPolarization2DParams(forward_only=forward_only).input_file - - with self.workspace.open(mode="r+"): - self._window = InversionWindow(self.workspace, self.pseudo3d_params) - self._inversion_data = InversionData(self.workspace, self.pseudo3d_params) - - self._inversion_topography = InversionTopography( - self.workspace, self.pseudo3d_params - ) - - xyz_in = get_locations(self.workspace, self.pseudo3d_params.mesh) - models = { - "starting_model": self.pseudo3d_params.starting_model, - "conductivity_model": self.pseudo3d_params.conductivity_model, - } - if not forward_only: - models.update( - { - "reference_model": self.pseudo3d_params.reference_model, - "lower_bound": self.pseudo3d_params.lower_bound, - "upper_bound": self.pseudo3d_params.upper_bound, - } - ) - - for uuid, trial in lookup.items(): - if trial["status"] != "pending": - continue - - filepath = Path(self.working_directory) / f"{uuid}.ui.geoh5" - with Workspace(filepath) as iter_workspace: - receiver_entity = extract_dcip_survey( - iter_workspace, - self.inversion_data.entity, - self.pseudo3d_params.line_object.values, - trial["line_id"], - ) - 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.pseudo3d_params.u_cell_size, - self.pseudo3d_params.v_cell_size, - ], - self.pseudo3d_params.depth_core, - [self.pseudo3d_params.horizontal_padding] * 2 - + [self.pseudo3d_params.vertical_padding, 1], - self.pseudo3d_params.expansion_factor, - )[0] - - iter_workspace.remove_entity(receiver_entity.current_electrodes) - iter_workspace.remove_entity(receiver_entity) - - xyz_out = mesh.centroids - model_uids = deepcopy(models) - 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}}) - model_uids[name] = model_object.uid - - for key in ifile.data: - param = getattr(self.pseudo3d_params, key, None) - if key not in ["title", "inversion_type"]: - ifile.data[key] = param - - self.pseudo3d_params.topography_object.copy( - parent=iter_workspace, copy_children=True - ) - self.pseudo3d_params.data_object.copy( - parent=iter_workspace, copy_children=True - ) - - ifile.data.update( - dict( - **{ - "geoh5": iter_workspace, - "mesh": mesh, - "line_id": trial["line_id"], - "out_group": None, - }, - **model_uids, - ) - ) - - ifile.name = f"{uuid}.ui.json" - ifile.path = self.working_directory # pylint: disable=E1101 - ifile.write_ui_json() - lookup[uuid]["status"] = "written" - - _ = self.update_lookup(lookup) # pylint: disable=E1101 + _model_list = ["conductivity_model"] From d128da2ea184b7b7337fd78ad7432f2b5237d9b9 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 13:23:11 -0800 Subject: [PATCH 09/17] Major re-work of logic --- geoapps/inversion/components/data.py | 64 ++++--------- .../factories/directives_factory.py | 11 +-- .../components/factories/entity_factory.py | 93 +++++++------------ .../components/factories/survey_factory.py | 22 ++--- geoapps/inversion/components/locations.py | 46 +-------- geoapps/inversion/components/topography.py | 5 - geoapps/inversion/driver.py | 2 +- .../pseudo_three_dimensions/constants.py | 4 - geoapps/inversion/electricals/driver.py | 12 ++- .../pseudo_three_dimensions/constants.py | 4 - geoapps/inversion/electricals/params.py | 12 +-- geoapps/inversion/line_sweep/driver.py | 19 +++- geoapps/inversion/params.py | 11 +++ geoapps/shared_utils/utils.py | 14 ++- geoapps/utils/surveys.py | 20 ++-- geoapps/utils/testing.py | 4 +- tests/run_tests/driver_dc_2d_test.py | 9 +- tests/run_tests/driver_ip_2d_test.py | 10 +- tests/run_tests/driver_ip_p3d_test.py | 7 +- 19 files changed, 132 insertions(+), 237 deletions(-) diff --git a/geoapps/inversion/components/data.py b/geoapps/inversion/components/data.py index 3d6725db7..fc34d3232 100644 --- a/geoapps/inversion/components/data.py +++ b/geoapps/inversion/components/data.py @@ -87,7 +87,6 @@ def __init__(self, workspace: Workspace, params: InversionBaseParams): self.radar: np.ndarray | None = None self.locations: np.ndarray | None = None self.mask: np.ndarray | None = None - self.global_map: np.ndarray | None = None self.indices: np.ndarray | None = None self.vector: bool | None = None self.n_blocks: int | None = None @@ -113,12 +112,16 @@ def _initialize(self) -> None: self.offset, self.radar = self.params.offset() self.locations = super().get_locations(self.params.data_object) - if self.angle is not None and self.angle != 0: - raise ValueError("Mesh is rotated.") - self.mask = np.ones(len(self.locations), dtype=bool) - if self.radar is not None: - if any(np.isnan(self.radar)): - self.mask[np.isnan(self.radar)] = False + if ( + getattr(self.params, "line_id", None) is not None + and getattr(self.params, "line_object", None) is not None + ): + self.mask = self.params.line_object.values == self.params.line_id + else: + self.mask = np.ones(len(self.locations), dtype=bool) + + if self.radar is not None and any(np.isnan(self.radar)): + self.mask[np.isnan(self.radar)] = False self.observed = self.filter(self.observed) self.radar = self.filter(self.radar) @@ -127,8 +130,8 @@ def _initialize(self) -> None: self.normalizations = self.get_normalizations() self.observed = self.normalize(self.observed) self.uncertainties = self.normalize(self.uncertainties, absolute=True) - self.locations = self.apply_transformations(self.locations) self.entity = self.write_entity() + self.params.data_object = self.entity self.locations = super().get_locations(self.entity) self.survey, self.local_index, _ = self.create_survey() @@ -163,25 +166,8 @@ def drape_locations(self, locations: np.ndarray) -> np.ndarray: def filter(self, a): """Remove vertices based on mask property.""" - if ( - self.params.inversion_type - in [ - "direct current pseudo 3d", - "direct current 3d", - "direct current 2d", - "induced polarization 3d", - "induced polarization 2d", - "induced polarization pseudo 3d", - ] - and self.indices is None - ): - ab_ind = np.where(np.any(self.mask[self.params.data_object.cells], axis=1))[ - 0 - ] - self.indices = ab_ind - if self.indices is None: - self.indices = np.where(self.mask) + self.indices = np.where(self.mask)[0] a = super().filter(a, mask=self.indices) @@ -251,11 +237,10 @@ def save_data(self, entity): else: for component in data: dnorm = data[component] / self.normalizations[None][component] - if "2d" in self.params.inversion_type: - dnorm = self._embed_2d(dnorm) data_dict[component] = entity.add_data( {f"{basename}_{component}": {"values": dnorm}} ) + if not self.params.forward_only: self._observed_data_types[component] = data_dict[ component @@ -265,20 +250,14 @@ def save_data(self, entity): / self.normalizations[None][component] ) uncerts[np.isinf(uncerts)] = np.nan - if "2d" in self.params.inversion_type: - uncerts = self._embed_2d(uncerts) + uncert_dict[component] = entity.add_data( {f"Uncertainties_{component}": {"values": uncerts}} ) if "direct current" in self.params.inversion_type: apparent_property = data[component].copy() - apparent_property[self.global_map] *= self.transformations[ - "apparent resistivity" - ] - - if "2d" in self.params.inversion_type: - apparent_property = self._embed_2d(apparent_property) + apparent_property *= self.transformations["apparent resistivity"] data_dict["apparent_resistivity"] = entity.add_data( { @@ -299,8 +278,7 @@ def apply_transformations(self, locations: np.ndarray): locations = self.displace(locations, self.offset) if self.radar is not None: locations = self.drape(locations, self.radar) - if self.is_rotated: - locations = super().rotate(locations) + return locations def displace(self, locs: np.ndarray, offset: np.ndarray) -> np.ndarray: @@ -503,12 +481,6 @@ def observed_data_types(self): """ return self._observed_data_types - def _embed_2d(self, data): - ind = np.ones_like(data, dtype=bool) - ind[self.global_map] = False - data[ind] = np.nan - return data - @staticmethod def check_tensor(channels): tensor_components = ["xx", "xy", "xz", "yx", "zx", "yy", "zz", "zy", "yz"] @@ -531,5 +503,7 @@ def update_params(self, data_dict, uncert_dict): setattr(self.params, f"{comp}_uncertainty", uncert_dict[comp]) if getattr(self.params, "line_object", None) is not None: - new_line = self.params.line_object.copy(parent=self.entity) + new_line = self.params.line_object.copy( + parent=self.entity, values=self.params.line_object.values[self.mask] + ) self.params.line_object = new_line diff --git a/geoapps/inversion/components/factories/directives_factory.py b/geoapps/inversion/components/factories/directives_factory.py index 79579f06e..8978c8b62 100644 --- a/geoapps/inversion/components/factories/directives_factory.py +++ b/geoapps/inversion/components/factories/directives_factory.py @@ -435,18 +435,9 @@ def assemble_data_keywords_dcip( "association": "CELL", } - if sorting is not None and "2d" not in self.factory_type: + if sorting is not None: kwargs["sorting"] = np.hstack(sorting) - if "2d" in self.factory_type: - - def transform_2d(x): - expanded_data = np.array([np.nan] * len(inversion_object.indices)) - expanded_data[inversion_object.global_map] = x[sorting] - return expanded_data - - kwargs["transforms"].insert(0, transform_2d) - if is_dc and name == "Apparent Resistivity": kwargs["transforms"].insert( 0, diff --git a/geoapps/inversion/components/factories/entity_factory.py b/geoapps/inversion/components/factories/entity_factory.py index 88333b9d2..4e76a075f 100644 --- a/geoapps/inversion/components/factories/entity_factory.py +++ b/geoapps/inversion/components/factories/entity_factory.py @@ -16,7 +16,7 @@ from geoapps.inversion.components.data import InversionData import numpy as np -from geoh5py.objects import Curve, Grid2D +from geoh5py.objects import CurrentElectrode, Curve, Grid2D, Points, PotentialElectrode from geoapps.inversion.components.factories.abstract_factory import AbstractFactory @@ -35,13 +35,9 @@ def factory_type(self): def concrete_object(self): """Returns a geoh5py object to be constructed by the build method.""" if "current" in self.factory_type or "polarization" in self.factory_type: - from geoh5py.objects import CurrentElectrode, PotentialElectrode - - return (PotentialElectrode, CurrentElectrode) + return PotentialElectrode, CurrentElectrode elif isinstance(self.params.data_object, Grid2D): - from geoh5py.objects import Points - return Points else: @@ -50,50 +46,18 @@ def concrete_object(self): def build(self, inversion_data: InversionData): """Constructs geoh5py object for provided inversion type.""" - if "current" in self.factory_type or "polarization" in self.factory_type: - entity = self._build_dcip(inversion_data) - else: - entity = self._build(inversion_data) + entity = self._build(inversion_data) return entity - def _build_dcip(self, inversion_data: InversionData): - PotentialElectrode, CurrentElectrode = self.concrete_object - workspace = inversion_data.workspace - - # Trim down receivers - rx_obj = self.params.data_object - rcv_ind = np.where(np.any(inversion_data.mask[rx_obj.cells], axis=1))[0] - rcv_locations, rcv_cells = EntityFactory._prune_from_indices(rx_obj, rcv_ind) - uni_src_ids, src_ids = np.unique( - rx_obj.ab_cell_id.values[rcv_ind], return_inverse=True - ) - ab_cell_id = np.arange(1, uni_src_ids.shape[0] + 1)[src_ids] - entity = PotentialElectrode.create( - workspace, - name="Data", - parent=self.params.out_group, - vertices=inversion_data.apply_transformations(rcv_locations), - cells=rcv_cells, - ) - entity.ab_cell_id = ab_cell_id - # Trim down sources - tx_obj = rx_obj.current_electrodes - src_ind = np.hstack( - [np.where(tx_obj.ab_cell_id.values == ind)[0] for ind in uni_src_ids] - ) - src_locations, src_cells = EntityFactory._prune_from_indices(tx_obj, src_ind) - new_currents = CurrentElectrode.create( - workspace, - name="Data (currents)", - parent=self.params.out_group, - vertices=inversion_data.apply_transformations(src_locations), - cells=src_cells, - ) - new_currents.add_default_ab_cell_id() - entity.current_electrodes = new_currents - - return entity + # def _build_dcip(self, inversion_data: InversionData): + # entity = extract_dcip_survey( + # inversion_data.workspace, + # self.params.data_object, + # inversion_data.indices + # ) + # + # return entity def _build(self, inversion_data: InversionData): if isinstance(self.params.data_object, Grid2D): @@ -102,11 +66,27 @@ def _build(self, inversion_data: InversionData): ) else: - entity = self.params.data_object.copy( - parent=self.params.out_group, - copy_children=False, - vertices=inversion_data.locations, - ) + kwargs = { + "parent": self.params.out_group, + "copy_children": False, + } + + if np.any(~inversion_data.mask): + if isinstance(self.params.data_object, PotentialElectrode): + active_poles = np.zeros( + self.params.data_object.n_vertices, dtype=bool + ) + active_poles[ + self.params.data_object.cells[inversion_data.mask, :].ravel() + ] = True + kwargs.update( + {"mask": active_poles, "cell_mask": inversion_data.mask} + ) + else: + kwargs.update({"mask": inversion_data.mask}) + + entity = self.params.data_object.copy(**kwargs) + entity.vertices = inversion_data.apply_transformations(entity.vertices) if getattr(entity, "transmitters", None) is not None: entity.transmitters.vertices = inversion_data.apply_transformations( @@ -116,11 +96,10 @@ def _build(self, inversion_data: InversionData): if tx_freq: tx_freq[0].copy(parent=entity.transmitters) - if np.any(~inversion_data.mask): - entity.remove_vertices(~inversion_data.mask) - - if getattr(entity, "transmitters", None) is not None: - entity.transmitters.remove_vertices(~inversion_data.mask) + if getattr(entity, "current_electrodes", None) is not None: + entity.current_electrodes.vertices = inversion_data.apply_transformations( + entity.current_electrodes.vertices + ) return entity diff --git a/geoapps/inversion/components/factories/survey_factory.py b/geoapps/inversion/components/factories/survey_factory.py index 8d91ec127..0efa6fab6 100644 --- a/geoapps/inversion/components/factories/survey_factory.py +++ b/geoapps/inversion/components/factories/survey_factory.py @@ -19,8 +19,6 @@ import SimPEG.electromagnetics.time_domain as tdem from scipy.interpolate import interp1d -from geoapps.utils.surveys import extract_dcip_survey - from .receiver_factory import ReceiversFactory from .simpeg_factory import SimPEGFactory from .source_factory import SourcesFactory @@ -227,12 +225,7 @@ def _add_data(self, survey, data, local_index, channel): survey.std = uncertainty_vec else: - index_map = ( - data.global_map[local_index] - if data.global_map is not None - else local_index - ) - local_data = {k: v[index_map] for k, v in data.observed.items()} + local_data = {k: v[local_index] for k, v in data.observed.items()} local_uncertainties = { k: v[local_index] for k, v in data.uncertainties.items() } @@ -274,12 +267,13 @@ def _dcip_arguments(self, data=None, local_index=None): receiver_entity = data.entity if "2d" in self.factory_type: - receiver_entity = extract_dcip_survey( - self.params.geoh5, - receiver_entity, - self.params.line_object.values, - self.params.line_id, - ) + # receiver_entity = extract_dcip_survey( + # self.params.geoh5, + # receiver_entity, + # self.params.line_object.values, + # self.params.line_id, + # ) + self.local_index = np.arange(receiver_entity.n_cells) source_ids, order = np.unique( receiver_entity.ab_cell_id.values[self.local_index], return_index=True diff --git a/geoapps/inversion/components/locations.py b/geoapps/inversion/components/locations.py index 02e75387f..e50c0a34d 100644 --- a/geoapps/inversion/components/locations.py +++ b/geoapps/inversion/components/locations.py @@ -22,7 +22,6 @@ from scipy.spatial import cKDTree from geoapps.shared_utils.utils import get_locations as get_locs -from geoapps.shared_utils.utils import rotate_xyz class InversionLocations: @@ -33,12 +32,6 @@ class InversionLocations: ---------- mask : Mask that stores cumulative filtering actions. - origin : - Rotation origin. - angle : - Rotation angle. - is_rotated : - True if locations have been rotated. locations : xyz locations. @@ -48,8 +41,6 @@ class InversionLocations: Returns locations of data object centroids or vertices. filter() : Apply accumulated self.mask to array, or dict of arrays. - rotate() : - Un-rotate data using origin and angle assigned to inversion mesh. """ @@ -60,17 +51,8 @@ def __init__(self, workspace: Workspace, params: InversionBaseParams): """ self.workspace = workspace self._params: InversionBaseParams = params - self.mask: np.ndarray = None - self.origin: list[float] = None - self.angle: float = None - self.is_rotated: bool = False - self.locations: np.ndarray = None - - if params.mesh is not None: - if hasattr(params.mesh, "rotation"): - self.origin = np.asarray(params.mesh.origin.tolist()) - self.angle = -1 * params.mesh.rotation - self.is_rotated = True if np.abs(self.angle) != 0 else False + self.mask: np.ndarray | None = None + self.locations: np.ndarray | None = None @property def mask(self): @@ -90,14 +72,6 @@ def mask(self, v): def create_entity(self, name, locs: np.ndarray, geoh5_object=Points): """Create Data group and Points object with observed data.""" - - if self.is_rotated: - locs[:, :2] = rotate_xyz( - locs[:, :2], - self.origin, - -1 * self.angle, - ) - entity = geoh5_object.create( self.workspace, name=name, @@ -168,22 +142,6 @@ def filter(self, a: dict[str, np.ndarray] | np.ndarray, mask=None): else: return a[mask] - def rotate(self, locs: np.ndarray) -> np.ndarray: - """ - Rotate data using origin and angle assigned to inversion mesh. - - Since rotation attribute is stored with a negative sign the applied - rotation will restore locations to an East-West, North-South orientation. - - :param locs: Array of xyz locations. - """ - - if locs is None: - return None - - xy = rotate_xyz(locs[:, :2], self.origin, self.angle) - return np.c_[xy, locs[:, 2]] - def set_z_from_topo(self, locs: np.ndarray): """interpolate locations z data from topography.""" diff --git a/geoapps/inversion/components/topography.py b/geoapps/inversion/components/topography.py index f1f24c82e..bfa0e52af 100644 --- a/geoapps/inversion/components/topography.py +++ b/geoapps/inversion/components/topography.py @@ -67,14 +67,9 @@ def _initialize(self): self.mask = filter_xy( self.locations[:, 0], self.locations[:, 1], - angle=self.angle, ) - self.locations = super().filter(self.locations) - if self.is_rotated: - self.locations = super().rotate(self.locations) - def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: """ Return mask that restricts models to set of earth cells. diff --git a/geoapps/inversion/driver.py b/geoapps/inversion/driver.py index a5d3d8996..2416a9311 100644 --- a/geoapps/inversion/driver.py +++ b/geoapps/inversion/driver.py @@ -402,7 +402,7 @@ def get_regularization(self): def get_tiles(self): if "2d" in self.params.inversion_type: - tiles = [self.inversion_data.indices] + tiles = [np.arange(len(self.inversion_data.indices))] else: locations = self.inversion_data.locations diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py index 8b26b059d..c9a94e140 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py @@ -27,7 +27,6 @@ "topography": None, "data_object": None, "line_object": None, - "line_id": 1, "z_from_topo": False, "receivers_radar_drape": None, "receivers_offset_z": 0.0, @@ -104,7 +103,6 @@ "topography": None, "data_object": None, "line_object": None, - "line_id": 1, "z_from_topo": False, "receivers_radar_drape": None, "receivers_offset_z": 0.0, @@ -149,7 +147,6 @@ "parent": "data_object", "value": None, }, - "line_id": 1, "data_object": { "main": True, "group": "Data", @@ -318,7 +315,6 @@ "potential_channel": UUID("{502e7256-aafa-4016-969f-5cc3a4f27315}"), "potential_uncertainty": UUID("{62746129-3d82-427e-a84c-78cded00c0bc}"), "line_object": UUID("{d400e8f1-8460-4609-b852-b3b93f945770}"), - "line_id": 1, "mesh": UUID("{da109284-aa8c-4824-a647-29951109b058}"), "starting_model": 1e-1, "reference_model": 1e-1, diff --git a/geoapps/inversion/electricals/driver.py b/geoapps/inversion/electricals/driver.py index 69c330023..f93a87570 100644 --- a/geoapps/inversion/electricals/driver.py +++ b/geoapps/inversion/electricals/driver.py @@ -89,11 +89,15 @@ def write_files(self, lookup): filepath = Path(self.working_directory) / f"{uid}.ui.geoh5" with Workspace(filepath) as iter_workspace: + cell_mask: np.ndarray = ( + self.pseudo3d_params.line_object.values == trial["line_id"] + ) + + if not np.any(cell_mask): + continue + receiver_entity = extract_dcip_survey( - iter_workspace, - self.inversion_data.entity, - self.pseudo3d_params.line_object.values, - trial["line_id"], + iter_workspace, self.inversion_data.entity, cell_mask ) current_entity = receiver_entity.current_electrodes receiver_locs = np.vstack( diff --git a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py index ecbf23556..6c0e74141 100644 --- a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py @@ -27,7 +27,6 @@ "topography": None, "data_object": None, "line_object": None, - "line_id": 1, "z_from_topo": False, "receivers_radar_drape": None, "receivers_offset_z": 0.0, @@ -105,7 +104,6 @@ "topography": None, "data_object": None, "line_object": None, - "line_id": 1, "z_from_topo": False, "receivers_radar_drape": None, "receivers_offset_z": 0.0, @@ -151,7 +149,6 @@ "parent": "data_object", "value": None, }, - "line_id": 1, "data_object": { "main": True, "group": "Data", @@ -331,7 +328,6 @@ "chargeability_channel": UUID("{162320e6-2b80-4877-9ec1-a8f5b6a13673}"), "chargeability_uncertainty": 0.001, "line_object": UUID("{d400e8f1-8460-4609-b852-b3b93f945770}"), - "line_id": 1, "mesh": UUID("{da109284-aa8c-4824-a647-29951109b058}"), "starting_model": 1e-4, "conductivity_model": 0.1, diff --git a/geoapps/inversion/electricals/params.py b/geoapps/inversion/electricals/params.py index e95ba8ff4..9c41e55c6 100644 --- a/geoapps/inversion/electricals/params.py +++ b/geoapps/inversion/electricals/params.py @@ -11,7 +11,6 @@ from geoapps.inversion import InversionBaseParams from geoapps.utils.models import get_drape_model -from geoapps.utils.surveys import extract_dcip_survey class Base2DParams(InversionBaseParams): @@ -114,17 +113,10 @@ def expansion_factor(self, value): @property def mesh(self): if self._mesh is None and self.geoh5 is not None: - receiver_entity = extract_dcip_survey( - self.geoh5, - self.data_object, - self.line_object.values, - self.line_id, - ) - current_entity = receiver_entity.current_electrodes + current_entity = self.data_object.current_electrodes receiver_locs = np.vstack( - [receiver_entity.vertices, current_entity.vertices] + [self.data_object.vertices, current_entity.vertices] ) - self._mesh = get_drape_model( self.geoh5, "Models", diff --git a/geoapps/inversion/line_sweep/driver.py b/geoapps/inversion/line_sweep/driver.py index f3b925c72..175eced17 100644 --- a/geoapps/inversion/line_sweep/driver.py +++ b/geoapps/inversion/line_sweep/driver.py @@ -14,6 +14,7 @@ import numpy as np from geoh5py.data import FilenameData from geoh5py.groups import SimPEGGroup +from geoh5py.objects import DrapeModel, PotentialElectrode from geoh5py.shared.utils import fetch_active_workspace from geoh5py.ui_json import InputFile from geoh5py.workspace import Workspace @@ -121,7 +122,14 @@ def collect_results(self): drape_models = [] for line in np.unique(line_ids): with Workspace(f"{path / files[line]}.ui.geoh5") as ws: - survey = ws.get_entity("Data")[0] + out_group = [ + group for group in ws.groups if isinstance(group, SimPEGGroup) + ][0] + survey = [ + child + for child in out_group.children + if isinstance(child, PotentialElectrode) + ][0] line_data = survey.get_entity(self.pseudo3d_params.line_object.name) if not line_data: @@ -129,8 +137,11 @@ def collect_results(self): line_indices = line_ids == line data = self.collect_line_data(survey, line_indices, data) - - mesh = ws.get_entity("Models")[0] + mesh = [ + child + for child in out_group.children + if isinstance(child, DrapeModel) + ][0] filedata = [ k for k in mesh.parent.children if isinstance(k, FilenameData) ] @@ -141,8 +152,8 @@ def collect_results(self): ) for fdat in filedata: fdat.copy(parent=local_simpeg_group) + mesh = mesh.copy(parent=local_simpeg_group) - mesh.name = "models" drape_models.append(mesh) self.pseudo3d_params.data_object.add_data(data) diff --git a/geoapps/inversion/params.py b/geoapps/inversion/params.py index 1b795a671..034c19537 100644 --- a/geoapps/inversion/params.py +++ b/geoapps/inversion/params.py @@ -13,6 +13,7 @@ import numpy as np from geoh5py.data import NumericData from geoh5py.groups import SimPEGGroup +from geoh5py.objects import Octree from geoh5py.shared.utils import fetch_active_workspace from geoh5py.ui_json import InputFile @@ -346,6 +347,16 @@ def mesh(self): @mesh.setter def mesh(self, val): self.setter_validator("mesh", val, fun=self._uuid_promoter) + + if ( + isinstance(self._mesh, Octree) + and self._mesh.rotation is not None + and self._mesh.rotation != 0.0 + ): + raise ValueError( + "Rotated meshes are not supported. Please use a mesh with an angle of 0.0." + ) + self.update_group_options() @property diff --git a/geoapps/shared_utils/utils.py b/geoapps/shared_utils/utils.py index fb22a4e56..1519d17a2 100644 --- a/geoapps/shared_utils/utils.py +++ b/geoapps/shared_utils/utils.py @@ -14,6 +14,7 @@ import numpy as np from discretize import TensorMesh, TreeMesh from geoh5py.objects import Curve, DrapeModel +from geoh5py.objects.surveys.direct_current import BaseElectrode from geoh5py.shared import Entity from geoh5py.workspace import Workspace from scipy.interpolate import interp1d @@ -53,7 +54,18 @@ def get_locations(workspace: Workspace, entity: UUID | Entity): if hasattr(entity, "centroids"): locations = entity.centroids elif hasattr(entity, "vertices"): - locations = entity.vertices + if isinstance(entity, BaseElectrode): + potentials = entity.potential_electrodes + locations = np.mean( + [ + potentials.vertices[potentials.cells[:, 0], :], + potentials.vertices[potentials.cells[:, 1], :], + ], + axis=0, + ) + else: + locations = entity.vertices + elif getattr(entity, "parent", None) is not None and entity.parent is not None: locations = get_locations(workspace, entity.parent) diff --git a/geoapps/utils/surveys.py b/geoapps/utils/surveys.py index 98bb8349f..1c403da25 100644 --- a/geoapps/utils/surveys.py +++ b/geoapps/utils/surveys.py @@ -188,24 +188,15 @@ def slice_and_map(obj: np.ndarray, slicer: np.ndarray | Callable): def extract_dcip_survey( - workspace: Workspace, - survey: PotentialElectrode, - lines: np.ndarray, - line_id: int, + 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: lines: Line indexer for survey. - :param: line_id: Index of line to extract data from. + :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. """ - cell_mask = lines == line_id - - if not np.any(cell_mask): - raise ValueError(f"Line '{line_id}' not found in survey.") - 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) @@ -234,7 +225,8 @@ def split_dcip_survey( with ws.open(mode="r+") as ws: line_surveys = [] for line_id in np.unique(lines): - line_survey = extract_dcip_survey(ws, survey, lines, line_id) + cell_mask = lines == line_id + line_survey = extract_dcip_survey(ws, survey, cell_mask) line_surveys.append(line_survey) return line_surveys diff --git a/geoapps/utils/testing.py b/geoapps/utils/testing.py index 553284cd8..40a8cdb4e 100644 --- a/geoapps/utils/testing.py +++ b/geoapps/utils/testing.py @@ -396,11 +396,11 @@ def topo_drape(x, y): # Create a mesh if "2d" in inversion_type: - lines = survey.get_entity("line_id")[0].values + lines = survey.get_entity("line_ids")[0].values entity, mesh, _ = get_drape_model( # pylint: disable=W0632 geoh5, "Models", - survey.vertices[survey.cells[lines == 2, :], :], + survey.vertices[np.unique(survey.cells[lines == 101, :]), :], [cell_size[0], cell_size[2]], 100.0, [padding_distance] * 2 + [padding_distance] * 2, diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index 574b0eb90..b5c988928 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -19,7 +19,6 @@ DirectCurrent2DParams, ) from geoapps.shared_utils.utils import get_inversion_output -from geoapps.utils.surveys import survey_lines from geoapps.utils.testing import check_target, setup_inversion_workspace # To test the full run and validate the inversion. @@ -52,7 +51,6 @@ def test_dc_2d_fwr_run( drape_height=0.0, flatten=False, ) - _ = survey_lines(survey, [-100, -100], save="line_ids") params = DirectCurrent2DParams( forward_only=True, geoh5=geoh5, @@ -67,7 +65,7 @@ def test_dc_2d_fwr_run( data_object=survey.uid, starting_model=model.uid, line_object=geoh5.get_entity("line_ids")[0].uid, - line_id=2, + line_id=101, ) params.workpath = tmp_path fwr_driver = DirectCurrent2DDriver(params) @@ -82,7 +80,6 @@ def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): with Workspace(workpath) as geoh5: potential = geoh5.get_entity("Iteration_0_dc")[0] topography = geoh5.get_entity("topography")[0] - _ = survey_lines(potential.parent, [-100, 100], save="line_IDs") # Run the inverse np.random.seed(0) @@ -98,8 +95,8 @@ def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): data_object=potential.parent.uid, potential_channel=potential.uid, potential_uncertainty=1e-3, - line_object=geoh5.get_entity("line_IDs")[0].uid, - line_id=2, + line_object=geoh5.get_entity("line_ids")[0].uid, + line_id=101, starting_model=1e-2, reference_model=1e-2, s_norm=0.0, diff --git a/tests/run_tests/driver_ip_2d_test.py b/tests/run_tests/driver_ip_2d_test.py index 5da4c129e..9810fcba5 100644 --- a/tests/run_tests/driver_ip_2d_test.py +++ b/tests/run_tests/driver_ip_2d_test.py @@ -19,7 +19,6 @@ InducedPolarization2DDriver, ) from geoapps.shared_utils.utils import get_inversion_output -from geoapps.utils.surveys import survey_lines from geoapps.utils.testing import check_target, setup_inversion_workspace # To test the full run and validate the inversion. @@ -52,7 +51,6 @@ def test_ip_2d_fwr_run( flatten=False, drape_height=0.0, ) - _ = survey_lines(survey, [-100, -100], save="line_ids") params = InducedPolarization2DParams( forward_only=True, geoh5=geoh5, @@ -63,7 +61,7 @@ def test_ip_2d_fwr_run( starting_model=model.uid, conductivity_model=1e-2, line_object=geoh5.get_entity("line_ids")[0].uid, - line_id=2, + line_id=101, ) params.workpath = tmp_path fwr_driver = InducedPolarization2DDriver(params) @@ -83,8 +81,6 @@ def test_ip_2d_run( chargeability = geoh5.get_entity("Iteration_0_ip")[0] mesh = geoh5.get_entity("Models")[0] topography = geoh5.get_entity("topography")[0] - _ = survey_lines(chargeability.parent, [-100, 100], save="line_IDs") - # Run the inverse np.random.seed(0) params = InducedPolarization2DParams( @@ -94,8 +90,8 @@ def test_ip_2d_run( data_object=chargeability.parent.uid, chargeability_channel=chargeability.uid, chargeability_uncertainty=2e-4, - line_object=geoh5.get_entity("line_IDs")[0].uid, - line_id=2, + line_object=geoh5.get_entity("line_ids")[0].uid, + line_id=101, starting_model=1e-6, reference_model=1e-6, conductivity_model=1e-2, diff --git a/tests/run_tests/driver_ip_p3d_test.py b/tests/run_tests/driver_ip_p3d_test.py index ba413fb51..bfec0d336 100644 --- a/tests/run_tests/driver_ip_p3d_test.py +++ b/tests/run_tests/driver_ip_p3d_test.py @@ -21,7 +21,6 @@ InducedPolarizationPseudo3DParams, ) from geoapps.shared_utils.utils import get_inversion_output -from geoapps.utils.surveys import survey_lines from geoapps.utils.testing import check_target, setup_inversion_workspace # To test the full run and validate the inversion. @@ -55,7 +54,6 @@ def test_ip_p3d_fwr_run( flatten=False, ) - _ = survey_lines(survey, [-100, -100], save="line_ids") params = InducedPolarizationPseudo3DParams( forward_only=True, geoh5=geoh5, @@ -92,7 +90,6 @@ def test_ip_p3d_run( out_group = geoh5.get_entity("Line 1")[0].parent mesh = out_group.get_entity("mesh")[0] # Finds the octree mesh topography = geoh5.get_entity("topography")[0] - _ = survey_lines(chargeability.parent, [-100, 100], save="line_IDs") # Run the inverse np.random.seed(0) @@ -108,7 +105,7 @@ def test_ip_p3d_run( data_object=chargeability.parent.uid, chargeability_channel=chargeability.uid, chargeability_uncertainty=2e-4, - line_object=geoh5.get_entity("line_IDs")[0].uid, + line_object=geoh5.get_entity("line_ids")[0].uid, conductivity_model=1e-2, starting_model=1e-6, reference_model=1e-6, @@ -134,7 +131,7 @@ def test_ip_p3d_run( basepath = workpath.parent with open(basepath / "lookup.json", encoding="utf8") as f: lookup = json.load(f) - middle_line_id = [k for k, v in lookup.items() if v["line_id"] == 2][0] + middle_line_id = [k for k, v in lookup.items() if v["line_id"] == 101][0] with Workspace(basepath / f"{middle_line_id}.ui.geoh5", mode="r") as workspace: middle_inversion_group = [ From 9f6033e192470aa14b94e33e7cfffc42969a9d5b Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 13:50:04 -0800 Subject: [PATCH 10/17] Fix tests --- geoapps/utils/surveys.py | 4 ++++ tests/run_tests/driver_ip_p3d_test.py | 4 ++-- tests/utils_test.py | 8 ++++---- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/geoapps/utils/surveys.py b/geoapps/utils/surveys.py index 1c403da25..dff4a4f8f 100644 --- a/geoapps/utils/surveys.py +++ b/geoapps/utils/surveys.py @@ -197,6 +197,10 @@ def extract_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) diff --git a/tests/run_tests/driver_ip_p3d_test.py b/tests/run_tests/driver_ip_p3d_test.py index bfec0d336..67b24f0fd 100644 --- a/tests/run_tests/driver_ip_p3d_test.py +++ b/tests/run_tests/driver_ip_p3d_test.py @@ -28,8 +28,8 @@ target_run = { "data_norm": 0.08768, - "phi_d": 7267, - "phi_m": 0.1079, + "phi_d": 8239, + "phi_m": 0.1178, } np.random.seed(0) diff --git a/tests/utils_test.py b/tests/utils_test.py index 30da1a766..24988391b 100644 --- a/tests/utils_test.py +++ b/tests/utils_test.py @@ -268,14 +268,14 @@ def test_extract_dcip_survey(tmp_path: Path): line_id = potentials.get_data("line_ids")[0].values - with pytest.raises(ValueError, match="Line '3' not found in survey."): - extract_dcip_survey(workspace, potentials, line_id, 3) + surveys = extract_dcip_survey(workspace, potentials, line_id == 3) - surveys = extract_dcip_survey(workspace, potentials, line_id, 6) + with pytest.raises(ValueError, match="No cells found in the mask."): + extract_dcip_survey(workspace, potentials, line_id == 6) line_field = surveys.get_entity("line_ids") assert line_field - assert np.all(line_field[0].values == 6) + assert np.all(line_field[0].values == 3) def test_split_dcip_survey(tmp_path: Path): From 31568b38df7b5d6be7f51ef18466c1cc798a1ed2 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 13:55:56 -0800 Subject: [PATCH 11/17] Update test target --- tests/run_tests/driver_dc_2d_test.py | 6 +++--- tests/run_tests/driver_dc_p3d_test.py | 2 +- tests/run_tests/driver_ip_2d_test.py | 6 +++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index b5c988928..7097dc14d 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -25,9 +25,9 @@ # Move this file out of the test directory and run. target_run = { - "data_norm": 0.59762, - "phi_d": 1425, - "phi_m": 7.851, + "data_norm": 0.59563, + "phi_d": 1400, + "phi_m": 8.004, } np.random.seed(0) diff --git a/tests/run_tests/driver_dc_p3d_test.py b/tests/run_tests/driver_dc_p3d_test.py index 7cba3fcc7..48be340eb 100644 --- a/tests/run_tests/driver_dc_p3d_test.py +++ b/tests/run_tests/driver_dc_p3d_test.py @@ -26,7 +26,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.1003794750530917, "phi_d": 4086, "phi_m": 0.8238} +target_run = {"data_norm": 1.099, "phi_d": 4150, "phi_m": 0.7511} np.random.seed(0) diff --git a/tests/run_tests/driver_ip_2d_test.py b/tests/run_tests/driver_ip_2d_test.py index 9810fcba5..550525d6a 100644 --- a/tests/run_tests/driver_ip_2d_test.py +++ b/tests/run_tests/driver_ip_2d_test.py @@ -25,9 +25,9 @@ # Move this file out of the test directory and run. target_run = { - "data_norm": 0.091298, - "phi_d": 9581, - "phi_m": 0.08327, + "data_norm": 0.09141, + "phi_d": 9934, + "phi_m": 0.08341, } np.random.seed(0) From fbde70f5799d224f32bda8176d928dfbd7f8aca7 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 14:03:26 -0800 Subject: [PATCH 12/17] Clean out commented code --- geoapps/inversion/components/factories/entity_factory.py | 9 --------- geoapps/inversion/components/factories/survey_factory.py | 6 ------ 2 files changed, 15 deletions(-) diff --git a/geoapps/inversion/components/factories/entity_factory.py b/geoapps/inversion/components/factories/entity_factory.py index 4e76a075f..0881831c2 100644 --- a/geoapps/inversion/components/factories/entity_factory.py +++ b/geoapps/inversion/components/factories/entity_factory.py @@ -50,15 +50,6 @@ def build(self, inversion_data: InversionData): return entity - # def _build_dcip(self, inversion_data: InversionData): - # entity = extract_dcip_survey( - # inversion_data.workspace, - # self.params.data_object, - # inversion_data.indices - # ) - # - # return entity - def _build(self, inversion_data: InversionData): if isinstance(self.params.data_object, Grid2D): entity = inversion_data.create_entity( diff --git a/geoapps/inversion/components/factories/survey_factory.py b/geoapps/inversion/components/factories/survey_factory.py index 0efa6fab6..21587cec7 100644 --- a/geoapps/inversion/components/factories/survey_factory.py +++ b/geoapps/inversion/components/factories/survey_factory.py @@ -267,12 +267,6 @@ def _dcip_arguments(self, data=None, local_index=None): receiver_entity = data.entity if "2d" in self.factory_type: - # receiver_entity = extract_dcip_survey( - # self.params.geoh5, - # receiver_entity, - # self.params.line_object.values, - # self.params.line_id, - # ) self.local_index = np.arange(receiver_entity.n_cells) source_ids, order = np.unique( From f83a3974fce27b9a2a309f56e7e9902cb8bdfaf3 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 15:12:38 -0800 Subject: [PATCH 13/17] Fix for app tests --- .../inversion/components/factories/receiver_factory.py | 7 ------- .../inversion/components/factories/survey_factory.py | 2 +- geoapps/inversion/components/preprocessing.py | 5 +++++ geoapps/inversion/driver.py | 10 ---------- tests/locations_test.py | 9 --------- 5 files changed, 6 insertions(+), 27 deletions(-) diff --git a/geoapps/inversion/components/factories/receiver_factory.py b/geoapps/inversion/components/factories/receiver_factory.py index da440c481..6d99c28a9 100644 --- a/geoapps/inversion/components/factories/receiver_factory.py +++ b/geoapps/inversion/components/factories/receiver_factory.py @@ -86,13 +86,6 @@ def assemble_arguments( args = [] - if getattr(self.params.mesh, "rotation", None): - locations = rotate_xyz( - locations, - self.params.mesh.origin.tolist(), - -1 * self.params.mesh.rotation[0], - ) - if ( "direct current" in self.factory_type or "induced polarization" in self.factory_type diff --git a/geoapps/inversion/components/factories/survey_factory.py b/geoapps/inversion/components/factories/survey_factory.py index 21587cec7..d2be08c24 100644 --- a/geoapps/inversion/components/factories/survey_factory.py +++ b/geoapps/inversion/components/factories/survey_factory.py @@ -278,7 +278,7 @@ def _dcip_arguments(self, data=None, local_index=None): receiver_locations = data.drape_locations(receiver_entity.vertices) source_locations = data.drape_locations(currents.vertices) else: - receiver_locations = data.locations + receiver_locations = receiver_entity.vertices source_locations = currents.vertices # TODO hook up tile_spatial to handle local_index handling diff --git a/geoapps/inversion/components/preprocessing.py b/geoapps/inversion/components/preprocessing.py index 8b23a9e02..be00ecbd5 100644 --- a/geoapps/inversion/components/preprocessing.py +++ b/geoapps/inversion/components/preprocessing.py @@ -97,6 +97,11 @@ def window_data( distance=resolution, ) + if isinstance(data_object, PotentialElectrode): + vert_mask = np.zeros(data_object.n_vertices, dtype=bool) + vert_mask[data_object.cells[mask, :].ravel()] = True + mask = vert_mask + new_data_object = data_object.copy( parent=data_object.workspace, copy_children=True, diff --git a/geoapps/inversion/driver.py b/geoapps/inversion/driver.py index 2416a9311..65f3b3c3c 100644 --- a/geoapps/inversion/driver.py +++ b/geoapps/inversion/driver.py @@ -406,16 +406,6 @@ def get_tiles(self): else: locations = self.inversion_data.locations - # Use mid-point between M-N electrodes - if self.params.inversion_type in [ - "direct current 3d", - "induced polarization 3d", - ]: - cells = self.inversion_data.entity.cells - locations = ( - locations[cells[:, 0], :] + locations[cells[:, 1], :] - ) / 2.0 - tiles = tile_locations( locations, self.params.tile_spatial, diff --git a/tests/locations_test.py b/tests/locations_test.py index 802c260cb..9adc6c623 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -81,15 +81,6 @@ def test_filter(tmp_path: Path): assert np.all(filtered_data["key"] == [2, 3, 4]) -def test_rotate(tmp_path: Path): - # Basic test since rotate_xy already tested - ws, params = setup_params(tmp_path) - locations = InversionLocations(ws, params) - test_locs = np.array([[1.0, 2.0, 3.0], [3.0, 4.0, 5.0], [6.0, 7.0, 8.0]]) - locs_rot = locations.rotate(test_locs) - assert locs_rot.shape == test_locs.shape - - def test_z_from_topo(tmp_path: Path): ws, params = setup_params(tmp_path) locations = InversionLocations(ws, params) From ca0dd64a7e2a52323631c632f6c3ec25a9e55b2f Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 16:02:36 -0800 Subject: [PATCH 14/17] Adjust test --- tests/run_tests/driver_dc_test.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/run_tests/driver_dc_test.py b/tests/run_tests/driver_dc_test.py index 8fbebb055..bc3cd047d 100644 --- a/tests/run_tests/driver_dc_test.py +++ b/tests/run_tests/driver_dc_test.py @@ -24,7 +24,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.152105803389558, "phi_d": 31.56, "phi_m": 171.5} +target_run = {"data_norm": 0.15258, "phi_d": 31.85, "phi_m": 122.7} np.random.seed(0) @@ -59,6 +59,8 @@ def test_dc_3d_fwr_run( survey.ab_cell_id = tx_id survey.cells = cells + geoh5.close() + params = DirectCurrent3DParams( forward_only=True, geoh5=geoh5, From 290b6cab48d3797524be1799d2a4e4c119318cc3 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Wed, 7 Feb 2024 19:57:01 -0800 Subject: [PATCH 15/17] fix IP test --- tests/run_tests/driver_ip_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/run_tests/driver_ip_test.py b/tests/run_tests/driver_ip_test.py index 09993d1e6..811366220 100644 --- a/tests/run_tests/driver_ip_test.py +++ b/tests/run_tests/driver_ip_test.py @@ -24,7 +24,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.00838140484272608, "phi_d": 1.88, "phi_m": 0.4663} +target_run = {"data_norm": 0.008494, "phi_d": 1.734, "phi_m": 0.3202} np.random.seed(0) From cb42ba9466ff8382adfac4312344702e8d181e94 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 13 Feb 2024 10:29:17 -0800 Subject: [PATCH 16/17] Reduce code duplication with refactor of params classes. Add check on line object for association --- .../pseudo_three_dimensions/constants.py | 2 - .../pseudo_three_dimensions/params.py | 103 +------------- .../pseudo_three_dimensions/constants.py | 2 - .../pseudo_three_dimensions/params.py | 130 +----------------- geoapps/inversion/electricals/params.py | 72 ++++++---- 5 files changed, 55 insertions(+), 254 deletions(-) diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py index c9a94e140..3fea63dbb 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/constants.py @@ -60,7 +60,6 @@ "length_scale_z": 1.0, "s_norm": 0.0, "x_norm": 2.0, - "y_norm": 2.0, "z_norm": 2.0, "gradient_type": "total", "max_irls_iterations": 25, @@ -320,7 +319,6 @@ "reference_model": 1e-1, "s_norm": 0.0, "x_norm": 2.0, - "y_norm": 2.0, "z_norm": 2.0, "upper_bound": 100.0, "lower_bound": 1e-5, diff --git a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/params.py b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/params.py index 979e3eafc..2216c33f7 100644 --- a/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/params.py +++ b/geoapps/inversion/electricals/direct_current/pseudo_three_dimensions/params.py @@ -9,115 +9,34 @@ from copy import deepcopy -from geoapps.inversion import InversionBaseParams from geoapps.inversion.electricals.direct_current.pseudo_three_dimensions.constants import ( default_ui_json, forward_defaults, inversion_defaults, validations, ) +from geoapps.inversion.electricals.params import BasePseudo3DParams -class DirectCurrentPseudo3DParams(InversionBaseParams): +class DirectCurrentPseudo3DParams(BasePseudo3DParams): """ Parameter class for electrical->conductivity inversion. """ _physical_property = "conductivity" + _inversion_type = "direct current 3d" def __init__(self, input_file=None, forward_only=False, **kwargs): self._default_ui_json = deepcopy(default_ui_json) self._forward_defaults = deepcopy(forward_defaults) self._inversion_defaults = deepcopy(inversion_defaults) - self._inversion_type = "direct current 3d" self._validations = validations self._potential_channel_bool = None self._potential_channel = None self._potential_uncertainty = None - self._line_object = None - self._u_cell_size: float = (25.0,) - self._v_cell_size: float = (25.0,) - self._depth_core: float = (500.0,) - self._horizontal_padding: float = (1000.0,) - self._vertical_padding: float = (1000.0,) - self._expansion_factor: float = (1.3,) - self._files_only = None - self._cleanup = None super().__init__(input_file=input_file, forward_only=forward_only, **kwargs) - @property - def u_cell_size(self): - return self._u_cell_size - - @u_cell_size.setter - def u_cell_size(self, value): - self.setter_validator("u_cell_size", value) - - @property - def v_cell_size(self): - return self._v_cell_size - - @v_cell_size.setter - def v_cell_size(self, value): - self.setter_validator("v_cell_size", value) - - @property - def depth_core(self): - return self._depth_core - - @depth_core.setter - def depth_core(self, value): - self.setter_validator("depth_core", value) - - @property - def horizontal_padding(self): - return self._horizontal_padding - - @horizontal_padding.setter - def horizontal_padding(self, value): - self.setter_validator("horizontal_padding", value) - - @property - def vertical_padding(self): - return self._vertical_padding - - @vertical_padding.setter - def vertical_padding(self, value): - self.setter_validator("vertical_padding", value) - - @property - def expansion_factor(self): - return self._expansion_factor - - @expansion_factor.setter - def expansion_factor(self, value): - self.setter_validator("expansion_factor", value) - - @property - def inversion_type(self): - return self._inversion_type - - @inversion_type.setter - def inversion_type(self, val): - self.setter_validator("inversion_type", val) - - @property - def line_object(self): - return self._line_object - - @line_object.setter - def line_object(self, val): - self._line_object = val - - @property - def line_id(self): - return self._line_id - - @line_id.setter - def line_id(self, val): - self._line_id = val - @property def potential_channel_bool(self): return self._potential_channel_bool @@ -141,19 +60,3 @@ def potential_uncertainty(self): @potential_uncertainty.setter def potential_uncertainty(self, val): self.setter_validator("potential_uncertainty", val, fun=self._uuid_promoter) - - @property - def files_only(self): - return self._files_only - - @files_only.setter - def files_only(self, val): - self.setter_validator("files_only", val) - - @property - def cleanup(self): - return self._cleanup - - @cleanup.setter - def cleanup(self, val): - self.setter_validator("cleanup", val) diff --git a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py index 6c0e74141..d2019e6a5 100644 --- a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py +++ b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/constants.py @@ -61,7 +61,6 @@ "length_scale_z": 1.0, "s_norm": 0.0, "x_norm": 2.0, - "y_norm": 2.0, "z_norm": 2.0, "gradient_type": "total", "max_irls_iterations": 25, @@ -333,7 +332,6 @@ "conductivity_model": 0.1, "s_norm": 0.0, "x_norm": 2.0, - "y_norm": 2.0, "z_norm": 2.0, "upper_bound": 100.0, "lower_bound": 1e-5, diff --git a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/params.py b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/params.py index 4f0626d2b..77580a024 100644 --- a/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/params.py +++ b/geoapps/inversion/electricals/induced_polarization/pseudo_three_dimensions/params.py @@ -9,87 +9,35 @@ from copy import deepcopy -from geoapps.inversion import InversionBaseParams from geoapps.inversion.electricals.induced_polarization.pseudo_three_dimensions.constants import ( default_ui_json, forward_defaults, inversion_defaults, validations, ) +from geoapps.inversion.electricals.params import BasePseudo3DParams -class InducedPolarizationPseudo3DParams(InversionBaseParams): +class InducedPolarizationPseudo3DParams(BasePseudo3DParams): """ Parameter class for electrical->chargeability inversion. """ _physical_property = "chargeability" + _inversion_type = "induced polarization pseudo 3d" def __init__(self, input_file=None, forward_only=False, **kwargs): self._default_ui_json = deepcopy(default_ui_json) self._forward_defaults = deepcopy(forward_defaults) self._inversion_defaults = deepcopy(inversion_defaults) - self._inversion_type = "induced polarization pseudo 3d" self._validations = validations - self._potential_channel_bool = None - self._potential_channel = None - self._potential_uncertainty = None - self._u_cell_size: float = (25.0,) - self._v_cell_size: float = (25.0,) - self._depth_core: float = (500.0,) - self._horizontal_padding: float = (1000.0,) - self._vertical_padding: float = (1000.0,) self._conductivity_model: float | None = 1e-3 self._chargeability_channel_bool: bool = True self._chargeability_channel = None self._chargeability_uncertainty = None - self._expansion_factor: float = (1.3,) - self._line_object = None - self._files_only = None - self._cleanup = None super().__init__(input_file=input_file, forward_only=forward_only, **kwargs) - @property - def u_cell_size(self): - return self._u_cell_size - - @u_cell_size.setter - def u_cell_size(self, value): - self.setter_validator("u_cell_size", value) - - @property - def v_cell_size(self): - return self._v_cell_size - - @v_cell_size.setter - def v_cell_size(self, value): - self.setter_validator("v_cell_size", value) - - @property - def depth_core(self): - return self._depth_core - - @depth_core.setter - def depth_core(self, value): - self.setter_validator("depth_core", value) - - @property - def horizontal_padding(self): - return self._horizontal_padding - - @horizontal_padding.setter - def horizontal_padding(self, value): - self.setter_validator("horizontal_padding", value) - - @property - def vertical_padding(self): - return self._vertical_padding - - @vertical_padding.setter - def vertical_padding(self, value): - self.setter_validator("vertical_padding", value) - @property def conductivity_model(self): return self._conductivity_model @@ -121,75 +69,3 @@ def chargeability_uncertainty(self): @chargeability_uncertainty.setter def chargeability_uncertainty(self, value): self.setter_validator("chargeability_uncertainty", value) - - @property - def expansion_factor(self): - return self._expansion_factor - - @expansion_factor.setter - def expansion_factor(self, value): - self.setter_validator("expansion_factor", value) - - @property - def inversion_type(self): - return self._inversion_type - - @inversion_type.setter - def inversion_type(self, val): - self.setter_validator("inversion_type", val) - - @property - def line_object(self): - return self._line_object - - @line_object.setter - def line_object(self, val): - self._line_object = val - - @property - def line_id(self): - return self._line_id - - @line_id.setter - def line_id(self, val): - self._line_id = val - - @property - def potential_channel_bool(self): - return self._potential_channel_bool - - @potential_channel_bool.setter - def potential_channel_bool(self, val): - self.setter_validator("potential_channel_bool", val) - - @property - def potential_channel(self): - return self._potential_channel - - @potential_channel.setter - def potential_channel(self, val): - self.setter_validator("potential_channel", val, fun=self._uuid_promoter) - - @property - def potential_uncertainty(self): - return self._potential_uncertainty - - @potential_uncertainty.setter - def potential_uncertainty(self, val): - self.setter_validator("potential_uncertainty", val, fun=self._uuid_promoter) - - @property - def files_only(self): - return self._files_only - - @files_only.setter - def files_only(self, val): - self.setter_validator("files_only", val) - - @property - def cleanup(self): - return self._cleanup - - @cleanup.setter - def cleanup(self, val): - self.setter_validator("cleanup", val) diff --git a/geoapps/inversion/electricals/params.py b/geoapps/inversion/electricals/params.py index 9c41e55c6..de40d4fda 100644 --- a/geoapps/inversion/electricals/params.py +++ b/geoapps/inversion/electricals/params.py @@ -8,27 +8,19 @@ from __future__ import annotations import numpy as np +from geoh5py.data import Data, DataAssociationEnum from geoapps.inversion import InversionBaseParams from geoapps.utils.models import get_drape_model -class Base2DParams(InversionBaseParams): +class Core2DParams(InversionBaseParams): """ - Parameter class for electrical->induced polarization (IP) inversion. + Core parameter class for 2D electrical->conductivity inversion. """ - _directive_list = [ - "UpdateSensitivityWeights", - "Update_IRLS", - "BetaEstimate_ByEig", - "UpdatePreconditioner", - "SaveIterationsGeoH5", - ] - def __init__(self, input_file=None, forward_only=False, **kwargs): self._line_object = None - self._line_id = None self._u_cell_size: float = 25.0 self._v_cell_size: float = 25.0 self._depth_core: float = 100.0 @@ -47,14 +39,8 @@ def line_object(self): def line_object(self, val): self._line_object = val - @property - def line_id(self): - """Line ID to invert.""" - return self._line_id - - @line_id.setter - def line_id(self, val): - self._line_id = val + if isinstance(val, Data) and val.association is not DataAssociationEnum.CELL: + raise ValueError("Line identifier must be associated with cells.") @property def u_cell_size(self): @@ -110,6 +96,26 @@ def expansion_factor(self): def expansion_factor(self, value): self.setter_validator("expansion_factor", value) + +class Base2DParams(Core2DParams): + """ + Parameter class for electrical->induced polarization (IP) inversion. + """ + + def __init__(self, input_file=None, forward_only=False, **kwargs): + self._line_id = None + + super().__init__(input_file=input_file, forward_only=forward_only, **kwargs) + + @property + def line_id(self): + """Line ID to invert.""" + return self._line_id + + @line_id.setter + def line_id(self, val): + self._line_id = val + @property def mesh(self): if self._mesh is None and self.geoh5 is not None: @@ -136,10 +142,30 @@ def mesh(self): def mesh(self, val): self.setter_validator("mesh", val, fun=self._uuid_promoter) + +class BasePseudo3DParams(Core2DParams): + """ + Base parameter class for pseudo electrical->conductivity inversion. + """ + + def __init__(self, input_file=None, forward_only=False, **kwargs): + self._files_only = None + self._cleanup = None + + super().__init__(input_file=input_file, forward_only=forward_only, **kwargs) + @property - def length_scale_y(self): - return None + def files_only(self): + return self._files_only + + @files_only.setter + def files_only(self, val): + self.setter_validator("files_only", val) @property - def y_norm(self): - return None + def cleanup(self): + return self._cleanup + + @cleanup.setter + def cleanup(self, val): + self.setter_validator("cleanup", val) From 9ae257e4835187c899339934ded51dc2f9b0dac8 Mon Sep 17 00:00:00 2001 From: dominiquef Date: Tue, 13 Feb 2024 13:27:07 -0800 Subject: [PATCH 17/17] Update FlinFlon_dcip, fix issue with locs --- geoapps-assets/FlinFlon_dcip.geoh5 | 4 ++-- geoapps/utils/surveys.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/geoapps-assets/FlinFlon_dcip.geoh5 b/geoapps-assets/FlinFlon_dcip.geoh5 index 6387bb25f..ae14309ed 100644 --- a/geoapps-assets/FlinFlon_dcip.geoh5 +++ b/geoapps-assets/FlinFlon_dcip.geoh5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:5cbde9e640dbba9f4b82c19ed0a990a042e429d22c7d901fb0de56ddb7e6c230 -size 2759487 +oid sha256:2ab8fe351e9aba06dbe7149307ab5169b6c1f6d08af27c2f8b294451ba6a7611 +size 2438535 diff --git a/geoapps/utils/surveys.py b/geoapps/utils/surveys.py index dff4a4f8f..de6ef5097 100644 --- a/geoapps/utils/surveys.py +++ b/geoapps/utils/surveys.py @@ -50,7 +50,7 @@ def get_containing_cells( inds = np.unique(np.r_[inds, np.hstack(line_ind)]) elif isinstance(mesh, TensorMesh): - locations = data.drape_locations(get_unique_locations(data.survey)) + locations = data.drape_locations(np.unique(data.locations, axis=0)) xi = np.searchsorted(mesh.nodes_x, locations[:, 0]) - 1 yi = np.searchsorted(mesh.nodes_y, locations[:, -1]) - 1 inds = xi * mesh.shape_cells[1] + yi