relmt.ls module#

Functions to set up and solve the linear systems

relmt.ls.bootstrap_lsmr(A, b, ev_scale, peq, seq, bootstrap_samples, seed=0, ncpu=1)[source]#

Draw bootstrap samples from A and solve the modified linear system Am = b

Rows of A a drawn with repetion, where: - Rows representing P observations are drawn individually. - Rows representing S observations are drawn as pairs that contain the same Bacb and Bacb relative amplitudes. - Rows representing the reference MT are always included

Parameters:

  • A (ndarray) – (N, M) martix of the linear system Am = b

  • b (ndarray) – (1, N) solution vector of the linear system Am = b

  • ev_scale (float) – (1, N) event-wise condition factors A was devided by

  • peq (int) – Number of equations in A representing P and S wave constraints, respectivley

  • seq (int) – Number of equations in A representing P and S wave constraints, respectivley

  • bootstrap_samples (int) – Number of random realizations of the linear system

  • seed (int) – Random number seed supplied to np.random.default_rng

  • ncpu (int) – Number of parallel process to launch

Returns:

ndarray – (bootstrap_samples * M) bootstrap model vectors

relmt.ls.condition_by_norm(mat, n_homogenous=None)[source]#

Weights to condition the homogenous part of the amplitude matrix

Each weight is the inverse L2-norm of the non-zero values of each equation

Parameters:

  • mat (ndarray) – (N, ev*nmt) amplitude matrix

  • n_homogenous (int | None) – Number of rows that constitute the homogenous part of the matrix. If given, only produce weights for the first ‘n_homogenous’ lines of the matrix.

Returns:

ndarray – (1, ev*nmt) column vector of equation normalization factors

relmt.ls.dircoeff_deviatoric_p(gamma)[source]#

5-element direction coefficeint vector for deviatoric moment tensor

Implementation of Ploudre & Bostock (2019, GJI), Eq. A8

Parameters:

gamma (ndarray) – Direction cosines produced by gamma()

Returns:

ndarray – (5,) directional coefficients

relmt.ls.dircoeff_deviatoric_s(gamma)[source]#

3 5-element direction coefficeint vectors for deviatoric moment tensor

Implementation of Ploudre & Bostock (2019, GJI), Eq. A9

Parameters:

gamma (ndarray) – Direction cosines produced by gamma()

Returns:

tuple[ndarray, ndarray, ndarray] – (5,) directional coefficients

relmt.ls.dircoeff_general_p(gamma)[source]#

6-element direction coefficeint vector for unconstrained moment tensor

Implementation of Ploudre & Bostock (2019, GJI), Eq. A2

Parameters:

gamma (ndarray) – Direction cosines produced by gamma()

Returns:

ndarray – (6,) directional coefficients

relmt.ls.dircoeff_general_s(gamma)[source]#

3 6-element direction coefficeint vectors for unconstrained moment tensor

Implementation of Ploudre & Bostock (2019, GJI), Eq. A3

Parameters:

gamma (ndarray) – Direction cosines produced by gamma()

Returns:

tuple[ndarray, ndarray, ndarray] – (6,) directional coefficients

relmt.ls.distance_ratio(event_a, event_b, station)[source]#

Ratio of the distance of event_a and event_b to station

R = distance(event_a, station) / distance(event_b, station)

Parameters:
Return type:

float | ndarray

relmt.ls.gamma(azimuth, plunge)[source]#

Produce 3-element direction cosine vector

Parameters:

  • azimuth (float) – Degree east of north

  • plunge (float) – Degree down from horizontal

Returns:

ndarray – (3,) direction cosines along north (x), east (y), down (z) axis

relmt.ls.homogenous_equations(p_amplitudes, s_amplitudes, in_events, station_dictionary, event_dict, phase_dictionary, constraint, s_coefficients=None, two_s_equations=True)[source]#

Homogenous part of the linear system Am = b

Some more description of the function

Parameters:

  • p_amplitudes (list[P_Amplitude_Ratio]) – P observations to include in the system

  • s_amplitudes (list[S_Amplitude_Ratios]) – S observations to include in the system

  • in_events (list[int]) – Included events. Columns blocks of A and row blocks of b and m will reference into this vector.

  • station_dictionary (dict[str, Station]) – Lookup table for station coordinates

  • event_dict (dict[int, Event]) – The seismic event catalog

  • phase_dictionary (dict[str, Phase]) – Lookup table for ray take-off angles

  • constraint (str) – Constraint on the moment tensor solution

  • s_coefficients (tuple[int, int] | None) – Indices (0, 1, or 2) of S-wave directional coefficients to use. If None, chose those with the largest norm

  • two_s_equations (bool) – Keep the second S equation with the lower norm of the polarization vector

Returns:

  • Ah (numpy.ndarray) – (p_amplitudes + 2 * s_amplitudes, event_dict * mt_elements) left-hand side of the homogenous part of linear system

  • bh (numpy.ndarray) – (event_dict * mt_elements, 1) zero column vector, right-hand side of the homogenous linear system

relmt.ls.homogenous_equations_sparse(p_amplitudes, s_amplitudes, in_events, station_dictionary, event_dict, phase_dictionary, constraint, s_coefficients=None, two_s_equations=True, ncpu=1)[source]#

Homogenous part of the linear system Am = b

This functions writes the matrices directly in scipy.sparse.coo_matrix format

Parameters:

  • p_amplitudes (list[P_Amplitude_Ratio]) – P observations to include in the system

  • s_amplitudes (list[S_Amplitude_Ratios]) – S observations to include in the system

  • in_events (list[int]) – Included events. Columns blocks of A and row blocks of b and m will reference into this vector.

  • station_dictionary (dict[str, Station]) – Lookup table for station coordinates

  • event_dict (dict[int, Event]) – The seismic event catalog

  • phase_dictionary (dict[str, Phase]) – Lookup table for ray take-off angles

  • constraint (str) – Constraint on the moment tensor solution

  • s_coefficients (tuple[int, int] | None) – Indices (0, 1, or 2) of S-wave directional coefficients to use. If None, chose those with the largest norm

  • ncpu (int) – Number of parallel processes

  • two_s_equations (bool) – Keep the second S equation with the lower norm of the polarization vector

Returns:

  • Ah (scipy.sparse.coo_matrix) – (p_amplitudes + 2 * s_amplitudes, event_dict * mt_elements) left-hand side of the homogenous part of linear system

  • bh (numpy.ndarray) – (event_dict * mt_elements, 1) zero column vector, right-hand side of the homogenous linear system

relmt.ls.mt_elements(constraint)[source]#

Number of elements of the moment tensor

Parameters:

constraint (str) –

  • ‘none’ - invert for full moment tensor

  • ’deviatoric’ - no isotropic component

Returns:

int – 6 for ‘none’, 5 for ‘deviatoric’

Raises:

ValueError: – If unknown constaint

relmt.ls.norm_event_median_amplitude(mat, nmt, n_homogenous=None)[source]#

Factors to normalize amplitude matrix by median value per event

Weight is the inverse median of all non-zero absolute vlues

Parameters:

  • mat (ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]) – (N, ev*nmt) amplitude matrix

  • nmt (int) – Number of free moment tensor paramters

  • n_homogenous (int | None) – Homogenous equations (=number of amplitude observations) in the matrix. If given, only consider the first ‘n_homogenous’ lines in the matrix to estimate median amplitude

Returns:

norms (numpy.ndarray) – (1, ev*nmt) column vector of event norms

Raises:

IndexError: – If number of columns in matrix does not match ‘nmt’

relmt.ls.p_equation(p_amplitude, in_events, station, event_dict, phase_dictionary, nmt, return_sparse=False)[source]#

Equation with P wave constraint on the liniear system

Relative amplitudes are corrected by direction coefficient g and the distance ratio from the event pair to the station.

Parameters:

  • p_amplitude (P_Amplitude_Ratio) – One relative P amplitude observation

  • in_events (list) – Events part of the linear system

  • station (Station) – Station on which the observation has been made

  • event_dict (dict[int, Event]) – The full seismic event catalog

  • phase_dictionary (dict[str, Phase]) – All phase observations

  • nmt (int) – Number of moment tensor elements

  • return_sparse (bool) – Return values data and column indices compatible with scipy.sparse.coo_matrix

Returns:

ndarray | tuple[ndarray, ndarray] – * (1, events * nmt) line of the left hand side of the linear system if * return_sparse = False – Tuple of (2 * nmt,) column indices and (2 * nmt,) values of the * sparse matrix if return_sparse = True

relmt.ls.reference_mt_equations(reference_events, refmt_dict, included_events, constraint)[source]#

Inhomogenous part of the linear system Am = b

Absolute moment tensor constraint on the linear system

Parameters:

  • reference_events (list[int]) – Indices to the reference events

  • refmt_dict (dict[int, MT]) – Lookup table from event index to moment tensor

  • included_events (list[int]) – Event IDs that form the columns of the linear system

  • constraint (str) – Constraint to impose on the moment tensor

Returns:

  • Ai (numpy.ndarray) – (reference_events * mt_elements, number_events * mt_elements) Left-hand side of the inhomogenous part of the linear system.

  • bi (numpy.ndarray) – (reference_events * mt_elements, 1) Elements of the reference moment tensor, right-hand side of the inhomogenous part of the linear system.

relmt.ls.reference_mt_event_norm(ev_norm, ref_mts, nmt)[source]#

Pull event normalization factor for reference moment tensor out of matrix event normalization

Parameters:
Return type:

ndarray

relmt.ls.s_equations(s_amplitude, in_events, station, event_dict, phase_dictionary, nmt, coefficient_indices=None, return_sparse=False, two_s_equations=True)[source]#

Two equations with S-wave constraints on the linear system

Relative amplitudes are corrected by direction coefficient g and the distance ratio from the event pair to the station.

Parameters:

  • s_amplitude (S_Amplitude_Ratios) – One pair of relative S amplitude observations

  • in_events (list) – Events part of the linear system

  • station (Station) – Station on which the observation has been made

  • event_dict (dict[int, Event]) – The full seismic event catalog

  • phase_dictionary (dict[str, Phase]) – All phase observations

  • nmt (int) – Number of moment tensor elements

  • coefficient_indices (tuple[int, int] | None) – Select two indices (0, 1, or 2) of directional coefficients to use. If None, chose those with the largest norm

  • return_sparse (bool) – Return values data and column indices compatible with scipy.sparse.coo_matrix

  • two_s_equations (bool) – Keep the second S equation with the lower norm of the polarization vector

Returns:

ndarray | tuple[ndarray, ndarray] – * (2, events * nmt) lines of the left hand side of the linear system if * return_sparse = False, and two_s_equations = True, or – (1, events * nmt) if return_sparse = False, and * two_s_equations = False – Tuple of (6 * nmt,) column indices and (6 * nmt,) values of the * sparse matrix if return_sparse = True and two_s_equations = False. * The first set of 3*nmt indices and values corresponds to one line, the * second set to the next. – Tuple of (3 * nmt,) column indices and (3 * nmt,) values of the * sparse matrix if return_sparse = True and two_s_equations = False.

relmt.ls.solve_irls_sparse(G, d, tolerance=1e-05, eps=1e-06, efac=1.3)[source]#

Iteratively reweighted least squares for sparse matrices

Solves for m in d = Gm using scipy.sparse.linalg.lsmr() and iteratiley re-scaling d

Parameters:

  • G (spmatrix) – Matrix G in the linear system

  • d (ndarray) – Vector d in the linear system (pair-wise time shifts)

  • tolerance (float) – Average absolute change in model update mean(abs(m - m0)) to stop iteration. Should fall below sampling intervall (Units of m)

  • eps (float) – Truncate smallest residuals to this value to stabilize 1/residual weighting (Units of d)

  • efac (float) – Multiply epsilon by this value when iterative solution stagnates

Returns:

  • m (numpy.ndarray) – Model vector (absolute time shifts per trace)

  • r (numpy.ndarray) – Data residuals (unaccounted pair wise time shifts)

relmt.ls.solve_lsmr(A, b, ev_scale)[source]#

Solve linear system using scipy.sparse.linalg.lsmr()

Parameters:

  • A (ndarray) – Matrix A in the linear system Ax = b of shape (M, N)

  • b (ndarray) – Data vector b (1,M) in the linear system Ax = b

  • ev_scale (ndarray) – Event-wise condition factors (1,M) A was devided by

Returns:

  • x (numpy.ndarray) – The resulting model vector

  • eps (numpy.ndarray) – Residuals Am - b -> 0. Residuals are computed before multiplying m with ev_scale

relmt.ls.unpack_equation_vector(eq_vector, p_lines, ref_lines, nmt, two_s_equations=True)[source]#

Split vector pertaining to equations of the amplitude matrix (e.g. residuals, equation norm) into P-, S-, and reference MT parts

Parameters:

  • eq_vector (ndarray) – Reisdual vector

  • p_lines (int) – Number of equations with P-information

  • ref_lines (int) – Number of equations with reference MT information

  • nmt (int) – Number of elements that parameterize an MT

  • two_s_equations (bool) – We kept the second S equation

Returns:

p_part, s_part, ref_part (numpy.ndarray) – Vector contianing the part pertaining to the P-, S-, and reference equations

relmt.ls.weight_misfit(amplitude, min_misfit, max_misfit, min_weight, phase, two_s_equations=True)[source]#

Weights [min_weight … 1] for each row of the amplitude matrix by misfit

Maps misfit values above max_misfit to min_weight and below min_misfit to unit weight.

weight = 1 - (mis - min_misfit) * (1 - min_weight) / (max_misfit - min_misfit)

Note

A misfit of 0 indicates that the reconstruction and the target waveform are identical. A misfit of 1 that the residual of the waveform reconstruction has the same norm that the waveform itself.

Parameters:

  • amplitude (S_Amplitude_Ratios | P_Amplitude_Ratio) – One amplitude ratio

  • min_misfit (float) – Minimum misfit that receives unit weight

  • max_misfit (float) – Maximum allowed misfit receives zero weight

  • min_weight (float) – Weight of the maximum misfit.

  • phase (str) – If ‘S’, return each weight twice to account for two equations per S wave observation

  • two_s_equations (bool)

Returns:

ndarray – Weight of shape (1,) if phase=P, or (2,1) if phase=S

Raises:

ValueError: – If ‘phase’ is not ‘P’ or ‘S’

relmt.ls.weight_s_amplitude(s_amplitudes, two_s_equations=True)[source]#

Weights for each row of the amplitude matrix by S-wave amplitdue

Weight is the inverse of the larger amplitude, but not more than 1.

weight = max(1, (1 / max(abs(ampl.amp_abc, ampl.amp_acb))))

Parameters:

  • s_amplitudes (S_Amplitude_Ratios) – One pair of S-amplitude ratios

  • two_s_equations (bool) – Keep the second S equation with the lower norm of the polarization vector

Returns:

ndarray – * (2, 1) column vector of weights if two_s_equations is True, * else * (1, 1) column vector of weights