18988228033

Committed 01 Nov 2025 12:12AM UTC coverage: 72.783% (-1.9%) from 74.718%

Build # 18988228033

Build Type

push

github

Committed by

pavelkomarov

Commit Message

trying to fix weird issue where cvxopt isn't installing properly

Run Details

837 of 1150 relevant lines covered (72.78%)

0.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

57.14

/pynumdiff/optimize/_optimize.py

"""Optimization"""
import scipy.optimize
import numpy as np
from itertools import product
from functools import partial
from warnings import filterwarnings, warn
from multiprocessing import Pool, Manager
from hashlib import sha1
from tqdm import tqdm

from ..utils import evaluate
from ..finite_difference import finitediff, first_order, second_order, fourth_order
from ..smooth_finite_difference import kerneldiff, mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff
from ..polynomial_fit import polydiff, savgoldiff, splinediff
from ..basis_fit import spectraldiff, rbfdiff
from ..total_variation_regularization import tvrdiff, velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding
from ..kalman_smooth import rtsdiff, constant_velocity, constant_acceleration, constant_jerk, robustdiff
from ..linear_model import lineardiff

# Map from method -> (search_space, bounds_low_hi)
method_params_and_bounds = {
    kerneldiff: ({'kernel': {'mean', 'median', 'gaussian', 'friedrichs'},
             'window_size': [5, 15, 30, 50],
          'num_iterations': [1, 5, 10]},
            {'window_size': (1, 1e6),
          'num_iterations': (1, 100)}),
    meandiff: ({'window_size': [5, 15, 30, 50], # Deprecated method
             'num_iterations': [1, 5, 10]},
               {'window_size': (1, 1e6),
             'num_iterations': (1, 100)}),
    butterdiff: ({'filter_order': set(i for i in range(1,11)), # categorical to save us from doing double work by guessing between orders
                   'cutoff_freq': [0.0001, 0.001, 0.005, 0.01, 0.1, 0.5],
                'num_iterations': [1, 5, 10]},
                  {'cutoff_freq': (1e-4, 1-1e-2),
                'num_iterations': (1, 1000)}),
    finitediff: ({'num_iterations': [5, 10, 30, 50],
                           'order': {2, 4}}, # order is categorical here, because it can't be 3
                 {'num_iterations': (1, 1000)}),
    first_order: ({'num_iterations': [5, 10, 30, 50]}, # Separated because optimizing over this one is rare due to shifted answer
                  {'num_iterations': (1, 1000)}),
    polydiff: ({'step_size': [1, 2, 5],
                   'kernel': {'friedrichs', 'gaussian'}, # categorical
                   'degree': [2, 3, 5, 7],
              'window_size': [11, 31, 51, 91, 131]},
               {'step_size': (1, 100),
                   'degree': (1, 8),
              'window_size': (10, 1000)}),
    savgoldiff: ({'degree': [2, 3, 5, 7, 10],
             'window_size': [3, 10, 30, 50, 90, 130, 200, 300],
           'smoothing_win': [3, 10, 30, 50, 90, 130, 200, 300]},
                 {'degree': (1, 12),
             'window_size': (3, 1000),
           'smoothing_win': (3, 1000)}),
    splinediff: ({'degree': {3, 4, 5}, # categorical, because degree is whole number, and there aren't many choices
                       's': [0.5, 0.9, 0.95, 1, 10, 100],
          'num_iterations': [1, 5, 10]},
                      {'s': (1e-2, 1e6),
          'num_iterations': (1, 10)}),
    spectraldiff: ({'even_extension': {True, False}, # give categorical params in a set
                  'pad_to_zero_dxdt': {True, False},
                  'high_freq_cutoff': [1e-3, 5e-2, 1e-2, 5e-2, 1e-1]}, # give numerical params in a list to scipy.optimize over them
                 {'high_freq_cutoff': (1e-5, 1-1e-5)}),
    rbfdiff: ({'sigma': [1e-3, 1e-2, 1e-1, 1],
                'lmbd': [1e-3, 1e-2, 1e-1]},
              {'sigma': (1e-3, 1e3),
                'lmbd': (1e-4, 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
              {'gamma': (1e-4, 1e7)}),
    velocity: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]}, # Deprecated method
               {'gamma': (1e-4, 1e7)}),
    iterative_velocity: ({'scale': 'small', # Rare to optimize this one, because it's longer-running than convex version
                 'num_iterations': [1, 5, 10],
                          'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]},
                {'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': (3, 1000)}),
    rtsdiff: ({'forwardbackward': {True, False},
                         'order': {1, 2, 3}, # for this few options, the optimization works better if this is categorical
                      'qr_ratio': [10**k for k in range(-9, 10, 2)] + [1e12, 1e16]},
                     {'qr_ratio': [1e-10, 1e20]}), # qr_ratio is usually >>1
    constant_velocity: ({'q': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8], # Deprecated method
                         'r': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8],
           'forwardbackward': {True, False}},
                        {'q': (1e-10, 1e10),
                         'r': (1e-10, 1e10)}),
    robustdiff: ({'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
                      'q': [1e-1, 1e1, 1e4, 1e8, 1e12],
                      'r': [1e-1, 1e1, 1e4, 1e8, 1e12],
            'proc_huberM': {0, 1, 1.345, 2, 6}, # 0 is l1 norm, 1.345 is Huber 95% "efficiency", 2 assumes about 5% outliers,
            'meas_huberM': {0, 1, 1.345, 2, 6}}, # and 6 assumes basically no outliers -> l2 norm. Try (1 - norm.cdf(M))*2 to see outlier portion
                     {'q': (1e-1, 1e18),
                      'r': (1e-5, 1e18),
                 'huberM': (0, 5)}), # really only want to use l2 norm when nearby
    lineardiff: ({'kernel': 'gaussian',
                   'order': 3,
                   'gamma': [1e-1, 1, 10, 100],
             'window_size': [10, 30, 50, 90, 130]},
                  {'order': (1, 5),
                   'gamma': (1e-3, 1000),
             'window_size': (15, 1000)})
} # Methods with nonunique parameter sets are aliased in the dictionary below
for method in [second_order, fourth_order]: # Deprecated, redundant methods
    method_params_and_bounds[method] = method_params_and_bounds[first_order]
for method in [mediandiff, gaussiandiff, friedrichsdiff]: # Deprecated methods
    method_params_and_bounds[method] = method_params_and_bounds[meandiff]
for method in [acceleration, jerk]: # Deprecated, redundant methods
    method_params_and_bounds[method] = method_params_and_bounds[velocity]
method_params_and_bounds[jerk_sliding] = method_params_and_bounds[smooth_acceleration]
for method in [constant_acceleration, constant_jerk]: # Deprecated, redundant methods
    method_params_and_bounds[method] = method_params_and_bounds[constant_velocity]


# This function has to be at the top level for multiprocessing but is only used by optimize.
def _objective_function(point, func, x, dt, singleton_params, categorical_params, search_space_types,
    dxdt_truth, metric, tvgamma, padding, cache):
    """Function minimized by scipy.optimize.minimize, needs to have the form: (point, *args) -> float
    This is mildly complicated, because "point" controls the settings of a differentiation function, but
    the method may have numerical and non-numerical parameters, and all such parameters are now passed by
    keyword arguments. So the encoded `point` has to be decoded to dict.

    :param np.array point: a numerical vector scipy chooses to try in the objective function
    :param dict singleton_params: maps parameter names to singleton values
    :param dict categorical_params: maps parameter names to values
    :param dict search_space_types: maps parameter names to types, for turning float vector point into dict point_params
    :param multiprocessing.manager.dict cache: available across processes to save results and work
    Other params documented in `optimize`

    :return: float, cost of this objective at the point
    """
    key = sha1((''.join(f"{v:.3e}" for v in point) + # This hash is stable across processes. Takes bytes
               ''.join(str(v) for k,v in sorted(categorical_params.items()))).encode()).digest()
    if key in cache: return cache[key] # short circuit if this hyperparam combo has already been queried, ~10% savings per #160

    point_params = {k:(v if search_space_types[k] == float else int(np.round(v)))
                        for k,v in zip(search_space_types, point)} # point -> dict
    # add back in the singletons we're not searching over
    try: x_hat, dxdt_hat = func(x, dt, **point_params, **singleton_params, **categorical_params) # estimate x and dxdt
    except np.linalg.LinAlgError: cache[key] = 1e10; return 1e10 # some methods can fail numerically

    # evaluate estimate
    if dxdt_truth is not None: # then minimize ||dxdt_hat - dxdt_truth||^2
        if metric == 'rmse':
            rms_rec_x, rms_x, rms_dxdt = evaluate.rmse(x, dt, x_hat, dxdt_hat, dxdt_truth=dxdt_truth, padding=padding)
            cache[key] = rms_dxdt; return rms_dxdt
        elif metric == 'error_correlation':
            ec = evaluate.error_correlation(dxdt_hat, dxdt_truth, padding=padding)
            cache[key] = ec; return ec
    else: # then minimize [ || integral(dxdt_hat) - x||^2 + gamma*TV(dxdt_hat) ]
        rms_rec_x, rms_x, rms_dxdt = evaluate.rmse(x, dt, x_hat, dxdt_hat, dxdt_truth=None, padding=padding)
        cost = rms_rec_x + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding)
        cache[key] = cost; return cost


def optimize(func, x, dt, dxdt_truth=None, tvgamma=1e-2, search_space_updates={}, metric='rmse',
    padding=0, opt_method='Nelder-Mead', maxiter=10, parallel=True):
    """Find the optimal parameters for a given differentiation method.

    :param function func: differentiation method to optimize parameters for, e.g. linear_model.savgoldiff
    :param np.array[float] x: data to differentiate
    :param float dt: step size
    :param np.array[float] dxdt_truth: actual time series of the derivative of x, if known
    :param float tvgamma: Only used if :code:`dxdt_truth` is given. Regularization value used to select for parameters
                    that yield a smooth derivative. Larger value results in a smoother derivative.
    :param dict search_space_updates: At the top of :code:`_optimize.py`, each method has a search space of parameters
                    settings structured as :code:`{param1:[values], param2:[values], param3:value, ...}`. The Cartesian
                    product of values are used as initial starting points in optimization. If left None, the default search
                    space is used.
    :param str metric: either :code:`'rmse'` or :code:`'error_correlation'`, only applies if :code:`dxdt_truth`
                    is not None, see _objective_function
    :param int padding: number of time steps to ignore at the beginning and end of the time series in the
                    optimization, or :code:`'auto'` to ignore 2.5% at each end. Larger value causes the
                    optimization to emphasize the accuracy of dxdt in the middle of the time series
    :param str opt_method: Optimization technique used by :code:`scipy.minimize`, the workhorse
    :param int maxiter: passed down to :code:`scipy.minimize`, maximum iterations
    :param bool parallel: whether to use multiple processes to optimize, typically faster for single optimizations, but
                    for experiments it is a better use of resources to pool at that higher level

    :return: tuple[dict, float] of\n
            - **opt_params** -- best parameter settings for the differentation method
            - **opt_value** -- lowest value found for objective function
    """
    if metric not in ['rmse','error_correlation']:
        raise ValueError('`metric` should either be `rmse` or `error_correlation`.')
    if metric == 'error_correlation' and dxdt_truth is None:
        raise ValueError('`metric` can only be `error_correlation` if `dxdt_truth` is given.')

    search_space, bounds = method_params_and_bounds[func]
    search_space.update(search_space_updates) # for things not given, use defaults

    # No need to optimize over singletons, just pass them through
    singleton_params = {k:v for k,v in search_space.items() if not isinstance(v, (list, set))}

    # To handle categoricals, find their combination, and then pass each set individually
    categorical_params = {k for k,v in search_space.items() if isinstance(v, set)}
    categorical_combos = [dict(zip(categorical_params, combo)) for combo in
        product(*[search_space[k] for k in categorical_params])] # ends up [{}] if there are no categorical params

    # The Nelder-Mead's search space is the dimensions where multiple numerical options are given in a list
    search_space_types = {k:type(v[0]) for k,v in search_space.items() if isinstance(v, list)} # map param name -> type, for converting to and from point
    if any(v not in [float, int] for v in search_space_types.values()):
        raise ValueError("To optimize over categorical strings or bools, put them in a tuple, not a list.")
    # Cast ints to floats, and we're good to go
    starting_points = list(product(*[np.array(search_space[k]).astype(float) for k in search_space_types]))
    # The numerical space should have bounds
    bounds = [bounds[k] if k in bounds else # pass these to minimize(). It should respect them.
            None for k,v in search_space_types.items()] # None means no bound on a dimension

    results = []
    filterwarnings("ignore", '', UserWarning) # An extra filtering call, because some worker work can actually be done in the main process
    if parallel:
        with Manager() as manager:
            cache = manager.dict() # cache answers to avoid expensive repeat queries
            with Pool(initializer=filterwarnings, initargs=["ignore", '', UserWarning]) as pool: # The heavy lifting
                for categorical_combo in categorical_combos:
                    # wrap the objective and scipy.optimize.minimize to pass kwargs and a host of other things that remain the same
                    _obj_fun = partial(_objective_function, func=func, x=x, dt=dt, singleton_params=singleton_params,
                        categorical_params=categorical_combo, search_space_types=search_space_types, dxdt_truth=dxdt_truth,
                        metric=metric, tvgamma=tvgamma, padding=padding, cache=cache)
                    _minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options={'maxiter':maxiter})
                    results += pool.map(_minimize, starting_points) # returns a bunch of OptimizeResult objects
    else: # For experiments, where I want to parallelize optimization calls and am not allowed to have each spawn further processes
        cache = dict()
        for categorical_combo in categorical_combos:
            _obj_fun = partial(_objective_function, func=func, x=x, dt=dt, singleton_params=singleton_params,
                categorical_params=categorical_combo, search_space_types=search_space_types, dxdt_truth=dxdt_truth,
                metric=metric, tvgamma=tvgamma, padding=padding, cache=cache)
            _minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options={'maxiter':maxiter})
            results = [_minimize(p) for p in starting_points]

    opt_idx = np.nanargmin([r.fun for r in results])
    opt_point = results[opt_idx].x
    # results are going to be floats, but that may not be allowed, so convert back to a dict
    opt_params = {k:(v if search_space_types[k] == float else 
                    int(np.round(v)) if search_space_types[k] == int else
                    v > 0.5) for k,v in zip(search_space_types, opt_point)}
    opt_params.update(singleton_params)
    opt_params.update(categorical_combos[opt_idx//len(starting_points)]) # there are |starting_points| results for each combo

    return opt_params, results[opt_idx].fun


def suggest_method(x, dt, dxdt_truth=None, cutoff_frequency=None):
    """This is meant as an easy-to-use, automatic way for users with some time on their hands to determine
    a good method and settings for their data. It calls the optimizer over (almost) all methods in the repo
    using default search spaces defined at the top of the :code:`pynumdiff/optimize/_optimize.py` file.
    This routine will take a few minutes to run, especially due to `robustdiff`.
    
    Excluded:
        - ``first_order``, because iterating causes drift
        - ``lineardiff``, ``iterative_velocity``, and ``jerk_sliding``, because they either take too long,
          can be fragile, or tend not to do best
        - all ``cvxpy``-based methods if it is not installed
        - ``velocity`` because it tends to not be best but dominates the optimization process by directly
          optimizing the second term of the metric :math:`L = \\text{RMSE} \\Big( \\text{trapz}(\\mathbf{
          \\hat{\\dot{x}}}(\\Phi)) + \\mu, \\mathbf{y} \\Big) + \\gamma \\Big({TV}\\big(\\mathbf{\\hat{
          \\dot{x}}}(\\Phi)\\big)\\Big)`

    :param np.array[float] x: data to differentiate
    :param float dt: step size, because most methods are not designed to work with variable step sizes
    :param np.array[float] dxdt_truth: if known, you can pass true derivative values; otherwise you must use
            :code: `cutoff_frequency`
    :param float cutoff_frequency: in Hz, the highest dominant frequency of interest in the signal,
            used to find parameter :math:`\\gamma` for regularization of the optimization process
            in the absence of ground truth. See https://ieeexplore.ieee.org/document/9241009.
            Estimate by (a) counting real number of peaks per second in the data, (b) looking at
            power spectrum and choosing a cutoff, or (c) making an educated guess.

    :return: tuple[callable, dict] of\n
            - **method** -- a reference to the function handle of the differentiation method that worked best
            - **opt_params** -- optimal parameter settings for the differentation method
    """
    tvgamma = None
    if dxdt_truth is None: # parameter checking
        if cutoff_frequency is None:
            raise ValueError('Either dxdt_truth or cutoff_frequency must be provided.')
        tvgamma = np.exp(-1.6*np.log(cutoff_frequency) -0.71*np.log(dt) - 5.1) # See https://ieeexplore.ieee.org/document/9241009

    methods = [meandiff, mediandiff, gaussiandiff, friedrichsdiff, butterdiff,
        polydiff, savgoldiff, splinediff, spectraldiff, rbfdiff, finitediff, rtsdiff]
    try: # optionally skip some methods
        import cvxpy
        methods += [tvrdiff, smooth_acceleration, robustdiff]
    except ImportError:
        warn("CVXPY not installed, skipping tvrdiff, smooth_acceleration, and robustdiff")

    best_value = float('inf') # core loop
    for func in tqdm(methods):
        p, v = optimize(func, x, dt, dxdt_truth=dxdt_truth, tvgamma=tvgamma, search_space_updates=(
            {'order':{2,3}} if func in [tvrdiff, robustdiff] else {})) # convex-based with order 1 hack the cost function
        if v < best_value:
            method = func
            best_value = v
            opt_params = p

    return method, opt_params

1	"""Optimization"""
2	import scipy.optimize	1✔
3	import numpy as np	1✔
4	from itertools import product	1✔
5	from functools import partial	1✔
6	from warnings import filterwarnings, warn	1✔
7	from multiprocessing import Pool, Manager	1✔
8	from hashlib import sha1	1✔
9	from tqdm import tqdm	1✔
10
11	from ..utils import evaluate	1✔
12	from ..finite_difference import finitediff, first_order, second_order, fourth_order	1✔
13	from ..smooth_finite_difference import kerneldiff, mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff	1✔
14	from ..polynomial_fit import polydiff, savgoldiff, splinediff	1✔
15	from ..basis_fit import spectraldiff, rbfdiff	1✔
16	from ..total_variation_regularization import tvrdiff, velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding	1✔
17	from ..kalman_smooth import rtsdiff, constant_velocity, constant_acceleration, constant_jerk, robustdiff	1✔
18	from ..linear_model import lineardiff	1✔
19
20	# Map from method -> (search_space, bounds_low_hi)
21	method_params_and_bounds = {	1✔
22	kerneldiff: ({'kernel': {'mean', 'median', 'gaussian', 'friedrichs'},
23	'window_size': [5, 15, 30, 50],
24	'num_iterations': [1, 5, 10]},
25	{'window_size': (1, 1e6),
26	'num_iterations': (1, 100)}),
27	meandiff: ({'window_size': [5, 15, 30, 50], # Deprecated method
28	'num_iterations': [1, 5, 10]},
29	{'window_size': (1, 1e6),
30	'num_iterations': (1, 100)}),
31	butterdiff: ({'filter_order': set(i for i in range(1,11)), # categorical to save us from doing double work by guessing between orders
32	'cutoff_freq': [0.0001, 0.001, 0.005, 0.01, 0.1, 0.5],
33	'num_iterations': [1, 5, 10]},
34	{'cutoff_freq': (1e-4, 1-1e-2),
35	'num_iterations': (1, 1000)}),
36	finitediff: ({'num_iterations': [5, 10, 30, 50],
37	'order': {2, 4}}, # order is categorical here, because it can't be 3
38	{'num_iterations': (1, 1000)}),
39	first_order: ({'num_iterations': [5, 10, 30, 50]}, # Separated because optimizing over this one is rare due to shifted answer
40	{'num_iterations': (1, 1000)}),
41	polydiff: ({'step_size': [1, 2, 5],
42	'kernel': {'friedrichs', 'gaussian'}, # categorical
43	'degree': [2, 3, 5, 7],
44	'window_size': [11, 31, 51, 91, 131]},
45	{'step_size': (1, 100),
46	'degree': (1, 8),
47	'window_size': (10, 1000)}),
48	savgoldiff: ({'degree': [2, 3, 5, 7, 10],
49	'window_size': [3, 10, 30, 50, 90, 130, 200, 300],
50	'smoothing_win': [3, 10, 30, 50, 90, 130, 200, 300]},
51	{'degree': (1, 12),
52	'window_size': (3, 1000),
53	'smoothing_win': (3, 1000)}),
54	splinediff: ({'degree': {3, 4, 5}, # categorical, because degree is whole number, and there aren't many choices
55	's': [0.5, 0.9, 0.95, 1, 10, 100],
56	'num_iterations': [1, 5, 10]},
57	{'s': (1e-2, 1e6),
58	'num_iterations': (1, 10)}),
59	spectraldiff: ({'even_extension': {True, False}, # give categorical params in a set
60	'pad_to_zero_dxdt': {True, False},
61	'high_freq_cutoff': [1e-3, 5e-2, 1e-2, 5e-2, 1e-1]}, # give numerical params in a list to scipy.optimize over them
62	{'high_freq_cutoff': (1e-5, 1-1e-5)}),
63	rbfdiff: ({'sigma': [1e-3, 1e-2, 1e-1, 1],
64	'lmbd': [1e-3, 1e-2, 1e-1]},
65	{'sigma': (1e-3, 1e3),
66	'lmbd': (1e-4, 0.5)}),
67	tvrdiff: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
68	'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
69	{'gamma': (1e-4, 1e7)}),
70	velocity: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]}, # Deprecated method
71	{'gamma': (1e-4, 1e7)}),
72	iterative_velocity: ({'scale': 'small', # Rare to optimize this one, because it's longer-running than convex version
73	'num_iterations': [1, 5, 10],
74	'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]},
75	{'num_iterations': (1, 100), # gets expensive with more iterations
76	'gamma': (1e-4, 1e7)}),
77	smooth_acceleration: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
78	'window_size': [3, 10, 30, 50, 90, 130]},
79	{'gamma': (1e-4, 1e7),
80	'window_size': (3, 1000)}),
81	rtsdiff: ({'forwardbackward': {True, False},
82	'order': {1, 2, 3}, # for this few options, the optimization works better if this is categorical
83	'qr_ratio': [10**k for k in range(-9, 10, 2)] + [1e12, 1e16]},
84	{'qr_ratio': [1e-10, 1e20]}), # qr_ratio is usually >>1
85	constant_velocity: ({'q': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8], # Deprecated method
86	'r': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8],
87	'forwardbackward': {True, False}},
88	{'q': (1e-10, 1e10),
89	'r': (1e-10, 1e10)}),
90	robustdiff: ({'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
91	'q': [1e-1, 1e1, 1e4, 1e8, 1e12],
92	'r': [1e-1, 1e1, 1e4, 1e8, 1e12],
93	'proc_huberM': {0, 1, 1.345, 2, 6}, # 0 is l1 norm, 1.345 is Huber 95% "efficiency", 2 assumes about 5% outliers,
94	'meas_huberM': {0, 1, 1.345, 2, 6}}, # and 6 assumes basically no outliers -> l2 norm. Try (1 - norm.cdf(M))*2 to see outlier portion
95	{'q': (1e-1, 1e18),
96	'r': (1e-5, 1e18),
97	'huberM': (0, 5)}), # really only want to use l2 norm when nearby
98	lineardiff: ({'kernel': 'gaussian',
99	'order': 3,
100	'gamma': [1e-1, 1, 10, 100],
101	'window_size': [10, 30, 50, 90, 130]},
102	{'order': (1, 5),
103	'gamma': (1e-3, 1000),
104	'window_size': (15, 1000)})
105	} # Methods with nonunique parameter sets are aliased in the dictionary below
106	for method in [second_order, fourth_order]: # Deprecated, redundant methods	1✔
107	method_params_and_bounds[method] = method_params_and_bounds[first_order]	1✔
108	for method in [mediandiff, gaussiandiff, friedrichsdiff]: # Deprecated methods	1✔
109	method_params_and_bounds[method] = method_params_and_bounds[meandiff]	1✔
110	for method in [acceleration, jerk]: # Deprecated, redundant methods	1✔
111	method_params_and_bounds[method] = method_params_and_bounds[velocity]	1✔
112	method_params_and_bounds[jerk_sliding] = method_params_and_bounds[smooth_acceleration]	1✔
113	for method in [constant_acceleration, constant_jerk]: # Deprecated, redundant methods	1✔
114	method_params_and_bounds[method] = method_params_and_bounds[constant_velocity]	1✔
115
116
117	# This function has to be at the top level for multiprocessing but is only used by optimize.
118	def _objective_function(point, func, x, dt, singleton_params, categorical_params, search_space_types,	1✔
119	dxdt_truth, metric, tvgamma, padding, cache):
120	"""Function minimized by scipy.optimize.minimize, needs to have the form: (point, *args) -> float
121	This is mildly complicated, because "point" controls the settings of a differentiation function, but
122	the method may have numerical and non-numerical parameters, and all such parameters are now passed by
123	keyword arguments. So the encoded `point` has to be decoded to dict.
124
125	:param np.array point: a numerical vector scipy chooses to try in the objective function
126	:param dict singleton_params: maps parameter names to singleton values
127	:param dict categorical_params: maps parameter names to values
128	:param dict search_space_types: maps parameter names to types, for turning float vector point into dict point_params
129	:param multiprocessing.manager.dict cache: available across processes to save results and work
130	Other params documented in `optimize`
131
132	:return: float, cost of this objective at the point
133	"""
134	key = sha1((''.join(f"{v:.3e}" for v in point) + # This hash is stable across processes. Takes bytes	×
135	''.join(str(v) for k,v in sorted(categorical_params.items()))).encode()).digest()
136	if key in cache: return cache[key] # short circuit if this hyperparam combo has already been queried, ~10% savings per #160	×
137
138	point_params = {k:(v if search_space_types[k] == float else int(np.round(v)))	×
139	for k,v in zip(search_space_types, point)} # point -> dict
140	# add back in the singletons we're not searching over
141	try: x_hat, dxdt_hat = func(x, dt, point_params, singleton_params, **categorical_params) # estimate x and dxdt	×
142	except np.linalg.LinAlgError: cache[key] = 1e10; return 1e10 # some methods can fail numerically	×
143
144	# evaluate estimate
145	if dxdt_truth is not None: # then minimize \|\|dxdt_hat - dxdt_truth\|\|^2	×
146	if metric == 'rmse':	×
147	rms_rec_x, rms_x, rms_dxdt = evaluate.rmse(x, dt, x_hat, dxdt_hat, dxdt_truth=dxdt_truth, padding=padding)	×
148	cache[key] = rms_dxdt; return rms_dxdt	×
149	elif metric == 'error_correlation':	×
150	ec = evaluate.error_correlation(dxdt_hat, dxdt_truth, padding=padding)	×
151	cache[key] = ec; return ec	×
152	else: # then minimize [ \|\| integral(dxdt_hat) - x\|\|^2 + gamma*TV(dxdt_hat) ]
153	rms_rec_x, rms_x, rms_dxdt = evaluate.rmse(x, dt, x_hat, dxdt_hat, dxdt_truth=None, padding=padding)	×
154	cost = rms_rec_x + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding)	×
155	cache[key] = cost; return cost	×
156
157
158	def optimize(func, x, dt, dxdt_truth=None, tvgamma=1e-2, search_space_updates={}, metric='rmse',	1✔
159	padding=0, opt_method='Nelder-Mead', maxiter=10, parallel=True):
160	"""Find the optimal parameters for a given differentiation method.
161
162	:param function func: differentiation method to optimize parameters for, e.g. linear_model.savgoldiff
163	:param np.array[float] x: data to differentiate
164	:param float dt: step size
165	:param np.array[float] dxdt_truth: actual time series of the derivative of x, if known
166	:param float tvgamma: Only used if :code:`dxdt_truth` is given. Regularization value used to select for parameters
167	that yield a smooth derivative. Larger value results in a smoother derivative.
168	:param dict search_space_updates: At the top of :code:`_optimize.py`, each method has a search space of parameters
169	settings structured as :code:`{param1:[values], param2:[values], param3:value, ...}`. The Cartesian
170	product of values are used as initial starting points in optimization. If left None, the default search
171	space is used.
172	:param str metric: either :code:`'rmse'` or :code:`'error_correlation'`, only applies if :code:`dxdt_truth`
173	is not None, see _objective_function
174	:param int padding: number of time steps to ignore at the beginning and end of the time series in the
175	optimization, or :code:`'auto'` to ignore 2.5% at each end. Larger value causes the
176	optimization to emphasize the accuracy of dxdt in the middle of the time series
177	:param str opt_method: Optimization technique used by :code:`scipy.minimize`, the workhorse
178	:param int maxiter: passed down to :code:`scipy.minimize`, maximum iterations
179	:param bool parallel: whether to use multiple processes to optimize, typically faster for single optimizations, but
180	for experiments it is a better use of resources to pool at that higher level
181
182	:return: tuple[dict, float] of\n
183	- opt_params -- best parameter settings for the differentation method
184	- opt_value -- lowest value found for objective function
185	"""
186	if metric not in ['rmse','error_correlation']:	1✔
187	raise ValueError('`metric` should either be `rmse` or `error_correlation`.')	×
188	if metric == 'error_correlation' and dxdt_truth is None:	1✔
189	raise ValueError('`metric` can only be `error_correlation` if `dxdt_truth` is given.')	×
190
191	search_space, bounds = method_params_and_bounds[func]	1✔
192	search_space.update(search_space_updates) # for things not given, use defaults	1✔
193
194	# No need to optimize over singletons, just pass them through
195	singleton_params = {k:v for k,v in search_space.items() if not isinstance(v, (list, set))}	1✔
196
197	# To handle categoricals, find their combination, and then pass each set individually
198	categorical_params = {k for k,v in search_space.items() if isinstance(v, set)}	1✔
199	categorical_combos = [dict(zip(categorical_params, combo)) for combo in	1✔
200	product(*[search_space[k] for k in categorical_params])] # ends up [{}] if there are no categorical params
201
202	# The Nelder-Mead's search space is the dimensions where multiple numerical options are given in a list
203	search_space_types = {k:type(v[0]) for k,v in search_space.items() if isinstance(v, list)} # map param name -> type, for converting to and from point	1✔
204	if any(v not in [float, int] for v in search_space_types.values()):	1✔
205	raise ValueError("To optimize over categorical strings or bools, put them in a tuple, not a list.")	×
206	# Cast ints to floats, and we're good to go
207	starting_points = list(product(*[np.array(search_space[k]).astype(float) for k in search_space_types]))	1✔
208	# The numerical space should have bounds
209	bounds = [bounds[k] if k in bounds else # pass these to minimize(). It should respect them.	1✔
210	None for k,v in search_space_types.items()] # None means no bound on a dimension
211
212	results = []	1✔
213	filterwarnings("ignore", '', UserWarning) # An extra filtering call, because some worker work can actually be done in the main process	1✔
214	if parallel:	1✔
215	with Manager() as manager:	1✔
216	cache = manager.dict() # cache answers to avoid expensive repeat queries	1✔
217	with Pool(initializer=filterwarnings, initargs=["ignore", '', UserWarning]) as pool: # The heavy lifting	1✔
218	for categorical_combo in categorical_combos:	1✔
219	# wrap the objective and scipy.optimize.minimize to pass kwargs and a host of other things that remain the same
220	_obj_fun = partial(_objective_function, func=func, x=x, dt=dt, singleton_params=singleton_params,	1✔
221	categorical_params=categorical_combo, search_space_types=search_space_types, dxdt_truth=dxdt_truth,
222	metric=metric, tvgamma=tvgamma, padding=padding, cache=cache)
223	_minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options={'maxiter':maxiter})	1✔
224	results += pool.map(_minimize, starting_points) # returns a bunch of OptimizeResult objects	1✔
225	else: # For experiments, where I want to parallelize optimization calls and am not allowed to have each spawn further processes
226	cache = dict()	×
227	for categorical_combo in categorical_combos:	×
228	_obj_fun = partial(_objective_function, func=func, x=x, dt=dt, singleton_params=singleton_params,	×
229	categorical_params=categorical_combo, search_space_types=search_space_types, dxdt_truth=dxdt_truth,
230	metric=metric, tvgamma=tvgamma, padding=padding, cache=cache)
231	_minimize = partial(scipy.optimize.minimize, _obj_fun, method=opt_method, bounds=bounds, options={'maxiter':maxiter})	×
232	results = [_minimize(p) for p in starting_points]	×
233
234	opt_idx = np.nanargmin([r.fun for r in results])	1✔
235	opt_point = results[opt_idx].x	1✔
236	# results are going to be floats, but that may not be allowed, so convert back to a dict
237	opt_params = {k:(v if search_space_types[k] == float else	1✔
238	int(np.round(v)) if search_space_types[k] == int else
239	v > 0.5) for k,v in zip(search_space_types, opt_point)}
240	opt_params.update(singleton_params)	1✔
241	opt_params.update(categorical_combos[opt_idx//len(starting_points)]) # there are \|starting_points\| results for each combo	1✔
242
243	return opt_params, results[opt_idx].fun	1✔
244
245
246	def suggest_method(x, dt, dxdt_truth=None, cutoff_frequency=None):	1✔
247	"""This is meant as an easy-to-use, automatic way for users with some time on their hands to determine
248	a good method and settings for their data. It calls the optimizer over (almost) all methods in the repo
249	using default search spaces defined at the top of the :code:`pynumdiff/optimize/_optimize.py` file.
250	This routine will take a few minutes to run, especially due to `robustdiff`.
251
252	Excluded:
253	- ``first_order``, because iterating causes drift
254	- ``lineardiff``, ``iterative_velocity``, and ``jerk_sliding``, because they either take too long,
255	can be fragile, or tend not to do best
256	- all ``cvxpy``-based methods if it is not installed
257	- ``velocity`` because it tends to not be best but dominates the optimization process by directly
258	optimizing the second term of the metric :math:`L = \\text{RMSE} \\Big( \\text{trapz}(\\mathbf{
259	\\hat{\\dot{x}}}(\\Phi)) + \\mu, \\mathbf{y} \\Big) + \\gamma \\Big({TV}\\big(\\mathbf{\\hat{
260	\\dot{x}}}(\\Phi)\\big)\\Big)`
261
262	:param np.array[float] x: data to differentiate
263	:param float dt: step size, because most methods are not designed to work with variable step sizes
264	:param np.array[float] dxdt_truth: if known, you can pass true derivative values; otherwise you must use
265	:code: `cutoff_frequency`
266	:param float cutoff_frequency: in Hz, the highest dominant frequency of interest in the signal,
267	used to find parameter :math:`\\gamma` for regularization of the optimization process
268	in the absence of ground truth. See https://ieeexplore.ieee.org/document/9241009.
269	Estimate by (a) counting real number of peaks per second in the data, (b) looking at
270	power spectrum and choosing a cutoff, or (c) making an educated guess.
271
272	:return: tuple[callable, dict] of\n
273	- method -- a reference to the function handle of the differentiation method that worked best
274	- opt_params -- optimal parameter settings for the differentation method
275	"""
276	tvgamma = None	×
277	if dxdt_truth is None: # parameter checking	×
278	if cutoff_frequency is None:	×
279	raise ValueError('Either dxdt_truth or cutoff_frequency must be provided.')	×
280	tvgamma = np.exp(-1.6np.log(cutoff_frequency) -0.71np.log(dt) - 5.1) # See https://ieeexplore.ieee.org/document/9241009	×
281
282	methods = [meandiff, mediandiff, gaussiandiff, friedrichsdiff, butterdiff,	×
283	polydiff, savgoldiff, splinediff, spectraldiff, rbfdiff, finitediff, rtsdiff]
284	try: # optionally skip some methods	×
285	import cvxpy	×
286	methods += [tvrdiff, smooth_acceleration, robustdiff]	×
287	except ImportError:	×
288	warn("CVXPY not installed, skipping tvrdiff, smooth_acceleration, and robustdiff")	×
289
290	best_value = float('inf') # core loop	×
291	for func in tqdm(methods):	×
292	p, v = optimize(func, x, dt, dxdt_truth=dxdt_truth, tvgamma=tvgamma, search_space_updates=(	×
293	{'order':{2,3}} if func in [tvrdiff, robustdiff] else {})) # convex-based with order 1 hack the cost function
294	if v < best_value:	×
295	method = func	×
296	best_value = v	×
297	opt_params = p	×
298
299	return method, opt_params	×

florisvb / PyNumDiff / 18988228033

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous