moved TVR to new optimizer

Author: pavelkomarov

pynumdiff/optimize/_optimize.py

@@ -10,6 +10,8 @@
 from ..finite_difference import first_order
 from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
 from ..linear_model import spectraldiff, polydiff, savgoldiff, lineardiff
+from ..total_variation_regularization import velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding
+from ..kalman_smooth import constant_velocity, constant_acceleration, constant_jerk
 
 
 # Map from method -> (init_conds, bounds_low_hi)
@@ -59,9 +61,27 @@
                  {'order': (3, 5),
                   's': (1e-2, 1e6),
                   'num_iterations': (1, 10)}),
+    velocity: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]},
+               {'gamma': (1e-4, 1e7)}),
+    iterative_velocity: ({'num_iterations': [1, 5, 10],
+                          'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
+                          'scale': 'small'},
+                         {'num_iterations': (1, 100), # gets expensive with more iterations
+                          'gamma': (1e-4, 1e7)})
+    smooth_acceleration: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
+                           'window_size': [3, 10, 30, 50, 90, 130]},
+                          {'gamma': (1e-4, 1e7),
+                           'window_size': (1, 1000)})
+    # jerk_sliding:
+    # constant_velocity:
+    # constant_acceleration:
+    # constant_jerk:
 }
 for method in [meandiff, gaussiandiff, friedrichsdiff]:
     method_params_and_bounds[method] = method_params_and_bounds[mediandiff]
+for method in [acceleration, jerk]:
+    method_params_and_bounds[method] = method_params_and_bounds[velocity]
+method_params_and_bounds[jerk_sliding] = method_params_and_bounds[smooth_acceleration]
 
 
 # This function to be at the top level for multiprocessing
@@ -97,7 +117,7 @@ def _objective_function(point, func, x, dt, singleton_params, search_space_types
 
 
 def optimize(func, x, dt, init_conds={}, dxdt_truth=None, tvgamma=1e-2, padding='auto', metric='rmse',
-    opt_method='Nelder-Mead', opt_kwargs={'maxiter': 10}):
+    opt_method='Nelder-Mead', maxiter=10):
     """Find the optimal parameters for a given differentiation method.
 
     :param function func: differentiation method to optimize parameters for, e.g. linear_model.savgoldiff
@@ -115,7 +135,7 @@ def optimize(func, x, dt, init_conds={}, dxdt_truth=None, tvgamma=1e-2, padding=
     :param str metric: either :code:`'rmse'` or :code:`'error_correlation'`, only applies if :code:`dxdt_truth`
                     is not None, see _objective_function
     :param str opt_method: Optimization technique used by :code:`scipy.minimize`, the workhorse
-    :param dict opt_kwargs: keyword arguments to pass down to :code:`scipy.minimize`
+    :param int maxiter: passed down to :code:`scipy.minimize`, maximum iterations
 
     :return: tuple[dict, float] of\n
             - **opt_params** -- best parameter settings for the differentation method
@@ -147,7 +167,7 @@ def optimize(func, x, dt, init_conds={}, dxdt_truth=None, tvgamma=1e-2, padding=
     _obj_fun = partial(_objective_function, func=func, x=x, dt=dt, singleton_params=singleton_params,
         search_space_types=search_space_types, dxdt_truth=dxdt_truth, metric=metric, tvgamma=tvgamma,
         padding=padding)
-    _minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options=opt_kwargs)
+    _minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options={'maxiter':maxiter})
 
     with Pool(initializer=filterwarnings, initargs=["ignore", '', UserWarning]) as pool: # The heavy lifting
         results = pool.map(_minimize, search_space) # returns a bunch of OptimizeResult objects

pynumdiff/optimize/total_variation_regularization/__init__.py

@@ -1,15 +1,12 @@
 import numpy as np
 from pytest import mark, skip
 
-#from pynumdiff.optimize.total_variation_regularization import *
-#from pynumdiff.optimize.kalman_smooth import constant_velocity, constant_acceleration, \
-#    constant_jerk
-from pynumdiff.utils.simulate import pi_control
-
 from ..finite_difference import first_order as iterated_finite_difference # actually second order
 from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
 from ..linear_model import spectraldiff, polydiff, savgoldiff
+from ..total_variation_regularization import velocity, acceleration, iterative_velocity
 from ..optimize import optimize
+from ..utils.simulate import pi_control
 
 
 # simulation
@@ -51,62 +48,49 @@ def test_friedrichsdiff():
     assert params2['window_size'] == 1
 
 def test_butterdiff():
-    params1, val1 = optimize(butterdiff, x, dt, init_conds={'num_iterations':1}, tvgamma=tvgamma, dxdt_truth=dxdt_truth, opt_kwargs={'maxiter': 20})
-    params2, val2 = optimize(butterdiff, x, dt, init_conds={'num_iterations':1}, tvgamma=0, dxdt_truth=None, opt_kwargs={'maxiter': 20})
+    params1, val1 = optimize(butterdiff, x, dt, init_conds={'num_iterations':1}, tvgamma=tvgamma, dxdt_truth=dxdt_truth, maxiter=20)
+    params2, val2 = optimize(butterdiff, x, dt, init_conds={'num_iterations':1}, tvgamma=0, dxdt_truth=None, maxiter=20)
 
     assert params1['filter_order'] == 8
     np.testing.assert_almost_equal(params1['cutoff_freq'], 0.161, decimal=3)
     assert params2['filter_order'] == 3
     np.testing.assert_almost_equal(params2['cutoff_freq'], 0.99, decimal=3)
 
 def test_splinediff():
-    params1, val1 = optimize(splinediff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth, opt_kwargs={'maxiter': 20})
-    params2, val2 = optimize(splinediff, x, dt, tvgamma=0, dxdt_truth=None, opt_kwargs={'maxiter': 20})
+    params1, val1 = optimize(splinediff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth, maxiter=20)
+    params2, val2 = optimize(splinediff, x, dt, tvgamma=0, dxdt_truth=None, maxiter=20)
 
     assert (params1['order'], params1['num_iterations']) == (4, 1)
     np.testing.assert_almost_equal(params1['s'], 0.0146, decimal=3)
     assert (params2['order'], params2['num_iterations']) == (4, 1)
     np.testing.assert_almost_equal(params2['s'], 0.01, decimal=3)
 
-# def test_iterative_velocity():
-#     params1, val1 = iterative_velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-#     params2, val2 = iterative_velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-#     np.testing.assert_array_less( np.abs(params1[0] - 2), 1.001)
-#     np.testing.assert_array_less( np.abs(params2[0] - 2), 1.001)
-    
-#     np.testing.assert_almost_equal(params1[1], 0.0001, decimal=4)
-#     np.testing.assert_almost_equal(params2[1], 0.0001, decimal=4)
+def test_iterative_velocity():
+    params1, val1 = optimize(iterative_velocity, x, dt, init_conds={'num_iterations':1}, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
+    params2, val2 = optimize(iterative_velocity, x, dt, init_conds={'num_iterations':1}, tvgamma=0, dxdt_truth=None)
 
-#     #self.assertListEqual(params1, [2, 0.0001])
-#     #self.assertListEqual(params2, [2, 0.0001])
+    np.testing.assert_almost_equal(params1['gamma'], 0.0001, decimal=4)
+    np.testing.assert_almost_equal(params2['gamma'], 0.0001, decimal=4)
 
-# def test_velocity():
-#     try:
-#         import cvxpy
-#     except:
-#         pytest.skip("could not import cvxpy, skipping test_velocity", allow_module_level=True)
+def test_velocity():
+    try: import cvxpy
+    except: skip("could not import cvxpy, skipping test_velocity")
 
-#     params1, val1 = velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-#     params2, val2 = velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-#     param1_error = np.abs(params1[0] - 0.07218)
-#     param2_error = np.abs(params2[0] - 0.0001)
+    params1, val1 = optimize(velocity, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth, maxiter=20)
+    params2, val2 = optimize(velocity, x, dt, tvgamma=0, dxdt_truth=None, maxiter=20)
 
-#     np.testing.assert_array_less(param1_error, 2)
-#     np.testing.assert_array_less(param2_error, 2)
+    np.testing.assert_almost_equal(params1['gamma'], 0.0721, decimal=3)
+    np.testing.assert_almost_equal(params2['gamma'], 1e-4, decimal=4)
 
-# def test_acceleration():
-#     try:
-#         import cvxpy
-#     except:
-#         pytest.skip("could not import cvxpy, skipping test_acceleration", allow_module_level=True)
+def test_acceleration():
+    try: import cvxpy
+    except: pytest.skip("could not import cvxpy, skipping test_acceleration")
 
-#     params1, val1 = acceleration(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-#     params2, val2 = acceleration(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-#     param1_error = np.abs(params1[0] - 0.1447)
-#     param2_error = np.abs(params2[0] - 0.0001)
+    params1, val1 = optimize(acceleration, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth, maxiter=20)
+    params2, val2 = optimize(acceleration, x, dt, tvgamma=0, dxdt_truth=None, maxiter=20)
 
-#     np.testing.assert_array_less(param1_error, 2)
-#     np.testing.assert_array_less(param2_error, 2)
+    np.testing.assert_almost_equal(params1['gamma'], 0.144, decimal=3)
+    np.testing.assert_almost_equal(params2['gamma'], 1e-4, decimal=4)
 
 def test_savgoldiff():
     params1, val1 = optimize(savgoldiff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
@@ -127,12 +111,7 @@ def test_polydiff():
     assert (params2['poly_order'], params2['window_size']) == (4, 10)
 
 # def test_chebydiff(self):
-#     try:
-#         import pychebfun
-#     except:
-#         pytest.skip("could not import pychebfun, skipping test_chebydiff", allow_module_level=True)
-
-#     params1, val1 = chebydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-#     params2, val2 = chebydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-#     self.assertListEqual(params1, [9, 108])
-#     self.assertListEqual(params2, [9, 94])
+#     params1, val1 = optimize(chebydiff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
+#     params2, val2 = optimize(chebydiff, x, dt, tvgamma=0, dxdt_truth=None)
+#     assert params1 == [9, 108]
+#     assert params2 == [9, 94]