Addressing #167 #168
Addressing #167 #168
Conversation
…mpler rmse function. Amended how constants of integration are found to be robust. Optimizer loss function now uses x_hat directly instead of trying to integrate the derivative, since all algorithms return an x_hat estimate associated with dxdt_hat anyway. SpectralDiff now does an IFFT to get x_hat instead of trapezoidal integration. Improved a few docstrings. Got tests to pass with robust change. Nixed a few tests of the optimizer, because they're redundant with each other and the notebooks
| :param bool pad_to_zero_dxdt: if True, extend the data with extra regions that smoothly force the derivative to | ||
| zero before taking FFT. | ||
| :return: tuple[np.array, np.array] of\n |
There was a problem hiding this comment.
I decided I like a slightly shorter form for these.
| filt = np.ones(k.shape); filt[discrete_cutoff:N-discrete_cutoff] = 0 | ||
|
|
||
| # Smoothed signal | ||
| X = np.fft.fft(x) |
There was a problem hiding this comment.
I think this code is a bit clearer and more like what people canonically do now. Integrating to get x_hat was always strange in this spectral case.
|
|
||
| if num_iterations > 1: # We've lost a constant of integration in the above | ||
| x_hat += utility.estimate_integration_constant(x, x_hat) # uses least squares | ||
| x_hat += utility.estimate_integration_constant(x, x_hat) |
There was a problem hiding this comment.
This now emphatically does not use least squares, because that's actually a terrible way to try to estimate the constant of integration if there are outliers. It now uses an outlier-robust fancy-shmansy thing under the hood, which justifies having a whole separate function to this a little more too.
| :param float or array[float] _t: This function supports variable step size. This parameter is either the constant | ||
| step size if given as a single float, or sample locations if given as an array of same length as the state histories. | ||
| :param np.array A: state transition matrix, in discrete time if constant dt, in continuous time if variable dt | ||
| :param list[np.array] xhat_pre: a priori estimates of xhat from a kalman_filter forward pass |
There was a problem hiding this comment.
I turned these into 3D arrays rather than lists of 2D arrays a while ago and hadn't updated the docstring.
| objective = 0.5*cvxpy.sum_squares(proc_resids) if huberM == float('inf') \ | ||
| else np.sqrt(2)*cvxpy.sum(cvxpy.abs(proc_resids)) if huberM < 1e-3 \ | ||
| else huber_const(huberM)*cvxpy.sum(cvxpy.huber(proc_resids, huberM)) # 1/2 l2^2, l1, or Huber | ||
| objective = 0.5*cvxpy.sum_squares(proc_resids) if proc_huberM == float('inf') \ |
There was a problem hiding this comment.
I have to try out whether splitting up proc_huberM and meas_huberM etc gives better results, but I'm reasonably sure we need Nelder-Mead to do a 4D search here, which means it makes a lot of queries. Thankfully, each query is not fast, thanks to vectorized CVXPY code and CLARABEL taking full advantage of sparse matrices.
| {'q': (1e-10, 1e10), | ||
| 'r': (1e-10, 1e10)}), | ||
| # robustdiff: ({'order': {1, 2, 3}, # warning: order 1 hacks the loss function when tvgamma is used, tends to win but is usually suboptimal choice in terms of true RMSE | ||
| # 'log_q': [1., 4, 8, 12], # decimal after first entry ensure this is treated as float type |
There was a problem hiding this comment.
I'm in the midst of experimenting here. Got derailed by #167
| rms_rec_x, rms_x, rms_dxdt = evaluate.rmse(x, dt, x_hat, dxdt_hat, dxdt_truth=None, padding=padding) | ||
| cost = rms_rec_x + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding) | ||
| else: # then minimize sqrt{2*Mean(Huber((x_hat- x)/sigma))}*sigma + gamma*TV(dxdt_hat) | ||
| cost = evaluate.robust_rme(x, x_hat, padding=padding) + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding) |
There was a problem hiding this comment.
This case is now a bit simplified in code here, although it's doing a fancier thing. The evaluation code improvements paying dividends.
| [(-25, -25), (0, -1), (0, 0), (1, 1)], | ||
| [(-25, -25), (1, 1), (0, 0), (1, 1)], | ||
| [(-25, -25), (3, 3), (0, 0), (3, 3)]], | ||
| iterated_second_order: [[(-9, -10), (-25, -25), (0, -1), (0, 0)], |
There was a problem hiding this comment.
Using the finding the integration constant robustly rather than with true least squares hurts performance on these fragile examples a tiny bit. But it's worth it for not living dangerously in a land where the constant can be easily corrupted in a much more significant way by a single outlier.
| assert params1['num_iterations'] == 5 | ||
| assert params2['num_iterations'] == 1 | ||
|
|
||
| def test_mediandiff(): |
There was a problem hiding this comment.
It's silly to run the optimizer this many times on this many different functions. Nixed a few.
| x0 = utility.estimate_integration_constant(x, x_hat, M=float('inf')) | ||
| assert 0.95 < x0 < 1.05 # The result should be close to 1.0, but not exactly due to noise | ||
|
|
||
| x[100] = 100 # outlier case |
There was a problem hiding this comment.
Added and played with an outlier case to get an intuition that my function was working as expected.
| [5.582, -0.31529832], | ||
| [7.135, -0.58327787], | ||
| [8.603, -1.71278265]]) | ||
| assert np.allclose(maxtab, [[0.447, 1.58575613], # these numbers validated by eye with --plot |
There was a problem hiding this comment.
Due to run order change, the random seed behaves a tiny bit differently here.
| better in the presence of outliers""" | ||
| u = np.sin(np.arange(100)*0.1) | ||
| v = u + np.random.randn(100) | ||
| assert np.allclose(evaluate.rmse(u, v), evaluate.robust_rme(u, v, M=6)) |
There was a problem hiding this comment.
Finally added an evaluation code test. It's mildly tricky to see how the robust_rme can equal rmse, but with big enough M, they now produce the same answer. And in the presence of outliers, the robust version is thrown off far less.
| if show_error: | ||
| _, _, rms_dxdt = rmse(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth) | ||
| R_sqr = error_correlation(dxdt_hat, dxdt_truth) | ||
| rms_dxdt = rmse(dxdt_truth, dxdt_hat) |
There was a problem hiding this comment.
rmse now returns only one thing
| _, _, rms_dxdt = rmse(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth) | ||
| R_sqr = error_correlation(dxdt_hat, dxdt_truth) | ||
| rms_dxdt = rmse(dxdt_truth, dxdt_hat) | ||
| R_sqr = error_correlation(dxdt_truth, dxdt_hat) |
There was a problem hiding this comment.
Changed the order of inputs to all functions to have known stuff first and hatted, computed stuff second.
| def rmse(x, dt, x_hat, dxdt_hat, x_truth=None, dxdt_truth=None, padding=0): | ||
| """Evaluate x_hat based on RMSE, calculating different ones depending on whether :code:`dxdt_truth` | ||
| and :code:`x_truth` are known. | ||
| def robust_rme(x, x_hat, padding=0, M=6): |
There was a problem hiding this comment.
I decided to make this a separate function, although it reduces to be the same as rmse when M is infinity, because the two are used for different purposes in our optimization and plotting.
| x0 = utility.estimate_integration_constant(x, rec_x) | ||
| rec_x = rec_x + x0 | ||
| rms_rec_x = np.linalg.norm(rec_x[s] - x[s]) / root | ||
| def rmse(dxdt_truth, dxdt_hat, padding=0): |
There was a problem hiding this comment.
Many fewer parameters
| if padding == 'auto': padding = max(1, int(0.025*len(dxdt_hat))) | ||
| s = slice(padding, len(dxdt_hat)-padding) # slice out data we want to measure | ||
| errors = (dxdt_hat[s] - dxdt_truth[s]) | ||
| r = stats.linregress(dxdt_truth[s] - np.mean(dxdt_truth[s]), errors) |
There was a problem hiding this comment.
Subtracting off the mean here shouldn't be necessary, doesn't change correlation coefficient.
| :return: **integration constant** (float) -- initial condition that best aligns x_hat with x | ||
| """ | ||
| return minimize(lambda x0, x, xhat: np.linalg.norm(x - (x_hat+x0)), # fn to minimize in 1st argument |
There was a problem hiding this comment.
We were formerly solving the least squares problem here, but the solution to the L2 problem is just the mean of x - x_hat, so now I'm just doing that in the M=infinity case. However, having the minimization problem set up was helpful for seeing how to stick the Huber in there.
improved evaluation code to have robust root mean error metric and simpler rmse function. Amended how constants of integration are found to be robust. Optimizer loss function now uses x_hat directly instead of trying to integrate the derivative, since all algorithms return an x_hat estimate associated with dxdt_hat anyway. SpectralDiff now does an IFFT to get x_hat instead of trapezoidal integration. Improved a few docstrings. Got tests to pass with robust change. Nixed a few tests of the optimizer, because they're redundant with each other and the notebooks