PyHydroGeophysX.inversion package

Submodules

PyHydroGeophysX.inversion.base module

Base classes for geophysical inversion frameworks.

This module defines:

• InversionResult: A base class to store and manage common results from an inversion process, including final model, predicted data, convergence history, and plotting utilities.

• TimeLapseInversionResult: A specialized version of InversionResult for time-lapse studies, handling multiple models over time and providing time-slice plotting and animation.

• InversionBase: An abstract base class outlining the common structure and interface for various geophysical inversion methods (e.g., ERT, SRT). It handles data, mesh, and basic parameter management.

class PyHydroGeophysX.inversion.base.InversionBase(data: pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: object

Abstract base class for geophysical inversion methods.

This class provides a foundational structure for specific inversion techniques (e.g., ERT, SRT). It manages observed data, mesh, common inversion parameters, and an InversionResult object. Subclasses must implement methods for setting up the inversion, running the inversion loop, computing Jacobians, and evaluating the objective function.

compute_jacobian(model: ndarray) → ndarray

Abstract method to compute the Jacobian matrix (sensitivity matrix).

The Jacobian J_ij = ∂d_i / ∂m_j relates changes in model parameters (m_j) to changes in observed data (d_i).

Parameters:

model (np.ndarray) – The current model parameter vector for which to compute the Jacobian.

Returns:,

np.ndarray: The computed Jacobian matrix, typically of shape (n_data, n_model_params).

objective_function(model: ndarray, data_to_fit: ndarray | None = None) → float

Abstract method to compute the value of the objective function.

The objective function typically includes a data misfit term and one or more regularization terms: Φ(m) = Φ_d(m) + λ * Φ_m(m).

Parameters:

• model (np.ndarray) – The current model parameter vector.

• data_to_fit (Optional[np.ndarray], optional) – The observed data to fit. If None, self.data (from DataContainer) is typically used. Defaults to None.

Returns:

The calculated value of the objective function.

Return type:

float

run() → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.base.InversionResult

Bases: object

Base class to store, save, load, and plot results from a geophysical inversion.

final_model

The final inverted model parameters (e.g., resistivity, velocity). Typically a 1D array corresponding to mesh cells.

Type:

Optional[np.ndarray]

predicted_data

The data predicted by the forward model using the final_model.

Type:

Optional[np.ndarray]

coverage

Coverage or sensitivity values for the model parameters, often derived from the Jacobian or resolution matrix.

Type:

Optional[np.ndarray]

mesh

The PyGIMLi mesh object used in the inversion.

Type:

Optional[pg.Mesh]

iteration_models

A list storing the model parameters at each iteration of the inversion.

Type:

List[np.ndarray]

iteration_data_errors

A list storing the data misfit (e.g., residuals) at each iteration.

Type:

List[np.ndarray]

iteration_chi2

A list storing the chi-squared (χ²) value or a similar misfit metric at each iteration.

Type:

List[float]

meta

A dictionary to store any additional metadata about the inversion run (e.g., inversion parameters, timings, comments).

Type:

Dict[str, Any]

classmethod load(filename: str) → InversionResult

Load inversion results from a file previously saved by the save method.

Parameters:

filename (str) – The base path to the saved results file. If saved as name.pkl and name.bms, provide name.pkl or name.

Returns:

An instance of InversionResult (or a subclass if called on one)

populated with the loaded data.

Return type:

InversionResult

Raises:

• FileNotFoundError – If the main data file or associated mesh file (if referenced) is not found.

• IOError – If there’s an error during file reading.

• pickle.UnpicklingError – If the file cannot be unpickled.

plot_convergence(ax: Axes | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the convergence curve (chi-squared misfit vs. iteration number).

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• **kwargs (Any) – Additional keyword arguments passed directly to ax.plot (e.g., color, marker, linestyle).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.iteration_chi2 is empty (no convergence data).

plot_model(ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the final inverted model on its associated mesh.

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• cmap (str, optional) – The colormap to use for visualizing model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – If provided, cells with coverage values below this threshold will be masked (made transparent or semi-transparent) in the plot. Requires self.coverage to be populated. Defaults to None (no masking).

• **kwargs (Any) – Additional keyword arguments passed directly to pygimli.show (e.g., cMin, cMax, orientation, logScale).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.final_model or self.mesh is None.

save(filename: str) → None

Save the inversion results to a file using Python’s pickle format. If a mesh is present, it is saved separately as a PyGIMLi binary mesh file (.bms).

Parameters:

filename (str) – The base path (including filename without extension) to save the results. The main data will be saved as filename.pkl (or just filename if user includes .pkl). The mesh will be saved as filename.bms or filename.pkl.bms. It’s recommended to provide filename without .pkl.

Raises:

• IOError – If there’s an error during file writing.

• pickle.PicklingError – If an object cannot be pickled.

class PyHydroGeophysX.inversion.base.TimeLapseInversionResult

Bases: InversionResult

Specialized class to store and manage results from time-lapse inversions. Inherits from InversionResult and adds attributes specific to time-lapse data.

final_models

A 2D NumPy array where each column represents the inverted model for a specific timestep (shape: num_cells x num_timesteps).

Type:

Optional[np.ndarray]

timesteps

An array or list of time values (e.g., hours, days) corresponding to each model slice in final_models.

Type:

Optional[np.ndarray]

all_coverage

A list where each element is the coverage array for the corresponding timestep’s model.

Type:

List[np.ndarray]

all_chi2

A list where each element is a list of chi-squared values per iteration for the inversion of that specific timestep or window. (Note: Original was List[float], if chi2 is per window, might need adjustment)

Type:

List[Any]

create_time_lapse_animation(output_filename: str, cmap: str = 'viridis', coverage_threshold: float | None = None, dpi: int = 100, fps: int = 2, **kwargs: Any) → None

Create and save an animation (e.g., MP4 video) of the time-lapse inversion results.

Requires ffmpeg or another Matplotlib-supported animation writer to be installed.

Parameters:

• output_filename (str) – The filename for the output animation (e.g., ‘timelapse_animation.mp4’).

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• dpi (int, optional) – Dots Per Inch for the output animation. Defaults to 100.

• fps (int, optional) – Frames Per Second for the animation. Defaults to 2.

• **kwargs (Any) – Additional keyword arguments passed to plot_time_slice for each frame.

Raises:

• ValueError – If final_models or mesh is missing.

• ImportError – If matplotlib.animation cannot be imported.

classmethod load(filename: str) → TimeLapseInversionResult

Load time-lapse inversion results. Overrides base class load.

plot_time_slice(timestep_idx: int, ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot a single time slice (inverted model at a specific timestep) from the results.

Parameters:

• timestep_idx (int) – The zero-based index of the timestep to plot.

• ax (Optional[plt.Axes], optional) – Matplotlib Axes to plot on. If None, creates new.

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• **kwargs (Any) – Additional arguments passed to pg.show.

Returns:

The matplotlib Figure and Axes objects.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If models or mesh are missing, or if timestep_idx is out of range.

save(filename: str) → None

Save time-lapse inversion results. Overrides base class save.

PyHydroGeophysX.inversion.cross_constraints module

Cross-method constraint utilities for joint/cooperative inversion.

These utilities enable information sharing between geophysical methods and provide hydrology-to-geophysics coupling helpers.

class PyHydroGeophysX.inversion.cross_constraints.PetrophysicalCoupling

Bases: object

Coupling helpers from hydrological state to multi-method geophysics.

static compare_inversions_to_hydro(hydro_wc, ert_result, srt_result=None, em_result=None, petro_params=None)

Compare inversion products back to hydrological water content.

Returns misfit statistics for available methods.

static water_content_to_all_geophysics(water_content, porosity, **params)

Convert water content to resistivity, velocity, and conductivity.

Uses existing petrophysical functions in the package.

class PyHydroGeophysX.inversion.cross_constraints.StructuralConstraint

Bases: object

Build structural constraint terms from one method for another.

static apply_structural_weights_to_Wm(Wm, boundary_weights)

Modify an existing smoothness matrix Wm to respect boundaries.

static build_cross_gradient_operator(mesh, model_a, model_b)

Build a cross-gradient operator surrogate.

Returns a diagonal sparse matrix weighted by local cross-gradient magnitude, where small values indicate better structural agreement.

static build_linearized_cross_gradient_blocks(RCM: ndarray, X: ndarray, model_a: ndarray, model_b: ndarray, mode: str = 'direct') → Tuple[ndarray, ndarray]

Build linearized cross-gradient blocks B1 and B2.

B1 multiplies model_a and B2 multiplies model_b. The resulting penalty terms are ||B1 m_a||^2 and ||B2 m_b||^2.

static build_local_design_matrix(mesh) → ndarray

Build the local design matrix X=[x, z, 1] used by cross-gradient.

This mirrors the legacy xyzpos[:, 2] = 1 construction.

static build_neighborhood_matrix(mesh, Wm: Any | None = None, source: str = 'smoothness', correlation_lengths: Sequence[float] = (4.0, 4.0), threshold: float = 0.01, binarize: bool = True) → ndarray

Build the neighborhood/correlation matrix used by cross-gradient.

Parameters:

• mesh – Mesh object.

• Wm – Smoothness matrix. Required when source='smoothness'.

• source – One of 'smoothness' (legacy direct mode) or 'covariance' (legacy spatial mode).

• correlation_lengths – Correlation lengths passed to pg.utils.covarianceMatrix when source='covariance'.

• threshold – Entries with absolute value below this threshold are zeroed.

• binarize – If True, convert nonzero entries to 1.

static from_conductivity_model(conductivity, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an EM conductivity model.

static from_velocity_model(velocity_model, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an SRT velocity model.

Returns a pairwise-boundary weight map that can reduce smoothing across detected structural interfaces.

PyHydroGeophysX.inversion.em1d module

Occam-style 1D FDEM/TDEM inversion helpers.

PyHydroGeophysX.inversion.em1d.build_sounding_block(data: Dict[str, Any], geom: Dict[str, Any], inv: Dict[str, Any], method: str = 'TDEM', *, position: float = 0.0, line: int = 0, label: str = '')

Package one sounding for the coupled line inversion.

The observed vector, uncertainty, and forward operator are built by the same code the per-sounding inversion uses, so a station fits the same data whether it is solved alone or as part of a line. The Jacobian comes from SimPEG’s analytic sensitivity rather than from finite differences.

PyHydroGeophysX.inversion.em1d.fdem_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert an FDEM sounding for a layered resistivity model (Occam 1D).

PyHydroGeophysX.inversion.em1d.tdem_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert a TDEM sounding for a layered resistivity model (Occam 1D).

PyHydroGeophysX.inversion.em1d.tdem_joint_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert all available LM/HM gates at one station for one shared 1D model.

PyHydroGeophysX.inversion.em1d.tdem_moment_blocks(data: Dict[str, Any], geom: Dict[str, Any], inv: Dict[str, Any], thick: ndarray) → List[Dict[str, Any]]

One forward block per usable moment at a TDEM station.

Stations carrying separate LM and HM gate sets produce one block each; a plain single-response station produces one block named TDEM. Returning the same shape for both is what lets the per-sounding inversion and the coupled line inversion share this assembly instead of each writing its own copy of the uncertainty and gate-selection rules.

PyHydroGeophysX.inversion.em1d_lci module

Simultaneous laterally-constrained inversion (LCI) for 1D EM soundings.

Every sounding on a line is solved in one system. The model vector stacks the per-sounding layer models, the forward operator is block diagonal (a sounding only sees its own layers), and an explicit coupling operator ties neighbouring soundings together. That is the same structure as the time-lapse inversion in PyHydroGeophysX.inversion.time_lapse, with along-line distance in place of time.

This replaces block-coordinate LCI, where each sounding is re-inverted on its own against a reference built from its neighbours’ models from the previous pass. Under block coordinates the lateral constraint is never enforced while a sounding is being solved, so the passes chase each other; here it is part of the system being solved.

Two properties make the simultaneous form affordable:

• SimPEG’s Simulation1DLayered.getJ returns the analytic sensitivity, and measures faster than a single forward call. Finite differencing the same Jacobian costs one forward per layer, so the analytic route is roughly an order of magnitude cheaper per Gauss-Newton iteration.

• The normal matrix is block tridiagonal, so a sparse factorization scales with the number of soundings rather than its cube.

Models are parameterized as x = log10(resistivity), matching the Occam routine in PyHydroGeophysX.inversion.em1d so that smoothness and lateral_smoothness keep the meaning they have there.

PyHydroGeophysX.inversion.em1d_lci.CHI2_EQUIVALENCE: float = 0.02

Relative chi-squared gain a rougher model has to earn to be preferred. When the target misfit is out of reach — noisy ground data with model error well above the assumed gate errors — the search keeps relaxing the smoothness for gains in the third decimal place, and hands back a railed model that fits no better than the smooth one. Trials within this margin of the best misfit are treated as the same fit, and the smoothest of them wins.

PyHydroGeophysX.inversion.em1d_lci.DOI_SENSITIVITY_THRESHOLD: float = 0.8

Cumulated sensitivity a depth has to carry to count as investigated.

The published value from Christiansen and Auken (2012), who fine-tuned it across ground conductivity meters, DC soundings and airborne TEM and report 0.6 to 1.2 as the range they considered, moving their example by about 15 %. Because their measure lives in logarithmic data and model space it carries no units of its own, which is what lets one number serve every system: at the depth of investigation, moving every layer below it by one e-fold in resistivity shifts the predicted response by 0.8 error bars in total.

The threshold is tied to the error model by construction. Doubling the assumed error halves the sensitivity and the section gets shallower, which is the honest response to noisier data. Raise it for a more conservative picture, lower it to see what the deeper part of the model looks like.

On sparse ground TDEM the published value can saturate. With a handful of gates per station and a model of twenty layers, the deepest layer often clears 0.8 on its own: the measure cumulates from the bottom up, and that layer is thick. The reported depth then collapses onto the bottom of the parameterisation for much of the survey, which is the measure meeting a coarse deep grid rather than a claim about resolution. Two symptoms identify it: a large share of stations reporting exactly the model bottom, and stations holding three gates reporting the same depth as stations holding ten. Values in the 6 to 8 range keep the reported depth inside the model on such data. Christiansen and Auken’s value is the default because it is the published one and it travels across systems.

class PyHydroGeophysX.inversion.em1d_lci.LCIResult(models: ndarray, chi2: float, chi2_per_sounding: ndarray, chi2_history: List[float] = <factory>, iterations: int = 0, stop_reason: str = '', smoothness_scale: float = 1.0, lambda_vertical: float = 0.0, lambda_lateral: float = 0.0, n_data: int = 0, seconds: float = 0.0, lambda_search: Dict[str, ~typing.Any] | None=None, diagnostics: Dict[str, ~typing.Any]=<factory>, chi2_median_history: List[float] = <factory>)

Bases: object

Outcome of one coupled line inversion.

as_dict() → Dict[str, Any]

chi2: float

chi2_history: List[float]

chi2_median_history: List[float]

chi2_per_sounding: ndarray

diagnostics: Dict[str, Any]

iterations: int = 0

lambda_lateral: float = 0.0

lambda_search: Dict[str, Any] | None = None

lambda_vertical: float = 0.0

models: ndarray

n_data: int = 0

seconds: float = 0.0

smoothness_scale: float = 1.0

stop_reason: str = ''

PyHydroGeophysX.inversion.em1d_lci.LOG_RESISTIVITY_BOUNDS: Tuple[float, float] = (0.0, 5.0)

log10 resistivity bounds, matching _occam_1d (1 to 1e5 ohm-m).

PyHydroGeophysX.inversion.em1d_lci.SMOOTHNESS_SCALE_BOUNDS: Tuple[float, float] = (0.0001, 10000.0)

Bounds on the smoothness scale the chi-squared search may visit.

class PyHydroGeophysX.inversion.em1d_lci.SoundingBlock(forward: Callable[[ndarray], ndarray], jacobian: Callable[[ndarray], ndarray], dobs: ndarray, uncertainty: ndarray, position: float = 0.0, line: int = 0, label: str = '', prior_lower: ndarray | None = None, prior_weights: ndarray | None = None)

Bases: object

One sounding’s contribution to the coupled system.

forward(sigma) returns the predicted data for a conductivity model and jacobian(sigma) its derivative with respect to conductivity, shaped (n_data, n_layers). Keeping both as callables lets FDEM, single-moment TDEM, and joint LM+HM stations share one solver: the caller decides how the response is assembled, the solver only needs the pair.

dobs: ndarray

forward: Callable[[ndarray], ndarray]

jacobian: Callable[[ndarray], ndarray]

label: str = ''

line: int = 0

position: float = 0.0

prior_lower: ndarray | None = None

prior_weights: ndarray | None = None

uncertainty: ndarray

PyHydroGeophysX.inversion.em1d_lci.cumulated_sensitivity(block: SoundingBlock, model: ndarray) → ndarray

Cumulated sensitivity after Christiansen and Auken (2012).

Their construction, equations 2, 3 and 5 of A global measure for depth of investigation, GEOPHYSICS 77(4), WB171-WB177:

• G_ij = d log(data_i) / d log(rho_j), the Jacobian of the final model in logarithmic data and model space. Working in logarithms on both sides is what makes the resulting number comparable between data types, and so lets one absolute threshold serve every system.

• s_j = sum_i |G_ij| / sigma_i, summed over all N data with sigma_i the standard deviation of the log datum, which is its relative error.

• S_j = sum_{k >= j} s_k, cumulated from the bottom layer upward. Entry j is therefore the total information the data carry about layer j and everything below it, counted in error bars.

Their equation 4 divides s by the layer thickness; the paper uses that only for plotting, and the cumulated quantity here is built from equation 3, which is also what keeps it independent of the layer grid: split a layer in two and its sensitivity splits with it, so the value at a given depth does not move (measured at 0.2 to 2.5 % across a 2x refinement on real ground TDEM).

Only the data part of the Jacobian takes part, so a depth that clears the threshold is one the measurements reach, not one the lateral or vertical constraint filled in. The Jacobian is SimPEG’s analytic sensitivity, the same one the coupled solver uses, so this costs about one forward evaluation per sounding.

Their step 2, sub-discretizing a few-layer model before differentiating, is unnecessary here: the paper skips it for smooth models, and these are solved on a fixed grid of a dozen layers or more.

PyHydroGeophysX.inversion.em1d_lci.invert_lci(blocks: ~typing.Sequence[~PyHydroGeophysX.inversion.em1d_lci.SoundingBlock], n_layers: int, *, auto_lambda: bool = True, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_lambda_trials: int = 5, smoothness_scale: float = 1.0, scale_bounds: ~typing.Tuple[float, float] = (0.0001, 10000.0), verbose: bool = True, log: ~typing.Callable[[str], None] = <function noop>, **kwargs: ~typing.Any) → LCIResult

Run the coupled line inversion, relaxing the smoothness if needed.

The fixed-smoothness run happens first and is always kept. When it misses target_chi2 by more than chi2_tolerance and auto_lambda is set, a bracket-and-bisect search over a single scale multiplying both smoothness terms looks for a better fit. Each trial warm-starts from the models of the nearest scale already solved, which is what keeps the search affordable.

The returned result is whichever run landed closest to the target. When no trial reached the target band, trials whose misfit is within CHI2_EQUIVALENCE of the best are treated as equally good fits and the smoothest of them is returned instead of the roughest. The trial record lives in lambda_search.

PyHydroGeophysX.inversion.em1d_lci.invert_lci_rejecting_outliers(blocks: ~typing.Sequence[~PyHydroGeophysX.inversion.em1d_lci.SoundingBlock], n_layers: int, *, threshold: float = 3.0, passes: int = 2, min_fraction: float = 0.5, min_gates: int = 3, log: ~typing.Callable[[str], None] = <function noop>, **kwargs: ~typing.Any) → Tuple[LCIResult, List[SoundingBlock], Dict[str, Any]]

Solve the line, then drop the gates the model cannot explain and re-solve.

This is the EM counterpart of the ERT outlier pass. Each cycle removes the gates whose weighted residual exceeds threshold and solves again, warm started from the model just found and at the smoothness the first solve settled on, so a rejection pass costs far less than the first one. Two floors bound the cut. min_fraction is the survey-wide one: when more gates exceed the threshold than it allows, the pass drops the worst offenders up to the limit rather than refusing, because a bad fit is exactly when rejection is wanted. min_gates is the per-sounding one: a station never loses so many gates that fewer than this remain (or all of them, if it arrived with fewer). Without it a station holding one or two gates loses them both on the first pass, and the section then has a hole where only the lateral constraint is left holding the model. Its best-fitting gates are the ones kept.

A single noisy gate on a TDEM station carries a lot of weight (a station may hold only a handful), so cutting is per gate rather than per sounding.

Returns (outcome, kept_blocks, info).

PyHydroGeophysX.inversion.em1d_lci.invert_lci_with_robust_errors(blocks: ~typing.Sequence[~PyHydroGeophysX.inversion.em1d_lci.SoundingBlock], n_layers: int, *, threshold: float = 3.0, passes: int = 3, max_error_factor: float = 10.0, min_unchanged_fraction: float = 0.0, error_target_chi2: float = 0.0, target_tolerance: float = 0.25, log: ~typing.Callable[[str], None] = <function noop>, **kwargs: ~typing.Any) → Tuple[LCIResult, List[SoundingBlock], Dict[str, Any]]

Warm-started error reweighting, preserving all LM/HM data and operators.

The returned blocks carry effective errors (also used for DOI). Input blocks are untouched. The solver outcome uses effective chi2; the report additionally supplies original-error scores for honest comparisons between runs.

PyHydroGeophysX.inversion.em1d_lci.lateral_edges(positions: Sequence[float], lines: Sequence[int] | None = None, *, reference_distance: float = 10.0, distance_power: float = 1.0) → List[Tuple[int, int, float]]

Return (i, j, weight) for each pair of neighbouring soundings.

Soundings are neighbours when they are adjacent in along-line order on the same survey line. The penalty a pair produces scales as (reference_distance / d) ** distance_power, so at the default power of one, two stations 10 m apart are tied ten times as tightly as two 100 m apart. The weight is the square root of that, because the penalty is the square of the weighted difference. It is capped at 1 so that a pair closer than reference_distance is not tied arbitrarily hard.

distance_power of 0 removes the distance dependence and ties every neighbouring pair alike. Values between 0 and 1 loosen the fall-off, which suits a line whose station spacing varies enough that the linear rule leaves the widely spaced pairs effectively unconstrained.

PyHydroGeophysX.inversion.em1d_lci.mask_sounding_block(block: SoundingBlock, keep: Sequence[bool]) → SoundingBlock

A copy of block carrying only the gates flagged in keep.

The forward and Jacobian are wrapped rather than rebuilt, so the SimPEG simulation behind them is reused as is and dropping a gate costs nothing. A block may end up with no data at all: the sounding then contributes no residual and its model comes entirely from its neighbours through the lateral constraint, which is the honest outcome when every one of its gates was rejected.

PyHydroGeophysX.inversion.em1d_lci.sensitivity_doi(block: SoundingBlock, model: ndarray, depth_edges: ndarray, *, threshold: float = 0.8) → float

Depth of investigation: the bottom of the deepest layer still resolved.

Returns 0.0 when even the shallowest layer misses the threshold, which is the honest answer for a station whose gates were all rejected: it has no depth of investigation, and its column carries only what its neighbours say.

PyHydroGeophysX.inversion.em1d_lci.solve_lci(blocks: ~typing.Sequence[~PyHydroGeophysX.inversion.em1d_lci.SoundingBlock], n_layers: int, *, smoothness: float = 0.3, lateral_smoothness: float = 0.3, smoothness_scale: float = 1.0, reference_distance: float = 10.0, lateral_distance_power: float = 1.0, initial_model: ~numpy.ndarray | None = None, starting_resistivity: float = 100.0, max_iterations: int = 20, convergence_tolerance: float = 0.02, convergence_metric: str = 'data', solver: str = 'trf', trf_max_nfev: int = 90, trf_ftol: float = 0.0001, trf_xtol: float = 1e-06, trf_gtol: float = 1e-05, min_iterations: int = 2, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, line_search_steps: int = 6, target_steps: int = 4, bounds: ~typing.Tuple[float, float] = (0.0, 5.0), parallel_workers: int = 0, verbose: bool = True, log: ~typing.Callable[[str], None] = <function noop>) → LCIResult

Solve one coupled line inversion at a fixed smoothness.

Bound-aware sparse trust-region least squares (solver='trf') is the formal default. solver='gauss_newton' selects the fast legacy path. Its budget is trf_max_nfev (forward evaluations, including rejected trials), not max_iterations. Its three tolerances apply to the full objective, step and gradient, respectively. The legacy gauss_newton path retains its iteration/target controls for backwards compatibility.

The objective is

||W (F(x) - d)||^2 + lam_v ||Dz x||^2 + lam_l ||Dx x||^2

with lam_v = (smoothness_scale * smoothness)^2 and lam_l likewise from lateral_smoothness. Squaring keeps the two knobs numerically equivalent to the residual-stacking convention used by the per-sounding Occam routine, so an existing smoothness=0.3 setting means the same amount of vertical damping here.

On the legacy Gauss-Newton path, iteration stops at target_chi2, when the relative chi-squared improvement falls below convergence_tolerance (the plateau rule used throughout this package), when the line search cannot find a descent step, or at max_iterations. TRF instead uses its full- objective tolerances and forward-evaluation budget. Either path reports the reason in stop_reason.

A Gauss-Newton step near the target routinely shoots past it, so a run that stopped at the first iterate below target_chi2 would report whatever overshoot the step happened to produce. That both fits noise and makes the misfit a jumpy function of the smoothness, which the search in invert_lci() then cannot bisect. Set target_steps above zero and an overshooting step is shortened until the misfit lands within chi2_tolerance of the target instead.

parallel_workers runs the per-sounding forward and Jacobian on a thread pool, 0 choosing a count from the machine. Every sounding owns its forward operator, and the workers only read the shared model vector, so the arithmetic is the same either way; measured against the serial path, the models come back bit for bit identical.

PyHydroGeophysX.inversion.em1d_lci.weighted_residuals(blocks: Sequence[SoundingBlock], models: ndarray) → ndarray

(predicted - observed) / uncertainty for every gate on the line.

PyHydroGeophysX.inversion.em1d_priors module

Empirical resistive-background regularisation, separate from data errors.

A sustained weak absolute LM response may be consistent with a resistive half-space/background, but it does not identify a shallow layer or its depth. The optional prior therefore applies a weak, one-sided tendency to the whole 1-D model. It is not independent geological evidence and may also be triggered by interference, coupling or acquisition changes. It never enters data chi2 or data-only DOI.

The historical shallow_prior_* option names remain readable so existing project files and notebooks do not break. Their current physical interpretation is the background tendency described above, not a depth-local constraint.

PyHydroGeophysX.inversion.em1d_priors.raw_lm_quality_rows(datasets, reference_gate=2)

Recover identical diagnostics from live projects and saved Qt inputs.

The stored STD gives an absolute uncertainty proxy, not an independent measurement of ambient noise. Flags do not change this diagnostic gate.

PyHydroGeophysX.inversion.em1d_priors.resistive_prior_target(inv)

Return (reference, target, factor, source) in ohm-m.

The automatic reference is the effective starting half-space supplied by the workflow. A warm-started model uses its geometric-median resistivity. Standalone callers fall back to starting_model and then starting_resistivity. An explicit positive historical shallow_prior_min_resistivity still overrides the automatic target.

PyHydroGeophysX.inversion.em1d_priors.shallow_prior_scores(datasets, positions, lines, inv, quality_rows=None, signal_thresholds=None)

Fixed, pre-fit spatial scores from early LM SNR, reset at survey lines.

Baseline is the median of the first full window in acquisition-distance order. Raw diagnostics additionally require signal decline at the same physical gate time, without a large uncertainty increase. Without raw diagnostics, missing imported LM counts as zero under the user’s empirical assumption; this weaker fallback is identified in the report. Unreadable stations remain unavailable.

PyHydroGeophysX.inversion.em1d_priors.shallow_prior_terms(inv, thicknesses)

Whole-model lower log10(rho) tendency with grid-normalised weight.

thicknesses is accepted to retain the established public function signature. Every layer, including the basal half-space, receives the same coefficient. Dividing by sqrt(n_layers) keeps the penalty on a uniform model unchanged when the grid is split into more layers.

PyHydroGeophysX.inversion.em1d_priors.shallow_signal_thresholds(datasets, geometries, inv)

Absolute raw-signal limits, fixed by a manual value or a forward model.

Automatic mode models a homogeneous half-space at a fixed signal-reference resistivity, with the same waveform, filters, gate window and station geometry as the fit. This calibrates a heuristic trigger, NOT the depth of a resistive layer. No inverse data_scale enters: both observed diagnostics and limits are in the project’s stored normalized-response units. Operator caches are reused.

PyHydroGeophysX.inversion.ert_inversion module

Single-time ERT inversion functionality.

PyHydroGeophysX.inversion.ert_inversion.ERROR_SOURCES = ('file', 'estimate', 'max')

How err is chosen. file trusts the err column the instrument wrote and only estimates where it is missing; estimate always recomputes from relative_error/absolute_error; max takes the larger of the two per datum, which is the conservative reading when the file’s errors look optimistic. file is the default because silently discarding measured errors makes chi2 report on an error model the data never had.

class PyHydroGeophysX.inversion.ert_inversion.ERTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Single-time ERT inversion class.

run(initial_model: ndarray | None = None, reference_model: ndarray | None = None) → InversionResult

Run ERT inversion.

Parameters:

• initial_model – Initial model parameters (if None, a homogeneous model is used)

• reference_model – Model the regularization pulls toward. Defaults to the initial model, which is the usual smoothness-from-homogeneous setup. Pass it explicitly when continuing an interrupted run. Leaving it to default makes the restart point the new reference, so the penalty becomes roughness of the change since the restart rather than roughness of the model; accumulated structure then goes unpenalized and the run drifts steadily under-regularized. Measured on the Ashton line over a 20 to 0.3 lambda ladder, that reached chi2 5.72 where the same ladder with the reference pinned reached 6.38.

Returns:

InversionResult with inversion results

setup()

Set up ERT inversion (create operators, matrices, etc.)

class PyHydroGeophysX.inversion.ert_inversion.ERTRun(lam: float, chi2: float, iterations: int, stop: str, convergence: ~typing.List[float], model: ~numpy.ndarray, response: ~numpy.ndarray, mesh: ~typing.Any, coverage: ~numpy.ndarray | None = None, manager: ~typing.Any = None, metrics: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: object

One inversion at one lambda, described the same way by either engine.

chi2: float

convergence: List[float]

coverage: ndarray | None = None

iterations: int

lam: float

manager: Any = None

mesh: Any

metrics: Dict[str, Any]

model: ndarray

response: ndarray

property result

A manager-like handle for the viewer and the VTK export.

stop: str

PyHydroGeophysX.inversion.ert_inversion.GEOMETRIC_FACTOR_POLICIES = ('off', 'check', 'fix')

What to do about geometric factors before anything else touches the data. fix validates and, if the check fails, recomputes k numerically on the inversion mesh and rebuilds rhoa from the measured transfer resistance.

PyHydroGeophysX.inversion.ert_inversion.MESH_SUFFIXES = ('.bms', '.msh', '.vtk', '.vtu', '.poly')

Mesh formats a user can hand the inversion. .bms is PyGIMLi’s own, .msh is Gmsh (the usual route for a complex 3D domain), and the rest are what PyGIMLi’s loader recognises.

class PyHydroGeophysX.inversion.ert_inversion.ModelResult(mesh, model, response=None, coverage=None, velocity=None)

Bases: object

Manager-shaped view of an inverted model, so both engines feed one viewer.

MeshResultView and the VTK export ask for paraDomain, model and an optional coverage(); the in-house engine returns arrays rather than a PyGIMLi manager, so this wraps them in the same shape.

velocity is set only by travel time, where a PyGIMLi manager exposes it under that name. It stays None for ERT, because a velocity attribute that quietly returned resistivity would be worse than a missing one.

coverage()

PyHydroGeophysX.inversion.ert_inversion.ensure_geometric_factors(container, mesh, *, policy: str = 'fix', tolerance: float = 0.05, rho0: float = 100.0, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Validate k before the data are used, and repair it when asked.

A repair is kept only if it actually makes the check pass; otherwise the original factors are restored, because a failed repair is worse than a reported problem. Either way the caller is told what happened.

PyHydroGeophysX.inversion.ert_inversion.load_inversion_mesh(mesh_path: str | ~pathlib.Path, data=None, log: ~typing.Callable[[str], None] = <function noop>)

Load a user-supplied inversion mesh and check it can hold this survey.

Building a mesh from the electrode line is fine for a 2D profile and hopeless for a 3D domain with topography, boreholes or known structure, so those are meshed externally (usually in Gmsh) and brought in here.

An imported mesh fails in ways a generated one cannot: electrodes outside the domain, or every cell marked background so nothing is inverted. Both surface deep inside the forward solver as errors that name nothing useful, so they are checked here where the message can say what is wrong.

PyHydroGeophysX.inversion.ert_inversion.repair_geometric_factors(container, mesh, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Recompute k numerically on the inversion mesh and rebuild rhoa from R.

The transfer resistance R = rhoa/k is what the instrument actually measured; k and rhoa are both derived from it. So the repair keeps R fixed, takes k from a forward run on the mesh the inversion will use, and rebuilds rhoa = R * k. Numerical factors also carry the topography that the analytic half-space formula cannot.

R has to be recovered with the factors that formed rhoa, which is the file’s own k_file when the loader kept it, and only otherwise the container’s k. Dividing by the wrong one would bake the discrepancy into R and leave the section scaled after the repair rather than before it.

Modifies container in place and returns what changed.

PyHydroGeophysX.inversion.ert_inversion.run_ert_manager_inversion(data_path: str | ~pathlib.Path, output_dir: str | ~pathlib.Path, *, relative_error: float = 0.03, absolute_error: float = 0.0, error_source: str = 'file', error_floor: float = 0.005, mesh_quality: float = 34.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, mesh_file: str = '', lam: float = 50.0, max_iterations: int = 20, plateau_tolerance: float = 0.005, max_total_iterations: int = 60, engine: str = 'pyhydro', model_constraints: ~typing.Tuple[float, float] = (0.01, 100000.0), solver: str = 'cgls', geometric_factor_policy: str = 'fix', geometric_factor_tolerance: float = 0.05, instrument: str | None = None, reject_outliers: bool = False, outlier_threshold: float = 3.0, outlier_passes: int = 2, min_data_fraction: float = 0.5, auto_lambda: bool = False, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_lambda_trials: int = 6, lambda_bounds: ~typing.Tuple[float, float] = (0.001, 100000.0), lambda_warm_start: bool = True, lambda_cold_retry_chi2: float = 15.0, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert one ERT dataset, in the order that actually lowers chi2.

The stages run back to back, each optional:

0. Geometric factors (geometric_factor_policy). A homogeneous forward run must return the model resistivity; if it does not, k disagrees with the geometry being modelled and the whole section is scaled by that factor with no trace in chi2. fix recomputes k numerically on the inversion mesh and rebuilds rhoa from the measured transfer resistance.

1. Error model. error_source decides whether the file’s own err column is trusted, recomputed from relative_error/absolute_error, or combined. Overwriting a measured error with an assumed one makes chi2 report on an error model the data never had.

2. Inversion at the requested lambda, iterated to a plateau. A run that exhausts max_iterations is continued (up to max_total_iterations) rather than being judged where it stopped. This run is always kept, under fixed_lambda.

3. Outlier rejection (reject_outliers). Measurements the converged model cannot explain are dropped and the inversion repeated, at most outlier_passes times and never below min_data_fraction of the data.

4. Lambda search (auto_lambda). Only now, with a plateaued misfit on the cleaned data, is lambda allowed to move; every trial is itself iterated to a plateau before its chi2 counts. With lambda_warm_start each trial continues from the nearest lambda already solved rather than restarting from a homogeneous model, which is why lam should start on the smooth side: the sweep is then a relaxation from an over-regularized model down to the data, the direction in which continuation is stable.

engine selects the solver: "pyhydro" for the in-house Gauss-Newton inversion (ERTInversion), "pygimli" for ert.ERTManager, or "adtlert" for the optional differentiable ADTLERT 2.5D backend when Torch, CuPy CUDA 12 and cuDSS are available. Linux and Windows both use cuDSS; ADTLERT’s slower SciPy forward solver is intentionally disabled. Linux remains the recommended platform for the best performance. Without CUDA 12 or cuDSS, "adtlert" falls back to "pyhydro".

ADTLERT matches the published real-data branch: fast normal sensitivity, no Robin-boundary derivative, line search, a maximum log step of one, and GPU CGLS when CUDA is available.

PyHydroGeophysX.inversion.ert_inversion.validate_geometric_factors(container, mesh, *, rho0: float = 100.0, tolerance: float = 0.05, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Check k by forward-modelling a homogeneous half space.

A homogeneous model of resistivity rho0 must return rho0 as the apparent resistivity of every configuration. The forward response is (U/I) * k, so the returned ratio is exactly k / k_true for the geometry actually being modelled, and it is independent of the field data.

This catches what chi2 cannot. A geometric factor that is uniformly wrong by a factor X forces the inversion to scale the model by 1/X to fit the same apparent resistivities, so the section is wrong by X while the fit looks perfect. Topography, wrong electrode spacing, and half-space versus full-space convention errors all show up here.

Returns the ratio statistics, an ok flag, and a message naming the likely cause and the consequence.

PyHydroGeophysX.inversion.fdem_inversion module

Inversion utilities for Frequency-Domain Electromagnetic (FDEM) data.

Provides 1D layered-earth FDEM inversion using SimPEG.

class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversion(frequencies: ndarray, dobs: ndarray, uncertainties: ndarray, thicknesses: ndarray | None = None, source_location: ndarray | None = None, source_radius: float = 10.0, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', waveform_type: str = 'dipole', n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

1D FDEM inversion using SimPEG.

Follows the same pattern used in TDEMInversion.

run(starting_model: ndarray | None = None) → FDEMInversionResult

setup() → None

class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: Any = None, thicknesses: ndarray = None, chi2: float = None, frequencies: ndarray = None)

Bases: object

Container for FDEM inversion results.

chi2: float = None

frequencies: ndarray = None

l2_conductivity: ndarray = None

l2_model: ndarray = None

mesh: Any = None

predicted_data: ndarray = None

recovered_conductivity: ndarray = None

recovered_model: ndarray = None

thicknesses: ndarray = None

PyHydroGeophysX.inversion.gravmag module

SimPEG gravity and magnetics inversion.

PyHydroGeophysX.inversion.gravmag.BETA_BOUNDS: Tuple[float, float] = (1e-12, 100000000.0)

Beta is not lambda. It scales against the data units and the mesh, so its useful range sits far below the ERT/SRT regularization bounds; on a typical gravity survey chi2 = 1 lands near 1e-3. The estimate below makes the search scale free, and these bounds only stop it running away.

exception PyHydroGeophysX.inversion.gravmag.InversionBackendUnavailable

Bases: BackendUnavailable

SimPEG / discretize / a usable solver could not be imported.

PyHydroGeophysX.inversion.gravmag.apply_sensitivity_weights(sim, dmis, reg, m: ~numpy.ndarray, *, power: float = 1.0, floor: float = 1e-08, log: ~typing.Callable[[str], None] = <function noop>) → ndarray | None

Weight the regularization by each cell’s sensitivity, and say if it worked.

Potential-field sensitivity falls off sharply with depth, so an unweighted smallness term buys its misfit reduction most cheaply at the surface: the recovered body ends up plastered against the top of the mesh whatever its real depth. Weighting the model term by sqrt(diag(J^T W^T W J)) is the standard correction (Li and Oldenburg, 1996) and is what SimPEG’s UpdateSensitivityWeights directive applies on the iterative path.

Measured on a 0.7 g/cc block buried 60-160 m: unweighted the peak sits at 12 m, weighted it sits at 88 m.

PyHydroGeophysX.inversion.gravmag.backend_status() → Dict[str, Any]

Report whether the SimPEG potential-field inversion stack is available.

PyHydroGeophysX.inversion.gravmag.estimate_beta0(dmis, reg, m, *, ratio: float = 1.0, seed: int = 42, n_power: int = 20) → float

Scale-free starting beta, as SimPEG’s BetaEstimate_ByEig does it.

Power-iterate both Hessians and take the ratio of their largest eigenvalues, so beta starts where the two objective terms are comparable regardless of the data units or the mesh size.

PyHydroGeophysX.inversion.gravmag.invert_gravmag(x, y, value, kind: str, *, z: ~numpy.ndarray | None = None, field: ~typing.Dict[str, ~typing.Any] | None = None, detrend: int = 0, n_xy: int = 22, n_z: int = 12, max_iterations: int = 20, beta0_ratio: float = 1.0, max_stations: int = 600, relative_error: float = 0.03, noise_floor: float | None = None, solver: str = 'simpeg', auto_beta: bool = True, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_beta_trials: int = 6, sensitivity_power: float = 1.0, out_dir: str | None = None, random_seed: int | None = 42, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Run a SimPEG 3D potential-field inversion under the survey.

gravity recovers a density-contrast model (g/cc); magnetics recovers a susceptibility model (SI) and needs field = {inclination, declination, strength_nT}. z is optional per-station elevation (m, positive upward); a missing value falls back to 1 m. detrend (0..3) removes a polynomial regional trend before inversion. The returned grid uses elevation increasing upward. random_seed makes SimPEG’s eigenvalue-based beta estimate reproducible. Raises InversionBackendUnavailable if SimPEG is missing.

PyHydroGeophysX.inversion.gravmag.solve_tikhonov(dmis, reg, m_ref: ndarray, beta: float, *, bounds: Tuple[float, float] | None = None, cg_maxiter: int = 400, cg_tol: float = 1e-08)

Minimize phi_d(m) + beta * phi_m(m) for a linear forward operator.

Potential-field sensitivities do not depend on the model, so the objective is quadratic and a single Newton step from any point is the exact minimizer; there is no Gauss-Newton loop and no line search to run. The step is taken by conjugate gradients on SimPEG’s own deriv/deriv2, which keeps every weighting convention identical to the directive-driven path.

bounds are applied by clipping. That is a projection onto the box, not a constrained optimum, so the returned model is only the true minimizer when it lies inside; clipped in the result says whether the bound bit.

PyHydroGeophysX.inversion.gravmag.sweep_beta_for_chi2(dmis, reg, m_ref: ~numpy.ndarray, n_data: int, *, beta0: float, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_trials: int = 6, bounds: ~typing.Tuple[float, float] | None = None, beta_bounds: ~typing.Tuple[float, float] = (1e-12, 100000000.0), cg_maxiter: int = 400, cg_tol: float = 1e-08, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Find the beta whose chi2 lands on target_chi2.

Each trial is one linear solve rather than a whole nonlinear inversion, so unlike the ERT and travel-time searches this one can afford to be thorough, and there is nothing to warm start: the minimizer for a given beta does not depend on where the solve began.

PyHydroGeophysX.inversion.joint module

Capability registry and public dispatch API for multi-method inversion.

class PyHydroGeophysX.inversion.joint.JointInversionRequest(method_a: str, method_b: str, strategy: str, data: ~typing.Dict[str, ~typing.Any], parameters: ~typing.Dict[str, ~typing.Any] = <factory>, output_dir: str | ~pathlib.Path = 'results/joint_inversion', run_baseline: bool = True)

Bases: object

Input contract for a registered joint inversion runner.

data: Dict[str, Any]

method_a: str

method_b: str

output_dir: str | Path = 'results/joint_inversion'

parameters: Dict[str, Any]

run_baseline: bool = True

strategy: str

class PyHydroGeophysX.inversion.joint.JointInversionResult(methods: ~typing.Tuple[str, str], strategy: str, models: ~typing.Dict[str, ~typing.Any] = <factory>, predicted: ~typing.Dict[str, ~typing.Any] = <factory>, coverage: ~typing.Dict[str, ~typing.Any] = <factory>, chi2: ~typing.Dict[str, float] = <factory>, history: ~typing.List[~typing.Dict[str, ~typing.Any]] = <factory>, baseline: ~typing.Dict[str, ~typing.Any] = <factory>, artifacts: ~typing.Dict[str, str] = <factory>, warnings: ~typing.List[str] = <factory>, meta: ~typing.Dict[str, ~typing.Any] = <factory>, status: str = 'success')

Bases: object

Method-neutral result returned by all registered joint runners.

artifacts: Dict[str, str]

baseline: Dict[str, Any]

chi2: Dict[str, float]

coverage: Dict[str, Any]

history: List[Dict[str, Any]]

meta: Dict[str, Any]

methods: Tuple[str, str]

models: Dict[str, Any]

predicted: Dict[str, Any]

status: str = 'success'

strategy: str

summary() → Dict[str, Any]

Return a compact JSON-ready representation without model arrays.

warnings: List[str]

class PyHydroGeophysX.inversion.joint.JointPairCapability(methods: Tuple[str, str], strategies: Mapping[str, str], dimension: str, model_parameter: str, implemented: bool, dependencies: Tuple[str, ...] = (), description: str = '', runner: str | None = None, backends: Tuple[str, ...] = ())

Bases: object

Describe the strategies available for one normalized method pair.

backends: Tuple[str, ...] = ()

dependencies: Tuple[str, ...] = ()

description: str = ''

dimension: str

implemented: bool

methods: Tuple[str, str]

model_parameter: str

runner: str | None = None

strategies: Mapping[str, str]

PyHydroGeophysX.inversion.joint.get_joint_capabilities(include_planned: bool = True) → List[JointPairCapability]

List implemented capabilities and, optionally, planned method pairs.

PyHydroGeophysX.inversion.joint.get_joint_capability(method_a: str, method_b: str) → JointPairCapability

Return one capability, including a planned placeholder if unsupported.

PyHydroGeophysX.inversion.joint.normalize_joint_pair(method_a: str, method_b: str) → Tuple[str, str]

Return a stable pair key and reject duplicate methods.

PyHydroGeophysX.inversion.joint.run_joint_inversion(request: JointInversionRequest | Mapping[str, Any], progress: Any | None = None) → JointInversionResult

Validate and execute a registered joint or cooperative inversion.

PyHydroGeophysX.inversion.joint_api module

Dependency-free public types and capability registry for joint inversion.

class PyHydroGeophysX.inversion.joint_api.JointInversionRequest(method_a: str, method_b: str, strategy: str, data: ~typing.Dict[str, ~typing.Any], parameters: ~typing.Dict[str, ~typing.Any] = <factory>, output_dir: str | ~pathlib.Path = 'results/joint_inversion', run_baseline: bool = True)

Bases: object

Input contract for a registered joint inversion runner.

data: Dict[str, Any]

method_a: str

method_b: str

output_dir: str | Path = 'results/joint_inversion'

parameters: Dict[str, Any]

run_baseline: bool = True

strategy: str

class PyHydroGeophysX.inversion.joint_api.JointInversionResult(methods: ~typing.Tuple[str, str], strategy: str, models: ~typing.Dict[str, ~typing.Any] = <factory>, predicted: ~typing.Dict[str, ~typing.Any] = <factory>, coverage: ~typing.Dict[str, ~typing.Any] = <factory>, chi2: ~typing.Dict[str, float] = <factory>, history: ~typing.List[~typing.Dict[str, ~typing.Any]] = <factory>, baseline: ~typing.Dict[str, ~typing.Any] = <factory>, artifacts: ~typing.Dict[str, str] = <factory>, warnings: ~typing.List[str] = <factory>, meta: ~typing.Dict[str, ~typing.Any] = <factory>, status: str = 'success')

Bases: object

Method-neutral result returned by all registered joint runners.

artifacts: Dict[str, str]

baseline: Dict[str, Any]

chi2: Dict[str, float]

coverage: Dict[str, Any]

history: List[Dict[str, Any]]

meta: Dict[str, Any]

methods: Tuple[str, str]

models: Dict[str, Any]

predicted: Dict[str, Any]

status: str = 'success'

strategy: str

summary() → Dict[str, Any]

Return a compact JSON-ready representation without model arrays.

warnings: List[str]

class PyHydroGeophysX.inversion.joint_api.JointMethodAdapter(*args, **kwargs)

Bases: Protocol

Interface implemented by a capability-specific request runner.

capability: JointPairCapability

run(request: JointInversionRequest) → JointInversionResult

Execute a validated request and return method-neutral results.

class PyHydroGeophysX.inversion.joint_api.JointPairCapability(methods: Tuple[str, str], strategies: Mapping[str, str], dimension: str, model_parameter: str, implemented: bool, dependencies: Tuple[str, ...] = (), description: str = '', runner: str | None = None, backends: Tuple[str, ...] = ())

Bases: object

Describe the strategies available for one normalized method pair.

backends: Tuple[str, ...] = ()

dependencies: Tuple[str, ...] = ()

description: str = ''

dimension: str

implemented: bool

methods: Tuple[str, str]

model_parameter: str

runner: str | None = None

strategies: Mapping[str, str]

PyHydroGeophysX.inversion.joint_api.get_joint_capabilities(include_planned: bool = True) → List[JointPairCapability]

List implemented capabilities and, optionally, planned method pairs.

PyHydroGeophysX.inversion.joint_api.get_joint_capability(method_a: str, method_b: str) → JointPairCapability

Return one capability, including a planned placeholder if unsupported.

PyHydroGeophysX.inversion.joint_api.normalize_joint_pair(method_a: str, method_b: str) → Tuple[str, str]

Return a stable pair key and reject duplicate methods.

PyHydroGeophysX.inversion.joint_api.pair_joint_soundings(f_value: Any, t_value: Any, parameters: Mapping[str, Any]) → List[Tuple[int, int, str]]

Pair FDEM/TDEM soundings by coordinate or, when safe, by index.

PyHydroGeophysX.inversion.joint_api.split_joint_soundings(value: Any) → Tuple[List[Mapping[str, Any]], ndarray | None]

Normalize one EM sounding or a line-of-soundings container.

PyHydroGeophysX.inversion.joint_api.validate_profile_interface(interface: Any, sensor_x: Any, sensor_z: Any) → Tuple[ndarray, ndarray]

Validate a 2-D structural interface against an ERT surface profile.

PyHydroGeophysX.inversion.joint_ert_srt module

Joint ERT-SRT inversion with structural and geostatistical constraints.

This module integrates the legacy alternating inversion strategy from cg_joint_inversion_hang into the package architecture.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTInversion(ert_data: str | PathLike | pygimli.DataContainer, srt_data: str | PathLike | pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Joint ERT-SRT inversion using alternating Gauss-Newton updates.

• ERT model parameter: log-resistivity

• SRT model parameter: log-slowness

• Coupling: linearized cross-gradient terms (B1, B2)

• Regularization: smoothness or geostatistical

run() → JointERTSRTResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTResult(ert_log_resistivity: ndarray | None = None, srt_log_slowness: ndarray | None = None, ert_resistivity: ndarray | None = None, srt_velocity: ndarray | None = None, ert_predicted: ndarray | None = None, srt_predicted: ndarray | None = None, ert_coverage: ndarray | None = None, srt_coverage: ndarray | None = None, chi2_ert: float | None = None, chi2_srt: float | None = None, iteration_history: list = <factory>, mesh: Any = None, meta: Dict[str, ~typing.Any]=<factory>)

Bases: object

Container for joint ERT-SRT inversion outputs.

chi2_ert: float | None = None

chi2_srt: float | None = None

ert_coverage: ndarray | None = None

ert_log_resistivity: ndarray | None = None

ert_predicted: ndarray | None = None

ert_resistivity: ndarray | None = None

iteration_history: list

mesh: Any = None

meta: Dict[str, Any]

srt_coverage: ndarray | None = None

srt_log_slowness: ndarray | None = None

srt_predicted: ndarray | None = None

srt_velocity: ndarray | None = None

PyHydroGeophysX.inversion.joint_fdem_tdem module

Shared-model joint inversion for collocated FDEM and TDEM soundings.

class PyHydroGeophysX.inversion.joint_fdem_tdem.JointFDEMTDEMInversion(fdem_data: Mapping[str, Any], tdem_data: Mapping[str, Any], fdem_geometry: Mapping[str, Any] | None = None, tdem_geometry: Mapping[str, Any] | None = None, thicknesses: ndarray | None = None, **kwargs: Any)

Bases: object

Invert one FDEM and one TDEM sounding for a common 1-D conductivity model.

Both datasets retain their own survey geometry and uncertainty model. Their normalized residuals are combined with a common first-difference smoothness term, so neither method dominates solely because it has more channels.

classmethod native_backend_status() → Dict[str, Any]

Report whether the tested native SimPEG joint path is available.

run() → JointFDEMTDEMResult

Run native SimPEG joint inversion or the compatible SciPy fallback.

class PyHydroGeophysX.inversion.joint_fdem_tdem.JointFDEMTDEMResult(resistivity: ndarray, conductivity: ndarray, thicknesses: ndarray, predicted_fdem: ndarray, predicted_tdem: ndarray, chi2_fdem: float, chi2_tdem: float, coverage_fdem: ndarray | None = None, coverage_tdem: ndarray | None = None, convergence: Dict[str, float]]=<factory>, baseline: Dict[str, ~typing.Any]=<factory>, meta: Dict[str, ~typing.Any]=<factory>)

Bases: object

Outputs from a shared-conductivity FDEM-TDEM inversion.

baseline: Dict[str, Any]

chi2_fdem: float

chi2_tdem: float

conductivity: ndarray

convergence: List[Dict[str, float]]

coverage_fdem: ndarray | None = None

coverage_tdem: ndarray | None = None

meta: Dict[str, Any]

predicted_fdem: ndarray

predicted_tdem: ndarray

resistivity: ndarray

thicknesses: ndarray

PyHydroGeophysX.inversion.joint_gravity_magnetics module

SimPEG cross-gradient joint inversion for gravity and magnetic data.

class PyHydroGeophysX.inversion.joint_gravity_magnetics.GravityMagneticsJointResult(density: ~numpy.ndarray, susceptibility: ~numpy.ndarray, predicted_gravity: ~numpy.ndarray, predicted_magnetics: ~numpy.ndarray, coverage_gravity: ~numpy.ndarray, coverage_magnetics: ~numpy.ndarray, chi2_gravity: float, chi2_magnetics: float, convergence: ~typing.List[~typing.Dict[str, ~typing.Any]], edges: ~typing.Tuple[~numpy.ndarray, ~numpy.ndarray, ~numpy.ndarray], model_shape: ~typing.Tuple[int, int, int], cross_gradient: ~numpy.ndarray, baseline: ~typing.Dict[str, ~typing.Any] = <factory>, meta: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: object

Result of a shared-mesh gravity–magnetics inversion.

baseline: Dict[str, Any]

chi2_gravity: float

chi2_magnetics: float

convergence: List[Dict[str, Any]]

coverage_gravity: ndarray

coverage_magnetics: ndarray

cross_gradient: ndarray

density: ndarray

edges: Tuple[ndarray, ndarray, ndarray]

meta: Dict[str, Any]

model_shape: Tuple[int, int, int]

predicted_gravity: ndarray

predicted_magnetics: ndarray

susceptibility: ndarray

class PyHydroGeophysX.inversion.joint_gravity_magnetics.JointGravityMagneticsInversion(gravity_data: Mapping[str, Any], magnetics_data: Mapping[str, Any], *, field: Mapping[str, Any] | None = None, n_xy: int = 12, n_z: int = 8, max_iterations: int = 10, max_stations: int = 600, gravity_relative_error: float = 0.03, magnetics_relative_error: float = 0.03, gravity_noise_floor: float = 0.5, magnetics_noise_floor: float = 2.0, gravity_weight: float = 1.0, magnetics_weight: float = 1.0, cross_gradient_weight: float = 2000000000000.0, beta0_ratio: float = 1.0, gravity_detrend: int = 0, magnetics_detrend: int = 0, run_baseline: bool = True, baseline_max_iterations: int | None = None, random_seed: int | None = 42, output_dir: str | None = None, progress_callback: Callable[[Dict[str, Any]], None] | None = None)

Bases: object

Jointly recover density and susceptibility with SimPEG cross-gradient coupling.

Inputs are mappings with x, y, value and optional z arrays. Gravity values use mGal and magnetic total-field values use nT. Both surveys are placed on one tensor mesh; the stacked inversion model is split by a simpeg.maps.Wires map into density (g/cc) and susceptibility (SI).

run() → GravityMagneticsJointResult

Run the native SimPEG similarity-measure inversion.

PyHydroGeophysX.inversion.lambda_search module

Chi-squared targeted regularization search, shared by every method.

ERT, SRT, time-lapse, potential fields, and the EM1D line inversion all face the same question: the fixed regularization weight the user asked for did not land the misfit near its target, so which weight would? The search below answers it without knowing anything about the physics, so it lives here rather than in any one method’s module. Keeping it free of pygimli and SimPEG imports is the point: the SimPEG-only paths must not pull in pygimli to adjust a scalar.

PyHydroGeophysX.inversion.lambda_search.search_lambda_for_chi2(evaluate: ~typing.Callable[[float], float], *, start_lambda: float, start_chi2: float, target_chi2: float = 1.0, tolerance: float = 0.2, max_trials: int = 6, bounds: ~typing.Tuple[float, float] = (0.001, 100000.0), log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Search for the lambda whose chi-squared lands closest to target_chi2.

evaluate(lam) runs one inversion and returns its chi-squared; it is called at most max_trials times, on top of the run that produced start_chi2.

The search assumes chi2 grows with lambda, because more smoothing fits the data less well. It brackets the target by multiplying or dividing lambda by _BRACKET_FACTOR, then bisects in log-lambda. Monotonicity only steers the next guess: the reported lambda is whichever trial came closest to the target, so a noisy or non-monotonic response degrades to “best of the trials” instead of failing.

Returns a dict with lam, chi2, trials (every visited pair), status, and reason.

PyHydroGeophysX.inversion.metrics module

Backend-neutral helpers for inversion result summaries.

PyHydroGeophysX.inversion.metrics.metrics_from_manager(mgr: Any, *, n_data: int | None = None, lam: float | None = None, method: str | None = None) → Tuple[Dict[str, Any], List[float]]

Extract common metrics from a PyGIMLi-style manager.

PyHydroGeophysX.inversion.multi_method module

Unified multi-method geophysical inversion interface.

class PyHydroGeophysX.inversion.multi_method.GeophysicalInversion(method: str, **kwargs)

Bases: object

Unified factory for multi-method geophysical inversion.

Dispatches to the correct inversion engine.

SUPPORTED = {'ert', 'fdem', 'joint', 'joint_ert_srt', 'srt', 'tdem'}

property engine

run(**kwargs)

setup(**kwargs)

PyHydroGeophysX.inversion.robust_errors module

Bounded Huber IRLS for EM data only; regularisation is never reweighted.

These are effective fitting uncertainties, not revised measurement errors. Every imported gate remains present. Recomputing each factor from the original uncertainty permits a gate to regain weight and prevents cumulative inflation.

PyHydroGeophysX.inversion.robust_errors.huber_error_factor(residual, threshold: float, max_error_factor: float)

sigma_eff/sigma_base, with inverse-variance weight min(1, k/|r|).

The maximum error factor bounds the loss of influence. Beyond that cap this is a bounded Huber-style reweighting, not an exact unbounded Huber loss.

PyHydroGeophysX.inversion.robust_errors.reweight_errors(observed, uncertainty, solve: Callable, predict: Callable, *, threshold: float = 3.0, passes: int = 3, max_error_factor: float = 10.0, min_unchanged_fraction: float = 0.0, target_chi2: float = 0.0, target_tolerance: float = 0.2, solver_ready: Callable = <function <lambda>>, stage_statistics: Callable = <function <lambda>>, history: Callable = <function <lambda>>, log: Callable = <function <lambda>>)

Solve, update uncertainties, and warm-start without ever masking a gate.

solve(sigma, previous) is responsible for freezing regularisation after the initial solve. The returned report always scores the final model against ALL original observations, both with original and with actually used errors. solver_ready can veto further error inflation after an incomplete or failed inner solve. Budget exhaustion is not evidence of contaminated data.

PyHydroGeophysX.inversion.robust_errors.robust_error_options(inv: Dict[str, Any]) → Dict[str, Any]

PyHydroGeophysX.inversion.robust_errors.select_error_factors(residual, threshold=3.0, max_error_factor=10.0, *, min_unchanged_fraction=0.0, target_chi2=0.0)

Limit changes to the worst residuals; optionally calibrate their errors.

This target mode is NOT a Huber likelihood or independently estimated noise. For the current model it searches a bounded error multiplier, then the caller must refit. Protected gates have factor exactly 1, including when the target is unreachable. Ranks may change between passes; ties have stable ordering. The quota is across the input vector (one LCI run, or one independent fit).

PyHydroGeophysX.inversion.srt_inversion module

Seismic Refraction Tomography (SRT) inversion functionality.

Uses PyGIMLi’s TravelTimeManager for forward modeling and Jacobian computation. Provides custom Gauss-Newton inversion with the same architecture as ERTInversion but for travel-time data.

class PyHydroGeophysX.inversion.srt_inversion.SRTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Seismic Refraction Tomography inversion class.

Inverts travel-time data for subsurface velocity via log-slowness, and uses a Gauss-Newton optimization loop.

run(initial_model: ndarray | None = None) → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

PyHydroGeophysX.inversion.srt_inversion.build_srt_mesh(data, *, mesh_quality: float = 32.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, log: ~typing.Callable[[str], None] = <function noop>)

Build the travel-time inversion mesh, mirroring the ERT pipeline.

PyGIMLi sizes the parameter domain from the array length when paraDepth is left at 0. For refraction that reaches well past where any ray turns, so the deep cells are unconstrained and only slow the run down; capping the depth removes them. Cell size and quality trade resolution against cost the same way they do for ERT.

PyHydroGeophysX.inversion.srt_inversion.run_srt_manager_inversion(travel_time_path: str | ~pathlib.Path, output_dir: str | ~pathlib.Path, *, engine: str = 'pygimli', lam: float = 50.0, max_iterations: int = 20, plateau_tolerance: float = 0.005, max_total_iterations: int = 60, mesh_quality: float = 32.0, para_depth: float = 0.0, para_max_cell_size: float = 0.0, secondary_nodes: int = 3, auto_lambda: bool = False, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_lambda_trials: int = 6, lambda_warm_start: bool = True, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert one travel-time dataset for velocity.

With auto_lambda the same machinery as the ERT pipeline applies: each lambda is iterated to a plateau before its chi2 counts, the sweep relaxes from the requested lambda downward, and each trial continues from the nearest lambda already solved. engine="pyhydro" uses the in-house Gauss-Newton solver, which is what the search can drive; "pygimli" runs TravelTimeManager once and is the historical default.

The mesh is built once here and handed to whichever engine runs, so the two invert the same domain. para_depth and para_max_cell_size take 0 to mean “let PyGIMLi size it from the array”; secondary_nodes refines the ray tracing without adding unknowns.

PyHydroGeophysX.inversion.srt_time_lapse module

Time-lapse Seismic Refraction Tomography (SRT) inversion functionality.

Jointly inverts multiple travel-time datasets with spatial and temporal regularization using log-slowness parameterization.

class PyHydroGeophysX.inversion.srt_time_lapse.TimeLapseSRTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Time-lapse Seismic Refraction Tomography inversion.

Architecture mirrors TimeLapseERTInversion while using TravelTimeManager and log-slowness model parameters.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

PyHydroGeophysX.inversion.tdem_inversion module

Inversion utilities for Time-Domain Electromagnetic (TDEM) data.

This module provides classes for 1D TDEM inversion using SimPEG, with support for both L2 and sparse (IRLS) regularization.

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

Class for 1D TDEM inversion using SimPEG.

This class provides functionality for inverting TDEM sounding data to recover 1D layered Earth conductivity models.

Example

>>> # Load or create data
>>> times = np.logspace(-5, -2, 31)
>>> dobs = ...  # observed data
>>> uncertainties = ...  # data uncertainties
>>>
>>> # Create inversion
>>> inv = TDEMInversion(
...     times=times,
...     dobs=dobs,
...     uncertainties=uncertainties,
...     source_radius=10.0
... )
>>>
>>> # Run inversion
>>> result = inv.run()

plot_result(result: TDEMInversionResult, true_model: Tuple[ndarray, ndarray] | None = None, save_path: str | None = None) → None

Plot inversion results.

Parameters:

• result – TDEMInversionResult from run()

• true_model – Tuple of (thicknesses, conductivity) for true model

• save_path – Path to save figure (optional)

run(starting_model: ndarray | None = None) → TDEMInversionResult

Run TDEM inversion.

Parameters:

starting_model – Initial log-conductivity model (optional)

Returns:

TDEMInversionResult with inversion results

setup() → None

Set up inversion components (survey, mesh, simulation).

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: discretize.TensorMesh = None, thicknesses: ndarray = None, chi2: float = None, iterations: int = 0, convergence_history: List[float] = <factory>)

Bases: object

Container for TDEM inversion results.

recovered_model

Final recovered model (log-conductivity)

Type:

numpy.ndarray

recovered_conductivity

Recovered conductivity (S/m)

Type:

numpy.ndarray

l2_model

L2 model before IRLS (log-conductivity)

Type:

numpy.ndarray

l2_conductivity

L2 conductivity (S/m)

Type:

numpy.ndarray

predicted_data

Predicted data from recovered model

Type:

numpy.ndarray

mesh

TensorMesh used for inversion

Type:

discretize.TensorMesh

thicknesses

Layer thicknesses used in inversion

Type:

numpy.ndarray

chi2

Final chi-squared misfit

Type:

float

iterations

Number of iterations

Type:

int

convergence_history

History of data misfit per iteration

Type:

List[float]

chi2: float = None

convergence_history: List[float]

iterations: int = 0

l2_conductivity: ndarray = None

l2_model: ndarray = None

mesh: discretize.TensorMesh = None

predicted_data: ndarray = None

recovered_conductivity: ndarray = None

recovered_model: ndarray = None

thicknesses: ndarray = None

PyHydroGeophysX.inversion.tdem_inversion.run_tdem_inversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, n_layers: int = 25, use_irls: bool = True, verbose: bool = True, **kwargs) → TDEMInversionResult

Convenience function to run TDEM inversion.

Parameters:

• times – Time channels (s)

• dobs – Observed data (T)

• uncertainties – Data uncertainties (T)

• source_radius – Loop radius (m)

• n_layers – Number of inversion layers

• use_irls – Use IRLS for sparse inversion

• verbose – Print progress

• **kwargs – Additional parameters for TDEMInversion

Returns:

TDEMInversionResult with inversion results

PyHydroGeophysX.inversion.time_lapse module

Time-lapse ERT inversion functionality.

class PyHydroGeophysX.inversion.time_lapse.TimeLapseERTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Time-lapse ERT inversion class.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Run time-lapse ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

TimeLapseInversionResult with inversion results

setup()

Set up time-lapse ERT inversion (load data, create operators, matrices, etc.)

PyHydroGeophysX.inversion.windowed module

Windowed time-lapse ERT inversion for handling large temporal datasets.

class PyHydroGeophysX.inversion.windowed.WindowedTimeLapseERTInversion(data_dir: str, ert_files: List[str], measurement_times: List[float], window_size: int = 3, mesh: pygimli.Mesh | str | None = None, engine: str = 'pyhydro', log: Callable[[str], None] | None = None, **kwargs)

Bases: object

Class for windowed time-lapse ERT inversion to handle large temporal datasets.

run(window_parallel: bool = False, max_window_workers: int | None = None) → TimeLapseInversionResult

Run windowed time-lapse ERT inversion.

Parameters:

• window_parallel – Whether to process windows in parallel

• max_window_workers – Maximum number of parallel workers (None for auto)

Returns:

TimeLapseInversionResult with stitched results

Module contents

Lazy inversion-framework exports.

class PyHydroGeophysX.inversion.ERTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Single-time ERT inversion class.

run(initial_model: ndarray | None = None, reference_model: ndarray | None = None) → InversionResult

Run ERT inversion.

Parameters:

• initial_model – Initial model parameters (if None, a homogeneous model is used)

• reference_model – Model the regularization pulls toward. Defaults to the initial model, which is the usual smoothness-from-homogeneous setup. Pass it explicitly when continuing an interrupted run. Leaving it to default makes the restart point the new reference, so the penalty becomes roughness of the change since the restart rather than roughness of the model; accumulated structure then goes unpenalized and the run drifts steadily under-regularized. Measured on the Ashton line over a 20 to 0.3 lambda ladder, that reached chi2 5.72 where the same ladder with the reference pinned reached 6.38.

Returns:

InversionResult with inversion results

setup()

Set up ERT inversion (create operators, matrices, etc.)

class PyHydroGeophysX.inversion.FDEMInversion(frequencies: ndarray, dobs: ndarray, uncertainties: ndarray, thicknesses: ndarray | None = None, source_location: ndarray | None = None, source_radius: float = 10.0, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', waveform_type: str = 'dipole', n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

1D FDEM inversion using SimPEG.

Follows the same pattern used in TDEMInversion.

run(starting_model: ndarray | None = None) → FDEMInversionResult

setup() → None

class PyHydroGeophysX.inversion.FDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: Any = None, thicknesses: ndarray = None, chi2: float = None, frequencies: ndarray = None)

Bases: object

Container for FDEM inversion results.

chi2: float = None

frequencies: ndarray = None

l2_conductivity: ndarray = None

l2_model: ndarray = None

mesh: Any = None

predicted_data: ndarray = None

recovered_conductivity: ndarray = None

recovered_model: ndarray = None

thicknesses: ndarray = None

class PyHydroGeophysX.inversion.GeophysicalInversion(method: str, **kwargs)

Bases: object

Unified factory for multi-method geophysical inversion.

Dispatches to the correct inversion engine.

SUPPORTED = {'ert', 'fdem', 'joint', 'joint_ert_srt', 'srt', 'tdem'}

property engine

run(**kwargs)

setup(**kwargs)

class PyHydroGeophysX.inversion.GravityMagneticsJointResult(density: ~numpy.ndarray, susceptibility: ~numpy.ndarray, predicted_gravity: ~numpy.ndarray, predicted_magnetics: ~numpy.ndarray, coverage_gravity: ~numpy.ndarray, coverage_magnetics: ~numpy.ndarray, chi2_gravity: float, chi2_magnetics: float, convergence: ~typing.List[~typing.Dict[str, ~typing.Any]], edges: ~typing.Tuple[~numpy.ndarray, ~numpy.ndarray, ~numpy.ndarray], model_shape: ~typing.Tuple[int, int, int], cross_gradient: ~numpy.ndarray, baseline: ~typing.Dict[str, ~typing.Any] = <factory>, meta: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: object

Result of a shared-mesh gravity–magnetics inversion.

baseline: Dict[str, Any]

chi2_gravity: float

chi2_magnetics: float

convergence: List[Dict[str, Any]]

coverage_gravity: ndarray

coverage_magnetics: ndarray

cross_gradient: ndarray

density: ndarray

edges: Tuple[ndarray, ndarray, ndarray]

meta: Dict[str, Any]

model_shape: Tuple[int, int, int]

predicted_gravity: ndarray

predicted_magnetics: ndarray

susceptibility: ndarray

exception PyHydroGeophysX.inversion.InversionBackendUnavailable

Bases: BackendUnavailable

SimPEG / discretize / a usable solver could not be imported.

class PyHydroGeophysX.inversion.InversionBase(data: pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: object

Abstract base class for geophysical inversion methods.

This class provides a foundational structure for specific inversion techniques (e.g., ERT, SRT). It manages observed data, mesh, common inversion parameters, and an InversionResult object. Subclasses must implement methods for setting up the inversion, running the inversion loop, computing Jacobians, and evaluating the objective function.

compute_jacobian(model: ndarray) → ndarray

Abstract method to compute the Jacobian matrix (sensitivity matrix).

The Jacobian J_ij = ∂d_i / ∂m_j relates changes in model parameters (m_j) to changes in observed data (d_i).

Parameters:

model (np.ndarray) – The current model parameter vector for which to compute the Jacobian.

Returns:,

np.ndarray: The computed Jacobian matrix, typically of shape (n_data, n_model_params).

objective_function(model: ndarray, data_to_fit: ndarray | None = None) → float

Abstract method to compute the value of the objective function.

The objective function typically includes a data misfit term and one or more regularization terms: Φ(m) = Φ_d(m) + λ * Φ_m(m).

Parameters:

• model (np.ndarray) – The current model parameter vector.

• data_to_fit (Optional[np.ndarray], optional) – The observed data to fit. If None, self.data (from DataContainer) is typically used. Defaults to None.

Returns:

The calculated value of the objective function.

Return type:

float

run() → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.InversionResult

Bases: object

Base class to store, save, load, and plot results from a geophysical inversion.

final_model

The final inverted model parameters (e.g., resistivity, velocity). Typically a 1D array corresponding to mesh cells.

Type:

Optional[np.ndarray]

predicted_data

The data predicted by the forward model using the final_model.

Type:

Optional[np.ndarray]

coverage

Coverage or sensitivity values for the model parameters, often derived from the Jacobian or resolution matrix.

Type:

Optional[np.ndarray]

mesh

The PyGIMLi mesh object used in the inversion.

Type:

Optional[pg.Mesh]

iteration_models

A list storing the model parameters at each iteration of the inversion.

Type:

List[np.ndarray]

iteration_data_errors

A list storing the data misfit (e.g., residuals) at each iteration.

Type:

List[np.ndarray]

iteration_chi2

A list storing the chi-squared (χ²) value or a similar misfit metric at each iteration.

Type:

List[float]

meta

A dictionary to store any additional metadata about the inversion run (e.g., inversion parameters, timings, comments).

Type:

Dict[str, Any]

classmethod load(filename: str) → InversionResult

Load inversion results from a file previously saved by the save method.

Parameters:

filename (str) – The base path to the saved results file. If saved as name.pkl and name.bms, provide name.pkl or name.

Returns:

An instance of InversionResult (or a subclass if called on one)

populated with the loaded data.

Return type:

InversionResult

Raises:

• FileNotFoundError – If the main data file or associated mesh file (if referenced) is not found.

• IOError – If there’s an error during file reading.

• pickle.UnpicklingError – If the file cannot be unpickled.

plot_convergence(ax: Axes | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the convergence curve (chi-squared misfit vs. iteration number).

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• **kwargs (Any) – Additional keyword arguments passed directly to ax.plot (e.g., color, marker, linestyle).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.iteration_chi2 is empty (no convergence data).

plot_model(ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot the final inverted model on its associated mesh.

Parameters:

• ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

• cmap (str, optional) – The colormap to use for visualizing model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – If provided, cells with coverage values below this threshold will be masked (made transparent or semi-transparent) in the plot. Requires self.coverage to be populated. Defaults to None (no masking).

• **kwargs (Any) – Additional keyword arguments passed directly to pygimli.show (e.g., cMin, cMax, orientation, logScale).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.final_model or self.mesh is None.

save(filename: str) → None

Save the inversion results to a file using Python’s pickle format. If a mesh is present, it is saved separately as a PyGIMLi binary mesh file (.bms).

Parameters:

filename (str) – The base path (including filename without extension) to save the results. The main data will be saved as filename.pkl (or just filename if user includes .pkl). The mesh will be saved as filename.bms or filename.pkl.bms. It’s recommended to provide filename without .pkl.

Raises:

• IOError – If there’s an error during file writing.

• pickle.PicklingError – If an object cannot be pickled.

class PyHydroGeophysX.inversion.JointERTSRTInversion(ert_data: str | PathLike | pygimli.DataContainer, srt_data: str | PathLike | pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Joint ERT-SRT inversion using alternating Gauss-Newton updates.

• ERT model parameter: log-resistivity

• SRT model parameter: log-slowness

• Coupling: linearized cross-gradient terms (B1, B2)

• Regularization: smoothness or geostatistical

RCM: ndarray | None

Wd_ert: scipy.sparse.csr_matrix | None

Wd_srt: scipy.sparse.csr_matrix | None

Wm_ert: scipy.sparse.csr_matrix | None

Wm_srt: scipy.sparse.csr_matrix | None

X: ndarray | None

dobs_ert: ndarray | None

dobs_srt: ndarray | None

ert_grid: pygimli.Mesh | None

mr: ndarray | None

mr_ref: ndarray | None

mv: ndarray | None

mv_ref: ndarray | None

run() → JointERTSRTResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.JointERTSRTResult(ert_log_resistivity: ndarray | None = None, srt_log_slowness: ndarray | None = None, ert_resistivity: ndarray | None = None, srt_velocity: ndarray | None = None, ert_predicted: ndarray | None = None, srt_predicted: ndarray | None = None, ert_coverage: ndarray | None = None, srt_coverage: ndarray | None = None, chi2_ert: float | None = None, chi2_srt: float | None = None, iteration_history: list = <factory>, mesh: Any = None, meta: Dict[str, ~typing.Any]=<factory>)

Bases: object

Container for joint ERT-SRT inversion outputs.

chi2_ert: float | None = None

chi2_srt: float | None = None

ert_coverage: ndarray | None = None

ert_log_resistivity: ndarray | None = None

ert_predicted: ndarray | None = None

ert_resistivity: ndarray | None = None

iteration_history: list

mesh: Any = None

meta: Dict[str, Any]

srt_coverage: ndarray | None = None

srt_log_slowness: ndarray | None = None

srt_predicted: ndarray | None = None

srt_velocity: ndarray | None = None

class PyHydroGeophysX.inversion.JointFDEMTDEMInversion(fdem_data: Mapping[str, Any], tdem_data: Mapping[str, Any], fdem_geometry: Mapping[str, Any] | None = None, tdem_geometry: Mapping[str, Any] | None = None, thicknesses: ndarray | None = None, **kwargs: Any)

Bases: object

Invert one FDEM and one TDEM sounding for a common 1-D conductivity model.

Both datasets retain their own survey geometry and uncertainty model. Their normalized residuals are combined with a common first-difference smoothness term, so neither method dominates solely because it has more channels.

classmethod native_backend_status() → Dict[str, Any]

Report whether the tested native SimPEG joint path is available.

parameters: Dict[str, Any]

run() → JointFDEMTDEMResult

Run native SimPEG joint inversion or the compatible SciPy fallback.

class PyHydroGeophysX.inversion.JointFDEMTDEMResult(resistivity: ndarray, conductivity: ndarray, thicknesses: ndarray, predicted_fdem: ndarray, predicted_tdem: ndarray, chi2_fdem: float, chi2_tdem: float, coverage_fdem: ndarray | None = None, coverage_tdem: ndarray | None = None, convergence: Dict[str, float]]=<factory>, baseline: Dict[str, ~typing.Any]=<factory>, meta: Dict[str, ~typing.Any]=<factory>)

Bases: object

Outputs from a shared-conductivity FDEM-TDEM inversion.

baseline: Dict[str, Any]

chi2_fdem: float

chi2_tdem: float

conductivity: ndarray

convergence: List[Dict[str, float]]

coverage_fdem: ndarray | None = None

coverage_tdem: ndarray | None = None

meta: Dict[str, Any]

predicted_fdem: ndarray

predicted_tdem: ndarray

resistivity: ndarray

thicknesses: ndarray

class PyHydroGeophysX.inversion.JointGravityMagneticsInversion(gravity_data: Mapping[str, Any], magnetics_data: Mapping[str, Any], *, field: Mapping[str, Any] | None = None, n_xy: int = 12, n_z: int = 8, max_iterations: int = 10, max_stations: int = 600, gravity_relative_error: float = 0.03, magnetics_relative_error: float = 0.03, gravity_noise_floor: float = 0.5, magnetics_noise_floor: float = 2.0, gravity_weight: float = 1.0, magnetics_weight: float = 1.0, cross_gradient_weight: float = 2000000000000.0, beta0_ratio: float = 1.0, gravity_detrend: int = 0, magnetics_detrend: int = 0, run_baseline: bool = True, baseline_max_iterations: int | None = None, random_seed: int | None = 42, output_dir: str | None = None, progress_callback: Callable[[Dict[str, Any]], None] | None = None)

Bases: object

Jointly recover density and susceptibility with SimPEG cross-gradient coupling.

Inputs are mappings with x, y, value and optional z arrays. Gravity values use mGal and magnetic total-field values use nT. Both surveys are placed on one tensor mesh; the stacked inversion model is split by a simpeg.maps.Wires map into density (g/cc) and susceptibility (SI).

run() → GravityMagneticsJointResult

Run the native SimPEG similarity-measure inversion.

class PyHydroGeophysX.inversion.JointInversionRequest(method_a: str, method_b: str, strategy: str, data: ~typing.Dict[str, ~typing.Any], parameters: ~typing.Dict[str, ~typing.Any] = <factory>, output_dir: str | ~pathlib.Path = 'results/joint_inversion', run_baseline: bool = True)

Bases: object

Input contract for a registered joint inversion runner.

data: Dict[str, Any]

method_a: str

method_b: str

output_dir: str | Path = 'results/joint_inversion'

parameters: Dict[str, Any]

run_baseline: bool = True

strategy: str

class PyHydroGeophysX.inversion.JointInversionResult(methods: ~typing.Tuple[str, str], strategy: str, models: ~typing.Dict[str, ~typing.Any] = <factory>, predicted: ~typing.Dict[str, ~typing.Any] = <factory>, coverage: ~typing.Dict[str, ~typing.Any] = <factory>, chi2: ~typing.Dict[str, float] = <factory>, history: ~typing.List[~typing.Dict[str, ~typing.Any]] = <factory>, baseline: ~typing.Dict[str, ~typing.Any] = <factory>, artifacts: ~typing.Dict[str, str] = <factory>, warnings: ~typing.List[str] = <factory>, meta: ~typing.Dict[str, ~typing.Any] = <factory>, status: str = 'success')

Bases: object

Method-neutral result returned by all registered joint runners.

artifacts: Dict[str, str]

baseline: Dict[str, Any]

chi2: Dict[str, float]

coverage: Dict[str, Any]

history: List[Dict[str, Any]]

meta: Dict[str, Any]

methods: Tuple[str, str]

models: Dict[str, Any]

predicted: Dict[str, Any]

status: str = 'success'

strategy: str

summary() → Dict[str, Any]

Return a compact JSON-ready representation without model arrays.

warnings: List[str]

class PyHydroGeophysX.inversion.JointPairCapability(methods: Tuple[str, str], strategies: Mapping[str, str], dimension: str, model_parameter: str, implemented: bool, dependencies: Tuple[str, ...] = (), description: str = '', runner: str | None = None, backends: Tuple[str, ...] = ())

Bases: object

Describe the strategies available for one normalized method pair.

backends: Tuple[str, ...] = ()

dependencies: Tuple[str, ...] = ()

description: str = ''

dimension: str

implemented: bool

methods: Tuple[str, str]

model_parameter: str

runner: str | None = None

strategies: Mapping[str, str]

class PyHydroGeophysX.inversion.PetrophysicalCoupling

Bases: object

Coupling helpers from hydrological state to multi-method geophysics.

static compare_inversions_to_hydro(hydro_wc, ert_result, srt_result=None, em_result=None, petro_params=None)

Compare inversion products back to hydrological water content.

Returns misfit statistics for available methods.

static water_content_to_all_geophysics(water_content, porosity, **params)

Convert water content to resistivity, velocity, and conductivity.

Uses existing petrophysical functions in the package.

class PyHydroGeophysX.inversion.SRTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Seismic Refraction Tomography inversion class.

Inverts travel-time data for subsurface velocity via log-slowness, and uses a Gauss-Newton optimization loop.

run(initial_model: ndarray | None = None) → InversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.StructuralConstraint

Bases: object

Build structural constraint terms from one method for another.

static apply_structural_weights_to_Wm(Wm, boundary_weights)

Modify an existing smoothness matrix Wm to respect boundaries.

static build_cross_gradient_operator(mesh, model_a, model_b)

Build a cross-gradient operator surrogate.

Returns a diagonal sparse matrix weighted by local cross-gradient magnitude, where small values indicate better structural agreement.

static build_linearized_cross_gradient_blocks(RCM: ndarray, X: ndarray, model_a: ndarray, model_b: ndarray, mode: str = 'direct') → Tuple[ndarray, ndarray]

Build linearized cross-gradient blocks B1 and B2.

B1 multiplies model_a and B2 multiplies model_b. The resulting penalty terms are ||B1 m_a||^2 and ||B2 m_b||^2.

static build_local_design_matrix(mesh) → ndarray

Build the local design matrix X=[x, z, 1] used by cross-gradient.

This mirrors the legacy xyzpos[:, 2] = 1 construction.

static build_neighborhood_matrix(mesh, Wm: Any | None = None, source: str = 'smoothness', correlation_lengths: Sequence[float] = (4.0, 4.0), threshold: float = 0.01, binarize: bool = True) → ndarray

Build the neighborhood/correlation matrix used by cross-gradient.

Parameters:

• mesh – Mesh object.

• Wm – Smoothness matrix. Required when source='smoothness'.

• source – One of 'smoothness' (legacy direct mode) or 'covariance' (legacy spatial mode).

• correlation_lengths – Correlation lengths passed to pg.utils.covarianceMatrix when source='covariance'.

• threshold – Entries with absolute value below this threshold are zeroed.

• binarize – If True, convert nonzero entries to 1.

static from_conductivity_model(conductivity, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an EM conductivity model.

static from_velocity_model(velocity_model, mesh, gradient_threshold: float = 0.3)

Extract structural boundaries from an SRT velocity model.

Returns a pairwise-boundary weight map that can reduce smoothing across detected structural interfaces.

class PyHydroGeophysX.inversion.TDEMInversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)

Bases: object

Class for 1D TDEM inversion using SimPEG.

This class provides functionality for inverting TDEM sounding data to recover 1D layered Earth conductivity models.

Example

>>> # Load or create data
>>> times = np.logspace(-5, -2, 31)
>>> dobs = ...  # observed data
>>> uncertainties = ...  # data uncertainties
>>>
>>> # Create inversion
>>> inv = TDEMInversion(
...     times=times,
...     dobs=dobs,
...     uncertainties=uncertainties,
...     source_radius=10.0
... )
>>>
>>> # Run inversion
>>> result = inv.run()

plot_result(result: TDEMInversionResult, true_model: Tuple[ndarray, ndarray] | None = None, save_path: str | None = None) → None

Plot inversion results.

Parameters:

• result – TDEMInversionResult from run()

• true_model – Tuple of (thicknesses, conductivity) for true model

• save_path – Path to save figure (optional)

run(starting_model: ndarray | None = None) → TDEMInversionResult

Run TDEM inversion.

Parameters:

starting_model – Initial log-conductivity model (optional)

Returns:

TDEMInversionResult with inversion results

setup() → None

Set up inversion components (survey, mesh, simulation).

class PyHydroGeophysX.inversion.TDEMInversionResult(recovered_model: ndarray = None, recovered_conductivity: ndarray = None, l2_model: ndarray = None, l2_conductivity: ndarray = None, predicted_data: ndarray = None, mesh: discretize.TensorMesh = None, thicknesses: ndarray = None, chi2: float = None, iterations: int = 0, convergence_history: List[float] = <factory>)

Bases: object

Container for TDEM inversion results.

recovered_model

Final recovered model (log-conductivity)

Type:

numpy.ndarray

recovered_conductivity

Recovered conductivity (S/m)

Type:

numpy.ndarray

l2_model

L2 model before IRLS (log-conductivity)

Type:

numpy.ndarray

l2_conductivity

L2 conductivity (S/m)

Type:

numpy.ndarray

predicted_data

Predicted data from recovered model

Type:

numpy.ndarray

mesh

TensorMesh used for inversion

Type:

discretize.TensorMesh

thicknesses

Layer thicknesses used in inversion

Type:

numpy.ndarray

chi2

Final chi-squared misfit

Type:

float

iterations

Number of iterations

Type:

int

convergence_history

History of data misfit per iteration

Type:

List[float]

chi2: float = None

convergence_history: List[float]

iterations: int = 0

l2_conductivity: ndarray = None

l2_model: ndarray = None

mesh: discretize.TensorMesh = None

predicted_data: ndarray = None

recovered_conductivity: ndarray = None

recovered_model: ndarray = None

thicknesses: ndarray = None

class PyHydroGeophysX.inversion.TimeLapseERTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs)

Bases: InversionBase

Time-lapse ERT inversion class.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Run time-lapse ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

TimeLapseInversionResult with inversion results

setup()

Set up time-lapse ERT inversion (load data, create operators, matrices, etc.)

class PyHydroGeophysX.inversion.TimeLapseInversionResult

Bases: InversionResult

Specialized class to store and manage results from time-lapse inversions. Inherits from InversionResult and adds attributes specific to time-lapse data.

final_models

A 2D NumPy array where each column represents the inverted model for a specific timestep (shape: num_cells x num_timesteps).

Type:

Optional[np.ndarray]

timesteps

An array or list of time values (e.g., hours, days) corresponding to each model slice in final_models.

Type:

Optional[np.ndarray]

all_coverage

A list where each element is the coverage array for the corresponding timestep’s model.

Type:

List[np.ndarray]

all_chi2

A list where each element is a list of chi-squared values per iteration for the inversion of that specific timestep or window. (Note: Original was List[float], if chi2 is per window, might need adjustment)

Type:

List[Any]

all_chi2: List[Any]

all_coverage: List[ndarray]

create_time_lapse_animation(output_filename: str, cmap: str = 'viridis', coverage_threshold: float | None = None, dpi: int = 100, fps: int = 2, **kwargs: Any) → None

Create and save an animation (e.g., MP4 video) of the time-lapse inversion results.

Requires ffmpeg or another Matplotlib-supported animation writer to be installed.

Parameters:

• output_filename (str) – The filename for the output animation (e.g., ‘timelapse_animation.mp4’).

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• dpi (int, optional) – Dots Per Inch for the output animation. Defaults to 100.

• fps (int, optional) – Frames Per Second for the animation. Defaults to 2.

• **kwargs (Any) – Additional keyword arguments passed to plot_time_slice for each frame.

Raises:

• ValueError – If final_models or mesh is missing.

• ImportError – If matplotlib.animation cannot be imported.

final_models: ndarray | None

classmethod load(filename: str) → TimeLapseInversionResult

Load time-lapse inversion results. Overrides base class load.

plot_time_slice(timestep_idx: int, ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) → Tuple[Figure, Axes]

Plot a single time slice (inverted model at a specific timestep) from the results.

Parameters:

• timestep_idx (int) – The zero-based index of the timestep to plot.

• ax (Optional[plt.Axes], optional) – Matplotlib Axes to plot on. If None, creates new.

• cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

• coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

• **kwargs (Any) – Additional arguments passed to pg.show.

Returns:

The matplotlib Figure and Axes objects.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If models or mesh are missing, or if timestep_idx is out of range.

save(filename: str) → None

Save time-lapse inversion results. Overrides base class save.

timesteps: ndarray | None

class PyHydroGeophysX.inversion.TimeLapseSRTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs: Any)

Bases: InversionBase

Time-lapse Seismic Refraction Tomography inversion.

Architecture mirrors TimeLapseERTInversion while using TravelTimeManager and log-slowness model parameters.

run(initial_model: ndarray | None = None) → TimeLapseInversionResult

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() → None

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.WindowedTimeLapseERTInversion(data_dir: str, ert_files: List[str], measurement_times: List[float], window_size: int = 3, mesh: pygimli.Mesh | str | None = None, engine: str = 'pyhydro', log: Callable[[str], None] | None = None, **kwargs)

Bases: object

Class for windowed time-lapse ERT inversion to handle large temporal datasets.

run(window_parallel: bool = False, max_window_workers: int | None = None) → TimeLapseInversionResult

Run windowed time-lapse ERT inversion.

Parameters:

• window_parallel – Whether to process windows in parallel

• max_window_workers – Maximum number of parallel workers (None for auto)

Returns:

TimeLapseInversionResult with stitched results

PyHydroGeophysX.inversion.fdem_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert an FDEM sounding for a layered resistivity model (Occam 1D).

PyHydroGeophysX.inversion.get_joint_capabilities(include_planned: bool = True) → List[JointPairCapability]

List implemented capabilities and, optionally, planned method pairs.

PyHydroGeophysX.inversion.invert_gravmag(x, y, value, kind: str, *, z: ~numpy.ndarray | None = None, field: ~typing.Dict[str, ~typing.Any] | None = None, detrend: int = 0, n_xy: int = 22, n_z: int = 12, max_iterations: int = 20, beta0_ratio: float = 1.0, max_stations: int = 600, relative_error: float = 0.03, noise_floor: float | None = None, solver: str = 'simpeg', auto_beta: bool = True, target_chi2: float = 1.0, chi2_tolerance: float = 0.2, max_beta_trials: int = 6, sensitivity_power: float = 1.0, out_dir: str | None = None, random_seed: int | None = 42, log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Run a SimPEG 3D potential-field inversion under the survey.

gravity recovers a density-contrast model (g/cc); magnetics recovers a susceptibility model (SI) and needs field = {inclination, declination, strength_nT}. z is optional per-station elevation (m, positive upward); a missing value falls back to 1 m. detrend (0..3) removes a polynomial regional trend before inversion. The returned grid uses elevation increasing upward. random_seed makes SimPEG’s eigenvalue-based beta estimate reproducible. Raises InversionBackendUnavailable if SimPEG is missing.

PyHydroGeophysX.inversion.run_joint_inversion(request: JointInversionRequest | Mapping[str, Any], progress: Any | None = None) → JointInversionResult

Validate and execute a registered joint or cooperative inversion.

PyHydroGeophysX.inversion.run_tdem_inversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, n_layers: int = 25, use_irls: bool = True, verbose: bool = True, **kwargs) → TDEMInversionResult

Convenience function to run TDEM inversion.

Parameters:

• times – Time channels (s)

• dobs – Observed data (T)

• uncertainties – Data uncertainties (T)

• source_radius – Loop radius (m)

• n_layers – Number of inversion layers

• use_irls – Use IRLS for sparse inversion

• verbose – Print progress

• **kwargs – Additional parameters for TDEMInversion

Returns:

TDEMInversionResult with inversion results

PyHydroGeophysX.inversion.tdem_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert a TDEM sounding for a layered resistivity model (Occam 1D).

PyHydroGeophysX.inversion.tdem_joint_invert(data: ~typing.Dict[str, ~typing.Any], geom: ~typing.Dict[str, ~typing.Any], inv: ~typing.Dict[str, ~typing.Any], log: ~typing.Callable[[str], None] = <function noop>) → Dict[str, Any]

Invert all available LM/HM gates at one station for one shared 1D model.