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