|
5 | 5 | from scipy.stats import median_abs_deviation |
6 | 6 |
|
7 | 7 |
|
8 | | -def hankel_matrix(x, num_delays, pad=False): # fixed delay step of 1 |
9 | | - """Unused throughout the repo |
10 | | -
|
11 | | - :param np.array[float] x: data |
12 | | - :param int num_delays: number of times to 1-step shift data |
13 | | - :param bool pad: if True, return width is len(x), else width is len(x) - num_delays + 1 |
14 | | -
|
15 | | - :return: a Hankel Matrix `m`. For example, if `x = [a, b, c, d, e]` and `num_delays = 3`:\n |
16 | | - With `pad = False`::\n |
17 | | - m = [[a, b, c], |
18 | | - [b, c, d], |
19 | | - [c, d, e]]\n |
20 | | - With `pad = True`::\n |
21 | | - m = [[a, b, c, d, e], |
22 | | - [b, c, d, e, 0], |
23 | | - [c, d, e, 0, 0]] |
24 | | - """ |
25 | | - m = x.copy() |
26 | | - for d in range(1, num_delays): |
27 | | - xi = x[d:] |
28 | | - xi = np.pad(xi, (0, len(x)-len(xi)), 'constant', constant_values=0) |
29 | | - m = np.vstack((m, xi)) |
30 | | - if not pad: |
31 | | - return m[:, 0:-1*num_delays+1] |
32 | | - return m |
33 | | - |
34 | | - |
35 | | -def matrix_inv(X, max_sigma=1e-16): |
36 | | - """Stable (pseudo) matrix inversion using singular value decomposition. Unused throughout the repo. |
37 | | -
|
38 | | - :param np.array X: matrix to invert |
39 | | - :param float max_sigma: smallest singular values to take into account. matrix will be truncated prior to inversion based on this value. |
40 | | -
|
41 | | - :return: (np.array) -- pseudo inverse |
42 | | - """ |
43 | | - U, Sigma, V = np.linalg.svd(X, full_matrices=False) |
44 | | - Sigma_inv = Sigma**-1 |
45 | | - Sigma_inv[np.where(Sigma < max_sigma)[0]] = 0 # helps reduce instabilities |
46 | | - return V.T.dot(np.diag(Sigma_inv)).dot(U.T) |
47 | | - |
48 | | - |
49 | | -def peakdet(x, delta, t=None): |
50 | | - """Find peaks and valleys of 1D array. A point is considered a maximum peak if it has the maximal |
51 | | - value, and was preceded (to the left) by a value lower by delta. Converted from MATLAB script at |
52 | | - http://billauer.co.il/peakdet.html Eli Billauer, 3.4.05 (Explicitly not copyrighted). This function |
53 | | - is released to the public domain; Any use is allowed. |
54 | | -
|
55 | | - :param np.array[float] x: array for which to find peaks and valleys |
56 | | - :param float delta: threshold for finding peaks and valleys. A point is considered a maximum peak |
57 | | - if it has the maximal value, and was preceded (to the left) by a value lower by delta. |
58 | | - :param np.array[float] t: optional domain points where data comes from, to make indices into locations |
59 | | -
|
60 | | - :return: tuple[np.array, np.array] of\n |
61 | | - - **maxtab** -- indices or locations (column 1) and values (column 2) of maxima |
62 | | - - **mintab** -- indices or locations (column 1) and values (column 2) of minima |
63 | | - """ |
64 | | - maxtab = [] |
65 | | - mintab = [] |
66 | | - if t is None: |
67 | | - t = np.arange(len(x)) |
68 | | - elif len(x) != len(t): |
69 | | - raise ValueError('Input vectors x and t must have same length') |
70 | | - if not (np.isscalar(delta) and delta > 0): |
71 | | - raise ValueError('Input argument delta must be a positive scalar') |
72 | | - |
73 | | - mn, mx = np.inf, -1*np.inf |
74 | | - mnpos, mxpos = np.nan, np.nan |
75 | | - lookformax = True |
76 | | - for i in np.arange(len(x)): |
77 | | - this = x[i] |
78 | | - if this > mx: |
79 | | - mx = this |
80 | | - mxpos = t[i] |
81 | | - if this < mn: |
82 | | - mn = this |
83 | | - mnpos = t[i] |
84 | | - if lookformax: |
85 | | - if this < mx-delta: |
86 | | - maxtab.append((mxpos, mx)) |
87 | | - mn = this |
88 | | - mnpos = t[i] |
89 | | - lookformax = False # now searching for a min |
90 | | - else: |
91 | | - if this > mn+delta: |
92 | | - mintab.append((mnpos, mn)) |
93 | | - mx = this |
94 | | - mxpos = t[i] |
95 | | - lookformax = True # now searching for a max |
96 | | - |
97 | | - return np.array(maxtab), np.array(mintab) |
98 | | - |
99 | | - |
100 | 8 | def huber(x, M): |
101 | 9 | """Huber loss function, for outlier-robust applications, `see here <https://www.cvxpy.org/api_reference/cvxpy.atoms.elementwise.html#huber>`_""" |
102 | 10 | absx = np.abs(x) |
|
0 commit comments