seismic.inversion.wavefield_decomp package


seismic.inversion.wavefield_decomp.call_count_decorator module

Decorator to count number of times a function is called.


Decorator to count calls to a function. The number of calls can be queryied from func.counter.


func – Function whose calls to count


func wrapper

seismic.inversion.wavefield_decomp.plot_nd_batch module

Batch plotting a MCMC solution for batch of stations to a single pdf file.

seismic.inversion.wavefield_decomp.plot_nd_batch.convert_Vs_to_k(soln, config)[source]

Transform Vs variable into k variable in MCMC solution. Modifies soln in-place.

  • soln (Customized scipy.optimize.OptimizeResult) – Solution container

  • config (dict) – Solution configuration



seismic.inversion.wavefield_decomp.plot_nd_batch.plot_aux_data(soln, config, log, scale)[source]

Plot auxiliary data such as energy distribution and receiver functions.

  • soln (Customized scipy.optimize.OptimizeResult) – Solution container

  • config (dict) – Solution configuration

  • log (logging.Logger) – Logging instance

  • scale (float) – Overall image scaling factor


Matplotlib figure containing the plotted data

Return type


seismic.inversion.wavefield_decomp.runners module

Batch execution interfaces for wavefield continuation methods and solvers.

seismic.inversion.wavefield_decomp.runners.curate_seismograms(data_all, curation_opts, logger, rotate_to_zrt=True)[source]

Curation function to remove bad data from streams. Note that this function will modify the input dataset during curation.


None, curation operates directly on data_all

seismic.inversion.wavefield_decomp.runners.load_mcmc_solution(h5_file, job_timestamp=None, logger=None)[source]

Load Monte Carlo Markov Chain solution from HDF5 file.

  • h5_file (str or pathlib.Path) – File from which to load solution

  • job_timestamp (str or NoneType) – Timestamp of job whose solution is to be loaded

  • logger (logging.Logger) – Output logging instance


(solution, job configuration), job timestamp

Return type

(solution, dict), str

seismic.inversion.wavefield_decomp.runners.mcmc_solver_wrapper(model, obj_fn, mantle, Vp, rho, flux_window)[source]

Wrapper callable for passing to MCMC solver in scipy style, which unpacks inputs into vector variables for solver.

  • model (numpy.array) – Per-layer model values (the vector being solved for) as flat array of (H, Vs) value pairs ordered by layer.

  • obj_fn – Callable to WfContinuationSuFluxComputer to compute SU flux.

  • mantle (seismic.model_properties.LayerProps) – Mantle properties stored in class LayerProps instance.

  • Vp (numpy.array) – Array of Vp values ordered by layer.

  • rho (numpy.array) – Array of rho values ordered by layer.

  • flux_window ((float, float)) – Pair of floats indicating the time window over which to perform SU flux integration


Integrated SU flux energy at top of mantle

seismic.inversion.wavefield_decomp.runners.run_mcmc(waveform_data, config, logger)[source]

Top level runner function for MCMC solver on SU flux minimization for given settings.

  • waveform_data (Iterable obspy.Stream objects) – Collection of waveform data to process

  • config (dict) – Dict of job settings. See example files for fields and format of settings.

  • logger (logging.Logger or NoneType) – Log message receiver. Optional, pass None for no logging.


Solution object based on scipy.optimize.OptimizeResult

Return type

scipy.optimize.OptimizeResult with additional custom attributes

seismic.inversion.wavefield_decomp.runners.run_station(config_file, waveform_file, network, station, location, logger)[source]

Runner for analysis of single station. For multiple stations, set up config file to run batch job using mpi_job CLI.

The output file is in HDF5 format. The configuration details are added to the output file for traceability.

  • config_file (dict) – Config filename specifying job settings

  • waveform_file (str or pathlib.Path) – Event waveform source file for seismograms, generated using script

  • network (str) – Network code of station to analyse

  • station (str) – Station code to analyse

  • location (str) – Location code of station to analyse. Can be ‘’ (empty string) if not set.

  • logger (logging.Logger) – Output logging instance


Pair containing (solution, configuration) containers. Configuration will have additional traceability information.

Return type

(solution, dict)

seismic.inversion.wavefield_decomp.runners.save_mcmc_solution(soln_configs, input_file, output_file, job_timestamp, job_tracking, logger=None)[source]

Save solution to HDF5 file. In general soln_configs will be an ensemble of per-station solutions.

  • soln_configs (list((solution, dict))) – List of (solution, configuration) pairs to save (one per station).

  • input_file (str) – Name of input file used for the job. Saved to job node for traceability

  • output_file (str or pathlib.Path) – Name of output file to create

  • job_timestamp (str) – Job timestamp that will be used to generate the top level job group

  • job_tracking (dict or dict-like) – Dict containing job identification information for traceability

  • logger (logging.Logger) – [OPTIONAL] Log message destination



seismic.inversion.wavefield_decomp.solvers module

Objective function minimization solvers.

class seismic.inversion.wavefield_decomp.solvers.AdaptiveStepsize(stepper, accept_rate=0.5, interval=50, factor=0.9, ar_tolerance=0.05)[source]

Bases: object

Class to implement adaptive stepsize. Pulled from scipy.optimize.basinhopping and adapted.

class seismic.inversion.wavefield_decomp.solvers.BoundedRandNStepper(bounds, initial_step=None)[source]

Bases: object

Step one dimensional at a time using normal distribution of steps within defined parameter bounds.

property stepsize
class seismic.inversion.wavefield_decomp.solvers.HistogramIncremental(bounds, nbins=20)[source]

Bases: object

Class to incrementally accumulate N-dimensional histogram stats at runtime.

property bins
property dims
property histograms
class seismic.inversion.wavefield_decomp.solvers.SolverGlobalMhMcmc[source]

Bases: object

Drop-in custom solver for scipy.optimize.minimize, based on Metrolpolis-Hastings Monte Carlo Markov Chain random walk with burn-in and adaptive acceptance rate, followed by N-dimensional clustering.

Rather than returning one global solution, the solver returns the N best ranked local solutions. It also returns a probability distribution of each unknown based on Monte Carlo statistics.

seismic.inversion.wavefield_decomp.solvers.optimize_minimize_mhmcmc_cluster(objective, bounds, args=(), x0=None, T=1, N=3, burnin=100000, maxiter=1000000, target_ar=0.4, ar_tolerance=0.05, cluster_eps=0.05, rnd_seed=None, collect_samples=None, logger=None)[source]

Minimize objective function and return up to N local minima solutions.

  • objective (Callable(*args) -> float) – Objective function to minimize. Takes unpacked args as function call arguments and returns a float.

  • bounds (scipy.optimize.Bounds) – Bounds of the parameter space.

  • args (tuple or list) – Any additional fixed parameters needed to completely specify the objective function.

  • x0 (numpy.array with same shape as elements of bounds) – Initial guess. If None, will be selected randomly and uniformly within the parameter bounds.

  • T (float) – The “temperature” parameter for the accept or reject criterion. To sample the domain well, should be in the order of the typical difference in local minima objective valuations.

  • N (int) – Maximum number of minima to return

  • burnin (int) – Number of random steps to discard before starting to accumulate statistics.

  • maxiter (int) – Maximum number of steps to take (including burnin).

  • target_ar (float between 0 and 1) – Target acceptance rate of point samples generated by stepping.

  • ar_tolerance (float) – Tolerance on the acceptance rate before actively adapting the step size.

  • cluster_eps (float) – Point proximity tolerance for DBSCAN clustering, in normalized bounds coordinates.

  • rnd_seed (int) – Random seed to force deterministic behaviour

  • collect_samples (int or NoneType) – If not None and integral type, collect collect_samples at regular intervals and return as part of solution.

  • logger – Logger instance for outputting log messages.


OptimizeResult containing solution(s) and solver data.

Return type

scipy.optimize.OptimizeResult with additional attributes

seismic.inversion.wavefield_decomp.wavefield_continuation_tao module

Class encapsulating algorithm for 1D inversion using wavefield continuation.

Based on reference:

  • Kai Tao, Tianze Liu, Jieyuan Ning, Fenglin Niu, “Estimating sedimentary and crustal structure using wavefield continuation: theory, techniques and applications”, Geophysical Journal International, Volume 197, Issue 1, April, 2014, Pages 443-457,

class seismic.inversion.wavefield_decomp.wavefield_continuation_tao.WfContinuationSuFluxComputer(station_event_dataset, f_s, time_window, cut_window)[source]

Bases: object

Implements computation of the upwards mean S-wave energy flux at the top of the mantle for an ensemble of events for one station.

Class instance can be loaded with a dataset and then evaluated for arbitrary 1D earth models.


  1. Load data from a station event dataset. Copy of the data is buffered by this class in efficient format for energy flux calculation.

  2. Define 1D earth model and mantle half-space material properties (in external code).

  3. Call instance with models and receive energy flux results.


  • station_event_dataset – Iterable container of obspy.Stream objects.

  • f_s – Processing sample rate. Usuually much less than the sampling rate of the input raw seismic traces.

  • time_window – Time window about onset to use for wave continuation processing.

  • cut_window – Shorter time segment within time_window to from which to extract primary arrival waveform and its multiples.

Compute SU energy flux over a grid of H and k values for a particular (single) layer. The layer to compute over is indicated by layer_index, which is a zero-based index into layer_props.

  • mantle_props (LayerProps) – Mantle bulk properties

  • layer_props (list(seismic.model_properties.LayerProps)) – Layer bulk properties as a list

  • layer_index (int) – Index of layer to vary.

  • H_vals (numpy.array) – Set of thickness (H) values to apply to the selected layer.

  • k_vals (numpy.array) – Set of k values to apply to the selected layer.

  • flux_window (pair of numeric values) – Time window in which to compute energy flux

  • ncpus (int) – Number of CPUs to use. Use -1 to use all.


tuple(H grid, k grid, energy values)

Return type

tuple(numpy.array, numpy.array, numpy.array)


Given a stack of layers in layer_props, propagate the surface seismograms down to the base of the bottom layer.


layer_props (list(seismic.model_properties.LayerProps)) – List of LayerProps through which to propagate the surface seismograms.


(Vr, Vz) velocity components time series per event at the bottom of the stack of layers.

Return type



Get array of discrete time values of the time axis


Array of evenly spaced time values

Return type


seismic.inversion.wavefield_decomp.wfd_plot module

Plotting helper functions for Wavefield Decomposition module.

seismic.inversion.wavefield_decomp.wfd_plot.plot_Esu_space(H, k, Esu, title=None, savefile_name=None, show=True, c_range=None, None, decorator=None)[source]

Plot SU energy as function of H-k.

  • H (numpy.array) – Grid of H coordinates corresponding to Esu point values.

  • k (numpy.array) – Grid of k coordinates corresponding to Esu point values.

  • Esu (numpy.array) – Energy values on H-k grid.

  • title (str) – Plot title [OPTIONAL]

  • savefile_name (str) – Output file in which to save plot [OPTIONAL]

  • show (bool) – If True, display the image and block until it is closed.

  • c_range (tuple(float, float)) – Custom range of Esu to contour (min, max values)

  • decorator (Callable) – Callback function to customize plot.



seismic.inversion.wavefield_decomp.wfd_plot.plot_Nd(soln, title='', scale=1.0, vars=None)[source]

Plotting routine for N-dimensional solution in grid format. Diagonal contains histograms of solution samples for each variable, and off-diagonal contains pair-wise scatter plots of sample points plus solution clusters. The distinct solutions coordinates are colour coded and labelled on the histograms.

This function is intended to be a domain-agnostic means of plotting an N-dimensional solution generated by optimize_minimize_mhmcmc_cluster(). It should not have seismology-specific fields or terminology added into it.

  • soln (scipy.optimize.OptimizeResult with additional custom fields) – Solution container returned from optimize_minimize_mhmcmc_cluster()

  • title (str, optional) – Overall plot title, defaults to ‘’ (no title)

  • scale (float, optional) – Scale size of overall plot, defaults to 1.0. Adjust this depending on dimensionality or desired size.

  • vars (list(str)) – List of names of variables, should be same length as dimensionality of solution


Tuple containing the PairGrid, list of axes for the secondary (histogram) axes in the diagonal of the grid, and list of text items representing solution labels.

Return type

tuple(seaborn.PairGrid, list(matplotlib.axes.Axes), list(matplotlib.text.Text))

Module contents