def svd_moments(u, s, v, stachans, event_list, n_svs=<NUM_LIT:2>):
<EOL>K_width = max([max(ev_list) for ev_list in event_list]) + <NUM_LIT:1><EOL>if len(stachans) == <NUM_LIT:1>:<EOL><INDENT>print('<STR_LIT>'<EOL>'<STR_LIT>')<EOL>return u[<NUM_LIT:0>][:, <NUM_LIT:0>], event_list[<NUM_LIT:0>]<EOL><DEDENT>for i, stachan in enumerate(stachans):<EOL><INDENT>k = []  <EOL>U_working = copy.deepcopy(v[i].T)<EOL>V_working = copy.deepcopy(u[i])<EOL>s_working = copy.deepcopy(s[i].T)<EOL>ev_list = event_list[i]<EOL>if len(ev_list) > len(U_working):<EOL><INDENT>print('<STR_LIT>' + str(U_working.shape))<EOL>print('<STR_LIT>' % str(len(ev_list)))<EOL>f_dump = open('<STR_LIT>', '<STR_LIT:wb>')<EOL>pickle.dump(U_working, f_dump)<EOL>f_dump.close()<EOL>raise IOError('<STR_LIT>')<EOL><DEDENT>s_working[n_svs:len(s_working)] = <NUM_LIT:0><EOL>s_working = np.diag(s_working)<EOL>U_working = np.matrix(U_working)<EOL>V_working = np.matrix(V_working)<EOL>s_working = np.matrix(s_working)<EOL>SVD_weights = U_working[:, <NUM_LIT:0>]<EOL>if np.all(SVD_weights < <NUM_LIT:0>):<EOL><INDENT>warnings.warn('<STR_LIT>')<EOL>SVD_weights = np.abs(SVD_weights)<EOL><DEDENT>SVD_weights = np.array(SVD_weights).reshape(-<NUM_LIT:1>).tolist()<EOL>random_SVD_weights = np.copy(SVD_weights)<EOL>random_SVD_weights = random_SVD_weights.tolist()<EOL>random_SVD_weights = [(random_SVD_weights[_i], _i)<EOL>for _i in range(len(random_SVD_weights))]<EOL>random.shuffle(random_SVD_weights)<EOL>random_SVD_weights.append(random_SVD_weights[<NUM_LIT:0>])<EOL>pairs = []<EOL>for pair in _pairwise(random_SVD_weights):<EOL><INDENT>pairs.append(pair)<EOL><DEDENT>for pairsIndex in range(len(pairs)):<EOL><INDENT>_weights = list(zip(*list(pairs[pairsIndex])))[<NUM_LIT:0>]<EOL>_indeces = list(zip(*list(pairs[pairsIndex])))[<NUM_LIT:1>]<EOL>min_weight = min(np.abs(_weights))<EOL>max_weight = max(np.abs(_weights))<EOL>min_index = _indeces[np.argmin(np.abs(_weights))]<EOL>max_index = _indeces[np.argmax(np.abs(_weights))]<EOL>row = []<EOL>for j in range(len(SVD_weights)):<EOL><INDENT>if j == max_index:<EOL><INDENT>result = -<NUM_LIT:1><EOL><DEDENT>elif j == min_index:<EOL><INDENT>normalised = max_weight / min_weight<EOL>result = float(normalised)<EOL><DEDENT>else:<EOL><INDENT>result = <NUM_LIT:0><EOL><DEDENT>row.append(result)<EOL><DEDENT>k.append(row)<EOL><DEDENT>k_filled = np.zeros([len(k), K_width])<EOL>for j in range(len(k)):<EOL><INDENT>for l, ev in enumerate(ev_list):<EOL><INDENT>k_filled[j, ev] = k[j][l]<EOL><DEDENT><DEDENT>if '<STR_LIT>' not in locals():<EOL><INDENT>K = k_filled<EOL><DEDENT>else:<EOL><INDENT>K = np.concatenate([K, k_filled])<EOL><DEDENT><DEDENT>K_nonempty = []<EOL>events_out = []<EOL>for i in range(<NUM_LIT:0>, K_width):<EOL><INDENT>if not np.all(K[:, i] == <NUM_LIT:0>):<EOL><INDENT>K_nonempty.append(K[:, i])<EOL>events_out.append(i)<EOL><DEDENT><DEDENT>K = np.array(K_nonempty).T<EOL>K = K.tolist()<EOL>K_width = len(K[<NUM_LIT:0>])<EOL>K.append(np.ones(K_width) * (<NUM_LIT:1.> / K_width))<EOL>print("<STR_LIT>")<EOL>del row<EOL>print('<STR_LIT:\n>'.join(['<STR_LIT>'.join([str(round(float(item), <NUM_LIT:3>)).ljust(<NUM_LIT:6>)<EOL>for item in row]) for row in K]))<EOL>Krounded = np.around(K, decimals=<NUM_LIT:4>)<EOL>W = np.matrix(np.identity(len(K)))<EOL>W[-<NUM_LIT:1>, -<NUM_LIT:1>] = len(K) - <NUM_LIT:1><EOL>K = np.matrix(K)<EOL>Kinv = np.array(np.linalg.inv(K.T * W * K) * K.T * W)<EOL>M = Kinv[:, -<NUM_LIT:1>]<EOL>return M, events_out<EOL>
Calculate relative moments/amplitudes using singular-value decomposition.

Convert basis vectors calculated by singular value \
decomposition (see the SVD functions in clustering) into relative \
moments.

For more information see the paper by \
`Rubinstein & Ellsworth (2010).
<http://www.bssaonline.org/content/100/5A/1952.short>`_

:type u: list
:param u:
    List of the :class:`numpy.ndarray` input basis vectors from the SVD,
    one array for each channel used.
:type s: list
:param s:
    List of the :class:`numpy.ndarray` of singular values, one array for
    each channel.
:type v: list
:param v:
    List of :class:`numpy.ndarray` of output basis vectors from SVD, one
    array per channel.
:type stachans: list
:param stachans: List of station.channel input
:type event_list: list
:param event_list: List of events for which you have data, such that \
    event_list[i] corresponds to stachans[i], U[i] etc. and \
    event_list[i][j] corresponds to event j in U[i].  These are a series \
    of indexes that map the basis vectors to their relative events and \
    channels - if you have every channel for every event generating these \
    is trivial (see example).
:type n_svs: int
:param n_svs: Number of singular values to use, defaults to 4.

:returns: M, array of relative moments
:rtype: :class:`numpy.ndarray`
:returns: events_out, list of events that relate to M (in order), \
    does not include the magnitude information in the events, see note.
:rtype: :class:`obspy.core.event.event.Event`

.. note:: M is an array of relative moments (or amplitudes), these cannot
    be directly compared to true moments without calibration.

.. note:: When comparing this method with the method used for creation
    of subspace detectors (Harris 2006) it is important to note that the
    input `design set` matrix in Harris contains waveforms as columns,
    whereas in Rubinstein & Ellsworth it contains waveforms as rows
    (i.e. the transpose of the Harris data matrix). The U and V matrices
    are therefore swapped between the two approaches. This is accounted
    for in EQcorrscan but may lead to confusion when reviewing the code.
    Here we use the Harris approach.

.. rubric:: Example

>>> from eqcorrscan.utils.mag_calc import svd_moments
>>> from obspy import read
>>> import glob
>>> import os
>>> from eqcorrscan.utils.clustering import svd
>>> import numpy as np
>>> # Do the set-up
>>> testing_path = 'eqcorrscan/tests/test_data/similar_events'
>>> stream_files = glob.glob(os.path.join(testing_path, '*'))
>>> stream_list = [read(stream_file) for stream_file in stream_files]
>>> event_list = []
>>> remove_list = [('WHAT2', 'SH1'), ('WV04', 'SHZ'), ('GCSZ', 'EHZ')]
>>> for i, stream in enumerate(stream_list):
...     st_list = []
...     for tr in stream:
...         if (tr.stats.station, tr.stats.channel) not in remove_list:
...             stream.remove(tr)
...             continue
...         tr.detrend('simple')
...         tr.filter('bandpass', freqmin=5.0, freqmax=15.0)
...         tr.trim(tr.stats.starttime + 40, tr.stats.endtime - 45)
...         st_list.append(i)
...     event_list.append(st_list) # doctest: +SKIP
>>> event_list = np.asarray(event_list).T.tolist()
>>> SVec, SVal, U, stachans = svd(stream_list=stream_list) # doctest: +SKIP
['GCSZ.EHZ', 'WV04.SHZ', 'WHAT2.SH1']
>>> M, events_out = svd_moments(u=U, s=SVal, v=SVec, stachans=stachans,
...                             event_list=event_list) # doctest: +SKIP