optimization now uses sets instead of tuples for categorical params, because that's more applicable notation; changed the way the t vector is calculated in splinediff to be more robust and always be same length as x; added outliers param to simulations

Author: pavelkomarov

pynumdiff/optimize/_optimize.py

@@ -17,18 +17,18 @@
 
 # Map from method -> (search_space, bounds_low_hi)
 method_params_and_bounds = {
-    spectraldiff: ({'even_extension': (True, False), # give boolean or numerical params in a list to scipy.optimize over them
-                   'pad_to_zero_dxdt': (True, False),
-                   'high_freq_cutoff': [1e-3, 5e-2, 1e-2, 5e-2, 1e-1]},
+    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)}),
     polydiff: ({'step_size': [1, 2, 5],
-                'kernel': ('friedrichs', 'gaussian'), # use a tuple of multiple options for categoricals
+                'kernel': {'friedrichs', 'gaussian'}, # categorical
                 'poly_order': [2, 3, 5, 7],
                 'window_size': [11, 31, 51, 91, 131]},
                {'step_size': (1, 100),
                 'poly_order': (1, 8),
                 'window_size': (10, 1000)}),
-    savgoldiff: ({'poly_order': [2, 3, 5, 7, 9, 11, 13],
+    savgoldiff: ({'poly_order': [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]},
                  {'poly_order': (1, 12),
@@ -42,26 +42,26 @@
                   'gamma': (1e-3, 1000),
                   'window_size': (15, 1000)}),
     finite_difference: ({'num_iterations': [5, 10, 30, 50],
-                         'order': (2, 4)}, # order is categorical here, because it can't be 3
+                         '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]},
                   {'num_iterations': (1, 1000)}),
     mediandiff: ({'window_size': [5, 15, 30, 50],
                   'num_iterations': [1, 5, 10]},
                 {'window_size': (1, 1e6),
                  'num_iterations': (1, 100)}),
-    butterdiff: ({'filter_order': tuple(i for i in range(1,11)), # categorical to save us from doing double work by guessing between orders
+    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)}),
-    splinediff: ({'order': (3, 4, 5), # categorical, because order is whole number, and there aren't many choices
+    splinediff: ({'order': {3, 4, 5}, # categorical, because order 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)}),
     tvrdiff: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
-               'order': (1, 2, 3)}, # categorical, because order is whole number, and there aren't many choices
+               'order': {1, 2, 3}}, # categorical, because order is whole number, and there aren't many choices
               {'gamma': (1e-4, 1e7)}),
     velocity: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000]},
                {'gamma': (1e-4, 1e7)}),
@@ -75,7 +75,7 @@
                           {'gamma': (1e-4, 1e7),
                            'window_size': (3, 1000)}),
     rts_const_deriv: ({'forwardbackward': (True, False),
-                       'order': (1, 2, 3), # for this few options, the optimization works better if this is categorical
+                       'order': {1, 2, 3}, # for this few options, the optimization works better if this is categorical
                        'qr_ratio': [1e-16, 1e-12] + [10**k for k in range(-9, 10, 2)] + [1e12, 1e16]},
                       {'qr_ratio': [1e-20, 1e20]}),
     constant_velocity: ({'forwardbackward': (True, False),
@@ -161,10 +161,10 @@ def optimize(func, x, dt, dxdt_truth=None, tvgamma=1e-2, search_space_updates={}
     params.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 params.items() if not isinstance(v, (list, tuple))}
+    singleton_params = {k:v for k,v in params.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 params.items() if isinstance(v, tuple)}
+    categorical_params = {k for k,v in params.items() if isinstance(v, set)}
     categorical_combos = [dict(zip(categorical_params, combo)) for combo in product(*[params[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

pynumdiff/smooth_finite_difference/_smooth_finite_difference.py

@@ -194,7 +194,7 @@ def splinediff(x, dt, params=None, options={}, order=3, s=None, num_iterations=1
         if 'iterate' in options and options['iterate']:
             num_iterations = params[2]
 
-    t = np.arange(0, len(x)*dt, dt)
+    t = np.arange(len(x))*dt
 
     x_hat = x
     for _ in range(num_iterations):

pynumdiff/utils/simulate.py

@@ -10,30 +10,35 @@
 
 
 # pylint: disable-msg=too-many-locals, too-many-arguments, no-member
-def _add_noise(x, random_seed, noise_type, noise_parameters):
+def _add_noise(x, random_seed, noise_type, noise_parameters, outliers=False):
     """Add synthetic noise to data
 
     :param np.array[float] x: data
     :param int random_seed: an integer seed used to initialize the random number generator
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
-    :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param array-like noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
 
     :return: (np.array) -- noisy time series data
     """
     np.random.seed(random_seed)
     noise = np.random.__getattribute__(noise_type)(*noise_parameters, x.shape[0])
+    if outliers:
+        ndxs = np.random.choice(len(noise), len(noise)//100, replace=False) # select 1% of locations
+        noise[ndxs] = (np.random.choice([-1, 1], size=len(ndxs)) + np.random.uniform(-0.5, 0.5, len(ndxs)))*(np.max(x) - np.min(x))
     return x + noise
 
 
-def sine(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_seed=1,
+def sine(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False, random_seed=1,
          dt=0.01, simdt=0.0001, frequencies=(1, 1.7), magnitude=1):
     """Create toy example of time series consisted of sinusoidal modes
 
     :param float duration: governs the length of the series, duration/dt
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -54,19 +59,20 @@ def sine(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_seed
         vel += magnitude/len(frequencies)*np.cos(t*2*np.pi*f)*2*np.pi*f
 
     idx = slice(0, len(t), int(dt/simdt)) # downsample so things are dt apart
-    noisy_pos = _add_noise(pos[idx], random_seed, noise_type, noise_parameters)
+    noisy_pos = _add_noise(pos[idx], random_seed, noise_type, noise_parameters, outliers)
 
     return noisy_pos, pos[idx], vel[idx]
 
 
-def triangle(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_seed=1,
+def triangle(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False, random_seed=1,
              dt=0.01, simdt=0.0001):
     """Create toy example of sharp-edged triangle wave with increasing frequencies
 
     :param float duration: governs the length of the series, duration/dt
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -103,13 +109,13 @@ def triangle(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_
 
     pos = np.interp(t, reversal_ts, reversal_vals)
     _, vel = finite_difference(pos, dt=simdt)
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters)
+    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
 
     idx = np.arange(0, len(t), int(dt/simdt))
     return noisy_pos[idx], pos[idx], vel[idx]
 
 
-def pop_dyn(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_seed=1,
+def pop_dyn(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False, random_seed=1,
             dt=0.01, simdt=0.0001):
     """Create toy example of bounded exponential growth:
     http://www.biologydiscussion.com/population/population-growth/population-growth-curves-ecology/51854
@@ -118,6 +124,7 @@ def pop_dyn(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_s
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -141,19 +148,20 @@ def pop_dyn(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_s
     idx = slice(0, len(t), int(dt/simdt)) # downsample so things are dt apart
     pos = np.array(pos[idx])
     vel = np.array(vel[idx])
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters)
+    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
 
     return noisy_pos, pos, vel
 
 
-def linear_autonomous(duration=4, noise_type='normal', noise_parameters=(0, 0.5),
+def linear_autonomous(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False,
                       random_seed=1, dt=0.01, simdt=0.0001):
     """Create toy example of time series from an autonomous linear system
 
     :param float duration: governs the length of the series, duration/dt
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -175,13 +183,13 @@ def linear_autonomous(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
     pos = xs[0,:]
 
     smooth_pos, vel = finite_difference(pos, simdt)
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters)
+    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
 
     idx = slice(0, len(t), int(dt/simdt)) # downsample so things are dt apart
     return noisy_pos[1:][idx], smooth_pos[1:][idx], vel[1:][idx]
 
 
-def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5),
+def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False,
                random_seed=1, dt=0.01):
     """Create a toy example of linear proportional integral controller with nonlinear control inputs.
     Simulate proportional integral control of a car attempting to maintain constant velocity while going
@@ -193,6 +201,7 @@ def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -242,19 +251,20 @@ def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
 
     pos = states[0, :]
     vel = states[1, :]
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters)
+    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
 
     return noisy_pos, pos, vel
 
 
-def lorenz_x(duration=4, noise_type='normal', noise_parameters=(0, 0.5),
+def lorenz_x(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False,
              random_seed=1, dt=0.01, simdt=0.0001):
     """Create toy example of x component from a lorenz attractor
 
     :param float duration: governs the length of the series, duration/dt
     :param str noise_type: type of noise, compatible with :code:`np.random` functions
                         (eg. 'normal', 'uniform', 'poisson')
     :param noise_parameters: parameters of the noise used in :code:`np.random`, leaving off :code:`size`
+    :param bool outliers: whether to corrupt 1% of the data points with out-of-distribution values
     :param int random_seed: an integer seed used to initialize the random number generator
     :param float dt: the step size
     :param float simdt: simulation step size used to generate the time series, typically smaller than
@@ -269,8 +279,7 @@ def lorenz_x(duration=4, noise_type='normal', noise_parameters=(0, 0.5),
     beta = 8/3
     rho = 45
 
-    def _lorenz_xyz_euler(duration=4, noise_type='normal', noise_parameters=(0, 0.5), random_seed=1,
-               dt=0.01, simdt=0.0001, x0=(5, 1, 3), normalize=True):
+    def _lorenz_xyz_euler(x0=(5, 1, 3), normalize=True):
         """Simulation of Lorenz system with Eular method"""
         t = np.arange(0, duration, simdt)
 
@@ -294,11 +303,10 @@ def _lorenz_xyz_euler(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
         xyz /= f # because xyz is one longer than xyz_dot, leave off last entry
         xyz_dot /= f
 
-        noisy_xyz = np.vstack([_add_noise(xyz[i,:], random_seed+i, noise_type, noise_parameters) for i in range(3)])
+        noisy_xyz = np.vstack([_add_noise(xyz[i,:], random_seed+i, noise_type, noise_parameters, outliers) for i in range(3)])
         return noisy_xyz, xyz, xyz_dot
 
-    def _lorenz_xyz_odeint(duration=4, noise_type='normal', noise_parameters=(0, 0.5),
-                   random_seed=1, dt=0.01, simdt=None, normalize=True):
+    def _lorenz_xyz_odeint(normalize=True):
         """Simulate the Lorenz system with scipy's ODE solver"""
         t = np.linspace(0, duration, int(duration/dt))
 
@@ -320,9 +328,8 @@ def dxyz_dt(xyz, t_): # t_ is the time, unused because the Lorenz system is auto
             xyz /= 20
             xyz_dot /= 20
 
-        noisy_xyz = np.vstack([_add_noise(xyz[i,:], random_seed+i, noise_type, noise_parameters) for i in range(3)])
+        noisy_xyz = np.vstack([_add_noise(xyz[i,:], random_seed+i, noise_type, noise_parameters, outliers) for i in range(3)])
         return noisy_xyz, xyz, xyz_dot
 
-    noisy_pos, pos, vel = _lorenz_xyz_euler(duration, noise_type, noise_parameters,
-                                                    random_seed, dt, simdt)
+    noisy_pos, pos, vel = _lorenz_xyz_euler()
     return noisy_pos[0,:], pos[0,:], vel[0,:]