realized dt was being used as simulation dt in pi_cruise_control and went on a quest to standardize the simulations a little more. All the others had a simdt, which means you can just downsample, but for cc I like how the function comes out with dt=0.01, so I introduced a way to upsample. There was then a consideration about when to add noise, which now happens *after* slicing for all sims, which has the added benefit of not needing to produce so many noise samples and of being sure the number of outliers is always 1%, if I want those. Prior to this change, outliers might be sampled in at a random 1% of 40000k locations, and then the signal would be downsampled after, which might leave behind more or fewer of them.

Author: pavelkomarov

pynumdiff/tests/test_utils.py

@@ -96,22 +96,23 @@ def test_simulations(request):
     """Just sprint through running them all to make sure they go. Optionally plot with flag."""
     if request.config.getoption("--plot"):
         from matplotlib import pyplot
-        fig, axes = pyplot.subplots(2, 3, figsize=(18,7), constrained_layout=True)
-
-    for i,(sim,title) in enumerate(zip([pi_cruise_control, sine, triangle, pop_dyn, linear_autonomous, lorenz_x],
-        ["Cruise Control", "Sum of Sines", "Triangles", "Logistic Growth", "Linear Autonomous", "Lorenz First Dimension"])):
-    
-        y, x, dxdt = sim(duration=4, dt=0.01, noise_type='normal', noise_parameters=[0,0.1])
-        assert len(y) == len(x) == len(dxdt)
-
-        if request.config.getoption("--plot"):
-            t = np.arange(len(y))*0.01
-            ax = axes[i//3, i%3]
-            ax.plot(t, x, 'k--', linewidth=3, label=r"true $x$")
-            ax.plot(t, y, '.', color='blue', zorder=-100, markersize=5, label="noisy data")
-            if i//3 == 0: ax.set_xticklabels([])
-            ax.set_title(title, fontsize=18)
-            if i == 5: ax.legend(loc='lower right', fontsize=12)
+        axes = [pyplot.subplots(2, 3, figsize=(18,7), constrained_layout=True)[1] for i in range(3)]
+
+    for j,dt in enumerate([0.005, 0.01, 0.02]):
+        for i,(sim,title) in enumerate(zip([pi_cruise_control, sine, triangle, pop_dyn, linear_autonomous, lorenz_x],
+            ["Cruise Control", "Sum of Sines", "Triangles", "Logistic Growth", "Linear Autonomous", "Lorenz First Dimension"])):
+
+            y, x, dxdt = sim(duration=4, dt=dt, noise_type='normal', noise_parameters=[0,0.1], outliers=True)
+            assert len(y) == len(x) == len(dxdt) == 4/dt # duration/dt
+
+            if request.config.getoption("--plot"):
+                t = np.arange(len(y))*dt
+                ax = axes[j][i//3, i%3]
+                ax.plot(t, x, 'k--', linewidth=3, label=r"true $x$")
+                ax.plot(t, y, '.', color='blue', zorder=-100, markersize=5, label="noisy data")
+                if i//3 == 0: ax.set_xticklabels([])
+                ax.set_title(title, fontsize=18)
+                if i == 5: ax.legend(loc='lower right', fontsize=12)
 
 
 def test_robust_rme():

pynumdiff/utils/simulate.py

@@ -109,10 +109,10 @@ def triangle(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outlier
 
     pos = np.interp(t, reversal_ts, reversal_vals)
     _, vel = finitediff(pos, dt=simdt)
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
+    idx = np.arange(0, len(t), int(dt/simdt)) # downsample
+    noisy_pos = _add_noise(pos[idx], random_seed, noise_type, noise_parameters, outliers)
 
-    idx = np.arange(0, len(t), int(dt/simdt))
-    return noisy_pos[idx], pos[idx], vel[idx]
+    return noisy_pos, pos[idx], vel[idx]
 
 
 def pop_dyn(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False, random_seed=1,
@@ -182,15 +182,15 @@ def linear_autonomous(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
     xs = np.vstack(xs).T
     pos = xs[0,:]
 
-    smooth_pos, vel = finitediff(pos, simdt)
-    noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
-
+    pos, vel = finitediff(pos, simdt)
     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]
+    noisy_pos = _add_noise(pos[idx], random_seed, noise_type, noise_parameters, outliers)
+
+    return noisy_pos, pos[idx], vel[idx]
 
 
 def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5), outliers=False,
-               random_seed=1, dt=0.01):
+               random_seed=1, dt=0.01, simdt=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
     up and down hills. We assume the car has arbitrary power and can achieve whatever acceleration it wants;
@@ -215,7 +215,7 @@ def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
             - **vel** -- a true derivative information of the time series
     """
     # disturbance
-    t = np.arange(0, duration, dt)
+    t = np.arange(0, duration, simdt)
     slope = 0.01*(np.sin(2*np.pi*t) + 0.3*np.sin(4*2*np.pi*t + 0.5) + 1.2*np.sin(1.7*2*np.pi*t + 0.5))
 
     # parameters
@@ -226,33 +226,33 @@ def pi_cruise_control(duration=4, noise_type='normal', noise_parameters=(0, 0.5)
     vd = 0.5 # desired velocity
 
     # Here state is [pos, vel, accel, cumulative pos error]
-    A = np.array([[1,  dt, (dt**2)/2, 0], # Taylor expand out to accel
-                  [0,   1,    dt,     0],
-                  [0, -fr,     0,    ki/(dt**2)], # (pos error) / dt^2 puts it in units of accel
-                  [0,   0,     0,     1]])
+    A = np.array([[1,  simdt, (simdt**2)/2,             0], # Taylor expand out to accel
+                  [0,      1,        simdt,             0],
+                  [0,    -fr,            0, ki/(simdt**2)], # (pos error) / dt^2 puts it in units of accel
+                  [0,      0,            0,             1]])
 
     # Here inputs are [slope, vel_desired - vel_estimated]
     B = np.array([[0,   0],
                   [0,   0],
-                  [-mg, kp/dt], # (vel error) / dt puts it in units of accel
-                  [0,   dt]])
+                  [-mg, kp/simdt], # (vel error) / dt puts it in units of accel
+                  [0,   simdt]])
 
     # run simulation
     states = [np.array([0, 0, 0, 0])] # x0 is all zeros
     controls = []
-    for i in range(len(slope)):
+    for i in range(len(t)):
         u = np.array([slope[i], vd - states[-1][1]]) # current vel is in 1st position of last state
-        xnew = A @ states[-1] + B @ u
-        states.append(xnew)
+        states.append(A @ states[-1] + B @ u)
         controls.append(u)
 
     states = np.vstack(states).T
     controls = np.vstack(controls).T
 
-    pos = states[0, :]
-    vel = states[1, :]
+    idx = slice(1, len(t)+1, int(dt/simdt)) if dt >= simdt else np.arange(1, len(t)+1, dt/simdt).astype(int) # upsampling allowed here
+    pos = states[0, idx]
+    vel = states[1, idx]
     noisy_pos = _add_noise(pos, random_seed, noise_type, noise_parameters, outliers)
-
+    
     return noisy_pos, pos, vel