moved rest of linear_model optimization

Author: pavelkomarov

pynumdiff/linear_model/_linear_model.py

@@ -13,16 +13,16 @@
 #########################
 # Savitzky-Golay filter #
 #########################
-def savgoldiff(x, dt, params=None, options=None, polynomial_order=None, window_size=None, smoothing_win=None):
+def savgoldiff(x, dt, params=None, options=None, poly_order=None, window_size=None, smoothing_win=None):
     """Use the Savitzky-Golay to smooth the data and calculate the first derivative. It wses
     scipy.signal.savgol_filter. The Savitzky-Golay is very similar to the sliding polynomial fit,
     but slightly noisier, and much faster
 
     :param np.array[float] x: data to differentiate
     :param float dt: step size
-    :param list params: (**deprecated**, prefer :code:`polynomial_order`, :code:`window_size`, and :code:`smoothing_win`)
+    :param list params: (**deprecated**, prefer :code:`poly_order`, :code:`window_size`, and :code:`smoothing_win`)
     :param dict options: (**deprecated**)
-    :param int polynomial_order: order of the polynomial
+    :param int poly_order: order of the polynomial
     :param int window_size: size of the sliding window, must be odd (if not, 1 is added)
     :param int smoothing_win: size of the window used for gaussian smoothing, a good default is
         window_size, but smaller for high frequnecy data
@@ -32,17 +32,17 @@ def savgoldiff(x, dt, params=None, options=None, polynomial_order=None, window_s
              - **dxdt_hat** -- estimated derivative of x
     """
     if params != None: # Warning to support old interface for a while. Remove these lines along with params in a future release.
-        warn("`params` and `options` parameters will be removed in a future version. Use `polynomial_order`, " +
+        warn("`params` and `options` parameters will be removed in a future version. Use `poly_order`, " +
             "`window_size`, and `smoothing_win` instead.", DeprecationWarning)
-        polynomial_order, window_size, smoothing_win = params
-    elif polynomial_order == None or window_size == None or smoothing_win == None:
-        raise ValueError("`polynomial_order`, `window_size`, and `smoothing_win` must be given.")
+        poly_order, window_size, smoothing_win = params
+    elif poly_order == None or window_size == None or smoothing_win == None:
+        raise ValueError("`poly_order`, `window_size`, and `smoothing_win` must be given.")
 
-    window_size = np.clip(window_size, polynomial_order + 1, len(x)-1)
+    window_size = np.clip(window_size, poly_order + 1, len(x)-1)
     if not window_size % 2: window_size += 1 # window_size needs to be odd
     smoothing_win = min(smoothing_win, len(x)-1)
 
-    dxdt_hat = scipy.signal.savgol_filter(x, window_size, polynomial_order, deriv=1)/dt
+    dxdt_hat = scipy.signal.savgol_filter(x, window_size, poly_order, deriv=1)/dt
 
     kernel = utility.gaussian_kernel(smoothing_win)
     dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel)
@@ -57,16 +57,16 @@ def savgoldiff(x, dt, params=None, options=None, polynomial_order=None, window_s
 ######################
 # Polynomial fitting #
 ######################
-def polydiff(x, dt, params=None, options=None, polynomial_order=None, window_size=None, step_size=1,
+def polydiff(x, dt, params=None, options=None, poly_order=None, window_size=None, step_size=1,
     kernel='friedrichs'):
     """Fit polynomials to the data, and differentiate the polynomials.
 
     :param np.array[float] x: data to differentiate
     :param float dt: step size
-    :param list[int] params: (**deprecated**, prefer :code:`polynomial_order` and :code:`window_size`)
+    :param list[int] params: (**deprecated**, prefer :code:`poly_order` and :code:`window_size`)
     :param dict options: (**deprecated**, prefer :code:`step_size` and :code:`kernel`)
             a dictionary consisting of {'sliding': (bool), 'step_size': (int), 'kernel_name': (str)}
-    :param int polynomial_order: order of the polynomial
+    :param int poly_order: order of the polynomial
     :param int window_size: size of the sliding window, if not given no sliding
     :param int step_size: step size for sliding
     :param str kernel: name of kernel to use for weighting and smoothing windows ('gaussian' or 'friedrichs')
@@ -76,24 +76,27 @@ def polydiff(x, dt, params=None, options=None, polynomial_order=None, window_siz
              - **dxdt_hat** -- estimated derivative of x
     """
     if params != None:
-        warn("`params` and `options` parameters will be removed in a future version. Use `polynomial_order` " +
+        warn("`params` and `options` parameters will be removed in a future version. Use `poly_order` " +
             "and `window_size` instead.", DeprecationWarning)
-        polynomial_order = params[0]
+        poly_order = params[0]
         if len(params) > 1: window_size = params[1]
         if options != None:
             if 'sliding' in options and not options['sliding']: window_size = None
             if 'step_size' in options: step_size = options['step_size']
             if 'kernel_name' in options: kernel = options['kernel_name']
-    elif polynomial_order == None or window_size == None:
-        raise ValueError("`polynomial_order` and `window_size` must be given.")
+    elif poly_order == None or window_size == None:
+        raise ValueError("`poly_order` and `window_size` must be given.")
 
-    if window_size < polynomial_order*3:
-        window_size = polynomial_order*3+1
+    if window_size < poly_order*3:
+        window_size = poly_order*3+1
+    if window_size % 2 == 0:
+        window_size += 1
+        warn("Kernel window size should be odd. Added 1 to length.")
 
-    def _polydiff(x, dt, polynomial_order, weights=None):
+    def _polydiff(x, dt, poly_order, weights=None):
         t = np.arange(len(x))*dt
 
-        r = np.polyfit(t, x, polynomial_order, w=weights) # polyfit returns highest order first
+        r = np.polyfit(t, x, poly_order, w=weights) # polyfit returns highest order first
         dr = np.polyder(r) # power rule already implemented for us
 
         dxdt_hat = np.polyval(dr, t) # evaluate the derivative and original polynomials at points t
@@ -102,10 +105,10 @@ def _polydiff(x, dt, polynomial_order, weights=None):
         return x_hat, dxdt_hat
 
     if not window_size:
-        return _polydiff(x, dt, polynomial_order)
+        return _polydiff(x, dt, poly_order)
 
     kernel = {'gaussian':utility.gaussian_kernel, 'friedrichs':utility.friedrichs_kernel}[kernel](window_size)
-    return utility.slide_function(_polydiff, x, dt, kernel, polynomial_order, stride=step_size, pass_weights=True)
+    return utility.slide_function(_polydiff, x, dt, kernel, poly_order, stride=step_size, pass_weights=True)
 
 
 #############

pynumdiff/optimize/_optimize.py

@@ -3,26 +3,43 @@
 import numpy as np
 from itertools import product
 from functools import partial
+from warnings import filterwarnings
 from multiprocessing import Pool
 
 from pynumdiff.utils import evaluate
 
-from ..linear_model import spectraldiff, polydiff
+from ..linear_model import spectraldiff, polydiff, savgoldiff, lineardiff
 
 
 # Map from method -> (init_conds, type_low_hi)
 method_params_and_bounds = {
     spectraldiff: ({'even_extension': True,
                    'pad_to_zero_dxdt': True,
                    'high_freq_cutoff': [1e-3, 5e-2, 1e-2, 5e-2, 1e-1]},
-                  {'high_freq_cutoff': [1e-5, 1-1e-5]}),
-    polydiff: ({'sliding': True,
-                'step_size': 1,
+                  {'high_freq_cutoff': (1e-5, 1-1e-5)}),
+    polydiff: ({'step_size': 1,
                 'kernel': 'friedrichs',
-                'order': [2, 3, 5, 7],
-                'window_size': [10, 30, 50, 90, 130]},
-               {'order': [1, 8],
-                'window_size': [10, 1000]})
+                'poly_order': [2, 3, 5, 7],
+                'window_size': [11, 31, 51, 91, 131]},
+               {'poly_order': (1, 8),
+                'window_size': (10, 1000)}),
+    savgoldiff: ({'poly_order': [2, 3, 5, 7, 9, 11, 13],
+                  'window_size': [3, 10, 30, 50, 90, 130, 200, 300],
+                  'smoothing_win': [3, 10, 30, 50, 90, 130, 200, 300]},
+                 {'poly_order': (1, 12),
+                  'window_size': (3, 1000),
+                  'smoothing_win': (3, 1000)}),
+    #chebydiff ({order: [3, 5, 7, 9],
+    #           window_size: [10, 30, 50, 90, 130],
+    #           kernel: 'friedrichs'},
+    #           {order: (1, 10),
+    #           window_size: (10, 1000)})
+    lineardiff: ({'kernel': 'gaussian',
+                  'order': 3,
+                  'gamma': [1e-1, 1, 10, 100],
+                  'window_size': [10, 30, 50, 90, 130]},
+                 {'gamma': (1e-3, 1000),
+                  'window_size': (15, 1000)}),
 }
 
 
@@ -110,7 +127,7 @@ def optimize(func, x, dt, init_conds={}, dxdt_truth=None, tvgamma=1e-2, padding=
         padding=padding)
     _minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options=opt_kwargs)
 
-    with Pool() as pool: # The heavy lifting
+    with Pool(initializer=filterwarnings, initargs=["ignore", '', UserWarning]) as pool: # The heavy lifting
         results = pool.map(_minimize, search_space) # returns a bunch of OptimizeResult objects
 
     opt_idx = np.nanargmin([r.fun for r in results])

pynumdiff/optimize/linear_model/__init__.py

@@ -34,8 +34,8 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
     (iterated_first_order, {'num_iterations':5}), (iterated_first_order, [5], {'iterate':True}),
     (second_order, {}),
     (lineardiff, {'order':3, 'gamma':5, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 5, 11], {'solver':'CLARABEL'}),
-    (polydiff, {'polynomial_order':2, 'window_size':3}), (polydiff, [2, 3]),
-    (savgoldiff, {'polynomial_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]),
+    (polydiff, {'poly_order':2, 'window_size':3}), (polydiff, [2, 3]),
+    (savgoldiff, {'poly_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]),
     (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]),
     (mediandiff, {'window_size':3, 'num_iterations':2}), (mediandiff, [3, 2], {'iterate':True}),
     (meandiff, {'window_size':3, 'num_iterations':2}), (meandiff, [3, 2], {'iterate':True}),

pynumdiff/optimize/linear_model/__linear_model__.py

@@ -1,9 +1,5 @@
-"""
-Unit tests for optimization module
-"""
-# pylint: skip-file
 import numpy as np
-import pytest 
+from pytest import mark, skip
 
 #from pynumdiff.optimize.finite_difference import first_order
 #from pynumdiff.optimize.smooth_finite_difference import mediandiff, meandiff, gaussiandiff, \
@@ -14,7 +10,7 @@
 #    constant_jerk
 from pynumdiff.utils.simulate import pi_control
 
-from ..linear_model import spectraldiff, polydiff
+from ..linear_model import spectraldiff, polydiff, savgoldiff
 from ..optimize import optimize
 
 
@@ -29,10 +25,6 @@
 log_gamma = -1.6 * np.log(cutoff_frequency) - 0.71 * np.log(dt) - 5.1
 tvgamma = np.exp(log_gamma)
 
-def get_err_msg(actual_params, desired_params):
-    err_msg = 'Actual params were: ' + ', '.join(map(str, actual_params)) + ' instead of: ' + ', '.join(map(str, desired_params))
-    return err_msg
-
 
 # def test_first_order():
 #     params_1, val_1 = first_order(x, dt, params=None, options={'iterate': True},
@@ -78,10 +70,10 @@ def get_err_msg(actual_params, desired_params):
 #     params_1, val_1 = butterdiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
 #     params_2, val_2 = butterdiff(x, dt, params=None, tvgamma=0, dxdt_truth=None)
 
-#     np.testing.assert_array_less( np.abs(params_1[0] - 9), 1.001, err_msg=get_err_msg(params_1, [9, 0.157]))
-#     np.testing.assert_array_less( np.abs(params_1[1] - 0.157), 0.01, err_msg=get_err_msg(params_1, [9, 0.157]))
-#     #np.testing.assert_almost_equal(params_1, [9, 0.157], decimal=3, err_msg=get_err_msg(params_1, [9, 0.157]))
-#     np.testing.assert_almost_equal(params_2, [3, 0.99], decimal=3, err_msg=get_err_msg(params_2, [3, 0.99]))
+#     np.testing.assert_array_less( np.abs(params_1[0] - 9), 1.001)
+#     np.testing.assert_array_less( np.abs(params_1[1] - 0.157), 0.01)
+#     #np.testing.assert_almost_equal(params_1, [9, 0.157], decimal=3)
+#     np.testing.assert_almost_equal(params_2, [3, 0.99], decimal=3)
 
 # def test_splinediff():
 #     params_1, val_1 = splinediff(x, dt, params=None, options={'iterate': True},
@@ -131,23 +123,23 @@ def get_err_msg(actual_params, desired_params):
 #     np.testing.assert_array_less(param_1_error, 2)
 #     np.testing.assert_array_less(param_2_error, 2)
 
-# def test_savgoldiff():
-#     params_1, val_1 = savgoldiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-#     params_2, val_2 = savgoldiff(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-#     assert params_1 == [10, 59, 3]
-#     assert params_2 == [9, 3, 3]
+def test_savgoldiff():
+    params_1, val_1 = optimize(savgoldiff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
+    params_2, val_2 = optimize(savgoldiff, x, dt, tvgamma=0, dxdt_truth=None)
+    assert (params_1['poly_order'], params_1['window_size'], params_1['smoothing_win']) == (10, 57, 3)
+    assert (params_2['poly_order'], params_2['window_size'], params_2['smoothing_win']) == (9, 4, 3)
 
 def test_spectraldiff():
     params_1, val_1 = optimize(spectraldiff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
     params_2, val_2 = optimize(spectraldiff, x, dt, tvgamma=0)
     np.testing.assert_almost_equal(params_1['high_freq_cutoff'], 0.0913, decimal=3)
     np.testing.assert_almost_equal(params_2['high_freq_cutoff'], 0.575, decimal=3)
 
-# def test_polydiff():
-    # params_1, val_1 = polydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
-    # params_2, val_2 = polydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None)
-    # assert params_1 == [2, 10]
-    # assert params_2 == [2, 10]
+def test_polydiff():
+    params_1, val_1 = optimize(polydiff, x, dt, tvgamma=tvgamma, dxdt_truth=dxdt_truth)
+    params_2, val_2 = optimize(polydiff, x, dt, tvgamma=0, dxdt_truth=None)
+    assert (params_1['poly_order'], params_1['window_size']) == (6, 50)
+    assert (params_2['poly_order'], params_2['window_size']) == (4, 10)
 
 # def test_chebydiff(self):
 #     try: