Inversion¶

Modules for waveform inversion

pytomo.inversion.*

fwi¶

class pytomo.inversion.fwi.FWI(model_ref, model_params, dataset_dict, windows_dict, n_phases, mode, phase_ref=None, buffer=10.0)¶

Implement full-waveform inversion.

Examples

>>> from dsmpy.seismicmodel import SeismicModel
>>> from dsmpy.modelparameters import ModelParameters, ParameterType
>>> from dsmpy.windowmaker import WindowMaker
>>> from dsmpy.dataset import Dataset, read_sac_meta
>>> model_params = ModelParameters(
...    types=[ParameterType.VSH],
...    radii=[3480. + 20 * i for i in range(21)],
...    mesh_type='boxcar')
>>> model_ref = SeismicModel.prem().boxcar_mesh(model_params)
>>> sac_files = [] # list of SAC files to use in inversion
>>> sac_meta, traces = read_sac_meta(sac_files)
>>> windows = WindowMaker.windows_from_obspy_traces(
...     traces, 'prem', ['ScS'], [Component.T],
...     t_before=10, t_after=30)
>>> buffer = 10
>>> WindowMaker.set_limit(
...     windows, t_before + buffer, t_after + buffer)
>>> freqs, freqs2 = [0.01, 0.01], [0.05, 0.08]
>>> dataset_dict = dict()
>>> for freq, freq2 in zip(freqs, freqs2):
...     ds = Dataset.dataset_from_sac_process(
...         sac_files, windows, freq, freq2)
...     dataset_dict[freqency_hash(freq, freq2)] = ds
...     windows_dict[freqency_hash(freq, freq2)] = windows
>>> fwi = FWI(
...    model_ref, model_params, dataset_dict,
...    windows_dict, n_phases=1, mode=2,
...    phase_ref='S', buffer=buffer)
>>> model_1 = fwi.step(
...     model_ref, freq=0.01, freq2=0.05, n_pca_components=[2, 4])
>>> model_2 = fwi.step(
...     model_1, 0.01, 0.08, n_pca_components=[4, 6, 8])

set_ignore_types(parameter_types, ignore=True)¶: Do not update the parameter types. Switch these parameters on again by calling it with ignore=True

set_selection(var: float, corr: float, ratio: float)¶: Set the variance, correlation, and ratio cutoffs used when deciding whether or not to keep a particular record.

step(model, freq, freq2, n_pca_components=[4], alphas=[1.0])¶

Advance one step in the FWI iteration.

Parameters

model (SeismicModel) – the current initial model
freq (float) – minimum frequency for the bandpass filter (in Hz)
freq2 (float) – maximum frequency
n_pca_components (list of int) – number of PCA components to test (default is [4,])

Returns

the updated model

Return type

SeismicModel

pytomo.inversion.fwi.freqency_hash(freq: float, freq2: float) → str¶: Return a hash to use as a key.

pytomo.inversion.fwi.frequencies_from_hash(freq_hash: str) -> (<class 'float'>, <class 'float'>)¶: Return freq_min, freq_max from a frequency hash obtained using frequency_hash().

umcutils¶

Utilities for the UniformMonteCarlo algorithm.

class pytomo.inversion.umcutils.UniformMonteCarlo(model, model_params, range_dict, mesh_type='boxcar', seed=None)¶

Implements the Uniform Monte Carlo method.

Parameters

model (SeismicModel) – reference seismic model.
model_params (ModelParameters) –
range_dict (dict) – entries are of type ParameterType:ndarray.
mesh_type (str) – ‘triangle’ or ‘boxcar’ (default is ‘triangle’).
seed (int) – seed for the random number generator (default is None).

sample_models(ns)¶

Sample ns seismic models using a uniform distribution.

Parameters

ns (int) – number of models.

Returns

list of sampled models. list of ndarray: model perturbations for each sampled

model. Entries are of shape [n_grid_params * ntype,].

Return type

list of SeismicModel

pytomo.inversion.umcutils.correlation(obs, syn)¶

Return the zero-lag cross-correlation misfit.

Parameters

obs (ndarray) – observed waveforms.
syn (ndarray) – synthetics.

Returns

cross-correlation misfit

Return type

float

pytomo.inversion.umcutils.get_best_models(misfit_dict, n_best, key='variance')¶

Get the n_best best models. Should be used together with an InversionResult object.

Parameters

misfit_dict (dict) – a dict with misfit values for each model, as given by InversionResult.misfit_dict.
n_best (int) – number of best models to return.
key (str) – misfit type. (‘variance’, ‘corr’, ‘rolling_variance’), (default is ‘variance’).

Returns

list of the indices of the n_best models in: InversionResult.models

Return type

list of int

pytomo.inversion.umcutils.process_outputs(outputs, dataset, models, windows, freq, freq2, filter_type, **kwargs)¶

Process the output of compute_models_parallel().

Parameters

outputs (list of list of PyDSMOutput) – (n_models, n_events)
dataset (Dataset) – dataset with observed data. Same as the one used for input to compute_models_parallel()
models (list of SeismicModel) – seismic models
windows (list of Window) – time windows. See windows_from_dataset().
freq (float) – frequency for the filter
freq2 (float) – maximum frequency for the filter (only if filter_type==’bandpass’)
filter_type (str) – filter type.
kwargs (**dict) – kwargs for the misfit functions.

Returns

values are ndarray((n_models, n_windows)): containing misfit values (corr, variance).

Return type

dict

pytomo.inversion.umcutils.rolling_variance(obs, syn, size, stride)¶

Return the variances of the obs-syn residual vectors summed over a rolling window. The obs-syn residual is scaled to max(abs(obs)) in each rolling window in order to give the same importance to signals of different amplitudes.

Parameters

obs (ndarray) – observed waveforms.
syn (ndarray) – synthetics.
size (int) – length of the rolling window.
stride (int) – stride length for the rolling window.

Returns

rolling_variance(obs - syn).

Return type

float

pytomo.inversion.umcutils.variance(obs, syn)¶

Return the variance of the obs-syn residual vector.

Parameters

obs (ndarray) – observed waveforms.
syn (ndarra) – synthetics.

Returns

variance(obs - syn)

Return type

float

na¶

Module for global optimization using the Neighbourhood algorithm.

class pytomo.inversion.na.InputFile(input_file)¶

Input file for NeighbourAlgorithm (NA) inversion.

Parameters: input_file (str) – path of NA input file

read()¶

Read from the input file into a dict.

Returns: input file parameters
Return type: dict

class pytomo.inversion.na.NeighbouhoodAlgorithm(dataset, model_ref, model_params, range_dict, tlen, nspc, sampling_hz, mode, n_mod, n_s, n_r, phases, components, t_before, t_after, filter_type, freq, freq2, distance_min, distance_max, convergence_threshold, stf_catalog, result_path, misfit_type, misfit_kwargs, seed, verbose, comm)¶

Implements the Neighbourhood Algorithm

Parameters

dataset (Dataset) – dataset.
model_ref (SeismicModel) – reference seismic model.
model_params (ModelParameters) – model parameters.
range_dict (dict) – range of sampled perturbations. Entries are of type ParameterType:ndarray of shape (n_nodes, 2).
tlen (float) – duration of the synthetics (in seconds) (better to be 2**n/10).
nspc (int) – number of frequency points in the synthetics (better to be 2**n).
sampling_hz (int) – sampling frequency of the synthetics.
mode (int) – computation mode. 0: both, 1: P-SV, 2: SH.
n_mod (int) – maximum number of models sampled.
n_s (int) – number of models at each step of the NA (must have n_mod % n_s = 0)
n_r (int) – number of best-fit models retained at each step of the NA (must have n_mod % n_s = 0)
phases (list of str) – list of seismic phases.
components (list of Component) – seismic components.
t_before (float) – time (in seconds) before each phase arrival time when creating time windows.
t_after (float) – time (in seconds) after each phase arrival time when creating time windows.
filter_type (str) – ‘bandpass’ or ‘lowpass’.
freq (float) – minimum frequency of the filter (in Hz)
freq2 (float) – maximum frequency of the filter (in Hz). Used only for bandpass filter.
distance_min (float) – minimum epicentral distance (in degree).
distance_max (float) – maximum epicentral distance (in degree).
convergence_threshold (float) – convergence threshold
stf_catalog (str) – path to a source time function catalog.
result_path (str) – path to the output folder.
misfit_type (str) – misfit used to select the best models. (‘variance’, ‘corr’, ‘rolling_variance’), (default is ‘variance’)
misfit_kwargs (dict) – kwargs for the misfit function. Used for rolling_variance. See umcutils.rolling_variance.
seed (int) – seed for the random generator.
verbose (int) – verbosity level.
comm (MPI_COMM_WORLD) – MPI communicator.

compute(comm, log=None)¶

Run the NA inversion.

Parameters

comm (COMM_WORLD) – MPI Communicator.
log (str) – path to a log file (default is None).

Returns

inversion result object

Return type

InversionResult

classmethod from_file(input_file, model_ref, model_params, range_dict, dataset, comm)¶

Build a NeighbourhoodAlgorithm object from an input file and key inputs.

Parameters

input_file (str) – path to an input file.
model_ref (SeismicModel) – reference seismic model.
model_params (ModelParameters) – model parameters.
range_dict (dict) – range of sampled perturbations. Entries are of type ParameterType:ndarray of shape (n_nodes, 2).
dataset (Dataset) – dataset
comm (MPI_COMM_WOLRD) – MPI communicator

Returns

NeighbourhoodAlgorithm object

Return type

NeighbourhoodAlgorithm

classmethod from_file_with_default(input_file_path, dataset, comm)¶

Build a NeighbourhoodAlgorithm from an input file and default values.

Parameters

input_file_path (str) – path of the input file
dataset (Dataset) – dataset.
comm (COMM_WORLD) – MPI communicator

Returns

NeighbourhoodAlgorithm object

Return type

NeighbourhoodAlgorithm

get_meta()¶

Return the meta parameters of the NeighbourhoodAlgorithm object.

Returns:
dict: dict with meta parameters

inversionresult¶

class pytomo.inversion.inversionresult.FWIResult(windows_dict)¶

Pack the inversion results for full-waveform inversion and the metadata needed to reproduce the results.

add(models, models_meta, misfits, freq_hash)¶

Add results of an iteration step.

Parameters

models (list of SeismicModel) – inverted models for different metaparameters.
models_meta (list of dict) – entries should correspond to the models list.
misfits (dict) – Keys are misfit names (corr, variance). Values are ndarray of shape (n_models, n_windows) with misfit values.
freq_hash (str) – frequency hash

static load(path)¶

Read file into self using pickle.load().

Parameters: path (str) – name of the file that contains self.

save(path)¶

Save self using pickle.dump().

Parameters: path (str) – name of the output file.

class pytomo.inversion.inversionresult.InversionResult(dataset, windows, meta=None)¶

Pack the inversion results, models and dataset needed to reproduce the results.

Parameters

dataset (Dataset) – dataset
windows (list of Window) – time windows. The order is the same as for axis 1 of misfit_dict.
meta (dict) – dict with meta parameters

add_result(models, misfit_dict, perturbations)¶

Add misfit_dict for new models to current inversion result. The dataset and windows must be the same.

Parameters

models (list of SeismicModel) – seismic models. The order
the same as for axis 0 of misfit_dict. (is) –
misfit_dict (dict) – values are ndarray((n_models, n_windows)) containing misfit values (corr, variance)
perturbations (ndarray) – perturbations

get_model_perturbations(smooth=True, n_s=None, in_percent=False)¶

Get perturbations from model_ref.

Parameters

smooth (bool) – smooth over n_s (default is True).
n_s (int) – number of models computed at each iteration (default is None).
in_percent (bool) – returns perturbations as percent (default is False).

Returns

array of perturbations.: If smooth, shape is (n_iteration,), else shape is (n_models,)

Return type

ndarray

get_variances(key='variance', smooth=True, n_s=None)¶

Get misfit array.

Parameters

key (object) – misfit type (‘variance’, ‘corr’, ‘rolling_variance’) (default is ‘variance’).
smooth (bool) – smooth over n_s (default is True).
n_s (int) – number of models computed at each iteration.

Returns

misfit values

Return type

ndarray

static load(path)¶

Read file into self using pickle.load().

Parameters: path (str) – name of the file that contains self.

plot_models(types, n_best=- 1, n_mod=- 1, ax=None, key='variance', **kwargs)¶

Plot models colored by misfit value.

Parameters

types (list of ParameterType) – parameter types (e.g., ParameterType.RHO).
n_best (int) – number of best models.
n_mod (int) – number of models.
ax (Axes) – matplotlib Axes object.
key (str) – ‘variance’, ‘corr’, or ‘rolling_variance’ (default is ‘variance’).

Returns

matplotlib figure object. Axes: matplotlib Axes object.

Return type

figures

save(path)¶

Save self using pickle.dump().

Parameters: path (str) – name of the output file.

save_convergence_curve(path, scale_arr, free_indices, smooth=True)¶

Save the convergence curve to path.

Parameters

path (str) – path of the output figure.
scale_arr (ndarray) –
free_indices (list of int) –
smooth (bool) – average the variance points over each iteration of the NA algorithm (n_s points) (default is True).

save_variance_curve(path, smooth=True)¶

Save the variance curve to path.

Parameters

path (str) – path to the output figure.
smooth (bool) – average the variance points over each iteration of the NA algorithm (n_s points) (default is True).

voronoi¶

pytomo.inversion.voronoi.compute_distances_to_points(points, p_arr, ips=slice(None, None, None))¶

Parameters

points (ndarray) –
p_arr (ndarray) – coordinates of anchor point
ips (list of int) – indices of target points (default is slice(None))

Returns

square of distances

Return type

ndarray

pytomo.inversion.voronoi.find_point_of_region(vor, ireg)¶

Return the index of point within region ireg.

Parameters

vor (Voronoi) – Voronoi object.
ireg (int) – index of region.

Returns

index of the Voronoi point for region ireg.

Return type

int

pytomo.inversion.voronoi.implicit_find_bound_for_dim(points, anchor_point, current_point, idim, n_nearest=10, min_bound=None, max_bound=None, step_size=0.01, n_step_max=1000, log=None)¶

Find the lower and upper boundaries of region of current_point along dimension idim, without explicitly computing the voronoi diagram, by performing a linear grid search.

Parameters

points (ndarray) – voronoi points
anchor_point (ndarray) – center of the voronoi cell of current_point.
current_point (ndarray) – point from which to compute distances to voronoi cell boundaries.
idim (int) – index of the dimension along which to compute the distance to boundaries.
n_nearest (int) – number of voronoi points used to compute the distance to boundaries.
min_bound (float) – minimum acceptable distance to boundaries. Used to avoid distances larger than largest perturbations defined in range_dict.
max_bound (float) – maximum acceptable distance.
step_size (float) – step size for the grid search (default is 0.01).
n_step_max (int) – maximum number of steps in the grid search (default is 1000).
log (file) – log file (default is None).

Returns

[lower, upper] bounds

Return type

ndarray

pytomo.inversion.voronoi.plot_voronoi_2d(points, misfits, xlim=None, ylim=None, ax=None, **kwargs)¶

Plot a voronoi diagram (only for 2-D points).

Parameters

points (ndarray) – voronoi points. (npoint, 2).
misfits (ndarray) – array of misfits for each voronoi point. (npoint,).
xlim (list of float) – x-axis range (default is None).
ylim (list of float) – y-axis range (default is None).
ax (Axes) – matplotlib Axes object (default is None).

Returns

matplotlib figure object ax (Axes): matplotlib ax object scalar_map (): color map

Return type

fig (figure)

modelutils¶

pytomo.inversion.modelutils.std_boxcar_mesh(n_upper_mantle=20, n_mtz=10, n_lower_mantle=12, n_dpp=6, types=[<ParameterType.VSH: 4>], verbose=0)¶: Boxcar mesh using ak135 as reference model for the structure of the upper mantle and transition zone down to 1000 km depth.

Inversion¶

fwi¶

umcutils¶

na¶

inversionresult¶

voronoi¶

modelutils¶

Table of Contents

Previous topic

Next topic

This Page