diff --git a/pynumdiff/kalman_smooth/_kalman_smooth.py b/pynumdiff/kalman_smooth/_kalman_smooth.py index febd68b..c27e618 100644 --- a/pynumdiff/kalman_smooth/_kalman_smooth.py +++ b/pynumdiff/kalman_smooth/_kalman_smooth.py @@ -84,7 +84,7 @@ def _constant_derivative(x, P0, A, C, R, Q, forwardbackward): """Helper for `constant_{velocity,acceleration,jerk}` functions, because there was a lot of repeated code. """ - xhat0 = np.zeros(A.shape[0]); xhat0[0] = x[0] + xhat0 = np.zeros(A.shape[0]); xhat0[0] = x[0] # See #110 for why this choice of xhat0 xhat_smooth = _RTS_smooth(xhat0, P0, x, A, C, Q, R) # noisy x are the "y" in Kalman-land x_hat_forward = xhat_smooth[:, 0] # first dimension is time, so slice first element at all times dxdt_hat_forward = xhat_smooth[:, 1] @@ -135,7 +135,7 @@ def constant_velocity(x, dt, params=None, options=None, r=None, q=None, forwardb C = np.array([[1, 0]]) # we measure only y = noisy x R = np.array([[r]]) Q = np.array([[1e-16, 0], [0, q]]) # uncertainty is around the velocity - P0 = np.array(100*np.eye(2)) # Why is this one magnitude 100 vs the other ones being 10? + P0 = np.array(100*np.eye(2)) # See #110 for why this choice of P0 return _constant_derivative(x, P0, A, C, R, Q, forwardbackward) @@ -174,7 +174,7 @@ def constant_acceleration(x, dt, params=None, options=None, r=None, q=None, forw Q = np.array([[1e-16, 0, 0], [0, 1e-16, 0], [0, 0, q]]) # uncertainty is around the acceleration - P0 = np.array(10*np.eye(3)) + P0 = np.array(100*np.eye(3)) # See #110 for why this choice of P0 return _constant_derivative(x, P0, A, C, R, Q, forwardbackward) @@ -215,7 +215,7 @@ def constant_jerk(x, dt, params=None, options=None, r=None, q=None, forwardbackw [0, 1e-16, 0, 0], [0, 0, 1e-16, 0], [0, 0, 0, q]]) # uncertainty is around the jerk - P0 = np.array(10*np.eye(4)) + P0 = np.array(100*np.eye(4)) # See #110 for why this choice of P0 return _constant_derivative(x, P0, A, C, R, Q, forwardbackward) diff --git a/pynumdiff/tests/conftest.py b/pynumdiff/tests/conftest.py index 08825eb..e36dbca 100644 --- a/pynumdiff/tests/conftest.py +++ b/pynumdiff/tests/conftest.py @@ -3,7 +3,9 @@ from matplotlib import pyplot from collections import defaultdict -def pytest_addoption(parser): parser.addoption("--plot", action="store_true", default=False) +def pytest_addoption(parser): + parser.addoption("--plot", action="store_true", default=False) # whether to show plots + parser.addoption("--bounds", action="store_true", default=False) # whether to print error bounds @pytest.fixture(scope="session", autouse=True) def store_plots(request): diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index d56fff5..2f142d9 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -48,6 +48,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) (constant_jerk, {'r':1e-4, 'q':10}), (constant_jerk, [1e-4, 10]), # TODO (known_dynamics), but presently it doesn't calculate a derivative ] +diff_methods_and_params = [(constant_jerk, {'r':1e-4, 'q':10})] # All the testing methodology follows the exact same pattern; the only thing that changes is the # closeness to the right answer various methods achieve with the given parameterizations. So index a @@ -133,14 +134,14 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(1, 1), (2, 2), (1, 1), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], constant_acceleration: [[(-25, -25), (-25, -25), (0, 0), (1, 0)], - [(-4, -4), (-2, -3), (0, 0), (1, 0)], - [(-3, -3), (-1, -2), (0, 0), (1, 0)], + [(-5, -5), (-3, -3), (0, 0), (1, 0)], + [(-4, -4), (-2, -2), (0, 0), (1, 0)], [(0, -1), (1, 0), (0, 0), (1, 0)], [(1, 1), (3, 2), (1, 1), (3, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], constant_jerk: [[(-25, -25), (-25, -25), (0, 0), (1, 0)], - [(-4, -4), (-2, -3), (0, 0), (1, 0)], - [(-3, -3), (-1, -2), (0, 0), (1, 0)], + [(-5, -5), (-3, -3), (0, 0), (1, 0)], + [(-4, -4), (-2, -2), (0, 0), (1, 0)], [(-1, -2), (1, 0), (0, 0), (1, 0)], [(1, 0), (2, 2), (1, 0), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]] @@ -174,18 +175,19 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re else diff_method(x_noisy, dt, params, options) # check x_hat and x_hat_noisy are close to x and that dxdt_hat and dxdt_hat_noisy are close to dxdt - #print("]\n[", end="") + if request.config.getoption("--bounds"): print("]\n[", end="") for j,(a,b) in enumerate([(x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy)]): - # l2_error = np.linalg.norm(a - b) - # linf_error = np.max(np.abs(a - b)) - # print(f"({l2_error},{linf_error})", end=", ") - # print(f"({int(np.ceil(np.log10(l2_error))) if l2_error > 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ") - - log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j] - assert np.linalg.norm(a - b) < 10**log_l2_bound - assert np.max(np.abs(a - b)) < 10**log_linf_bound - if 0 < np.linalg.norm(a - b) < 10**(log_l2_bound - 1) or 0 < np.max(np.abs(a - b)) < 10**(log_linf_bound - 1): - print(f"Improvement detected for method {diff_method.__name__}") + if request.config.getoption("--bounds"): + l2_error = np.linalg.norm(a - b) + linf_error = np.max(np.abs(a - b)) + #print(f"({l2_error},{linf_error})", end=", ") + print(f"({int(np.ceil(np.log10(l2_error))) if l2_error > 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ") + else: + log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j] + assert np.linalg.norm(a - b) < 10**log_l2_bound + assert np.max(np.abs(a - b)) < 10**log_linf_bound + if 0 < np.linalg.norm(a - b) < 10**(log_l2_bound - 1) or 0 < np.max(np.abs(a - b)) < 10**(log_linf_bound - 1): + print(f"Improvement detected for method {diff_method.__name__}") if request.config.getoption("--plot") and not isinstance(params, list): # Get the plot flag from pytest configuration fig, axes = request.config.plots[diff_method] # get the appropriate plot, set up by the store_plots fixture in conftest.py