last touch ups

Author: pavelkomarov

pynumdiff/optimize/_optimize.py

@@ -66,7 +66,7 @@
                 'lmbd': (1e-3, 0.5)}),
     tvrdiff: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
                '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
-              'huberM': [0., 1, 2, 3, 6]},
+              'huberM': [0., 1, 2, 6]}, # comb lower values more finely, because the scale of sigma is mad(x), bigger than mad(y-x) residuals
               {'gamma': (1e-4, 1e7),
               'huberM': (0, 6)}),
     velocity: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]}, # Deprecated method

pynumdiff/total_variation_regularization/_total_variation_regularization.py

@@ -89,7 +89,7 @@ def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None):
         y = cvxpy.cumsum(y) + integration_constants[i]
 
     # Compare the recursively integrated position to the noisy position. \ell_2 doesn't get scaled by 1/2 here,
-    # so cvxpy Huber is already the right scale, and \ell_1 should be scaled by 2\sqrt{2} to match.
+    # so cvxpy's doubled Huber is already the right scale, and \ell_1 should be scaled by 2\sqrt{2} to match.
     fidelity_cost = cvxpy.sum_squares(y - x) if huberM == float('inf') \
             else np.sqrt(8)*cvxpy.norm(y - x, 1) if huberM == 0 \
             else utility.huber_const(huberM)*cvxpy.sum(cvxpy.huber(y - x, huberM*sigma))

pynumdiff/utils/utility.py

@@ -138,9 +138,9 @@ def estimate_integration_constant(x, x_hat, M=6):
     :return: **integration constant** (float) -- initial condition that best aligns x_hat with x
     """
     sigma = median_abs_deviation(x - x_hat, scale='normal') # M is in units of this robust scatter metric
-    if M == float('inf') or sigma < 1e-6: # If no scatter, then no outliers, so use L2
+    if M == float('inf') or sigma < 1e-3: # If no scatter, then no outliers, so use L2
         return np.mean(x - x_hat) # Solves the L2 distance minimization, argmin_{x0} ||x_hat + x0 - x||_2^
-    elif M < 1e-2: # small M looks like L1 loss, and Huber gets too flat to work well
+    elif M < 1e-3: # small M looks like L1 loss, and Huber gets too flat to work well
         return np.median(x - x_hat) # Solves the L1 distance minimization
     else:
         return minimize(lambda x0: np.sum(huber(x - (x_hat+x0), M*sigma)), # fn to minimize in 1st argument