Merge pull request #111 from florisvb/kf-P0

addressing #110

Author: pavelkomarov

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)
 

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):

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