|
19 | 19 |
|
20 | 20 | # Analytic (function, derivative) pairs on which to test differentiation methods. |
21 | 21 | test_funcs_and_derivs = [ |
22 | | - (r"$x(t)=1$", lambda t: np.ones(t.shape), lambda t: np.zeros(t.shape)), # constant |
23 | | - (r"$x(t)=2t+1$", lambda t: 2*t + 1, lambda t: 2*np.ones(t.shape)), # affine |
24 | | - (r"$x(t)=t^2-t+1$", lambda t: t**2 - t + 1, lambda t: 2*t - 1), # quadratic |
25 | | - (r"$x(t)=\sin(3t)+1/2$", lambda t: np.sin(3*t) + 1/2, lambda t: 3*np.cos(3*t)), # sinuoidal |
26 | | - (r"$x(t)=e^t\sin(5t)$", lambda t: np.exp(t)*np.sin(5*t), # growing sinusoidal |
27 | | - lambda t: np.exp(t)*(5*np.cos(5*t) + np.sin(5*t))), |
28 | | - (r"$x(t)=\frac{\sin(8t)}{(t+0.1)^{3/2}}$", lambda t: np.sin(8*t)/((t + 0.1)**(3/2)), # steep challenger |
29 | | - lambda t: ((0.8 + 8*t)*np.cos(8*t) - 1.5*np.sin(8*t))/(0.1 + t)**(5/2))] |
| 22 | + (0, r"$x(t)=1$", lambda t: np.ones(t.shape), lambda t: np.zeros(t.shape)), # constant |
| 23 | + (1, r"$x(t)=2t+1$", lambda t: 2*t + 1, lambda t: 2*np.ones(t.shape)), # affine |
| 24 | + (2, r"$x(t)=t^2-t+1$", lambda t: t**2 - t + 1, lambda t: 2*t - 1), # quadratic |
| 25 | + (3, r"$x(t)=\sin(3t)+1/2$", lambda t: np.sin(3*t) + 1/2, lambda t: 3*np.cos(3*t)), # sinuoidal |
| 26 | + (4, r"$x(t)=e^t\sin(5t)$", lambda t: np.exp(t)*np.sin(5*t), # growing sinusoidal |
| 27 | + lambda t: np.exp(t)*(5*np.cos(5*t) + np.sin(5*t))), |
| 28 | + (5, r"$x(t)=\frac{\sin(8t)}{(t+0.1)^{3/2}}$", lambda t: np.sin(8*t)/((t + 0.1)**(3/2)), # steep challenger |
| 29 | + lambda t: ((0.8 + 8*t)*np.cos(8*t) - 1.5*np.sin(8*t))/(0.1 + t)**(5/2))] |
30 | 30 |
|
31 | 31 | # Call both ways, with kwargs (new) and with params list with default options dict (legacy), to ensure both work |
32 | 32 | diff_methods_and_params = [ |
|
81 | 81 | [(1, 1), (3, 3), (1, 1), (3, 3)]] |
82 | 82 | } |
83 | 83 |
|
84 | | - |
85 | 84 | @mark.filterwarnings("ignore::DeprecationWarning") # I want to test the old and new functionality intentionally |
86 | 85 | @mark.parametrize("diff_method_and_params", diff_methods_and_params) |
87 | | -def test_diff_method(diff_method_and_params, request): |
| 86 | +@mark.parametrize("test_func_and_deriv", test_funcs_and_derivs) |
| 87 | +def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # request gives access to context |
88 | 88 | diff_method, params = diff_method_and_params # unpack |
| 89 | + i, latex_name, f, df = test_func_and_deriv |
89 | 90 |
|
90 | 91 | # some methods rely on cvxpy, and we'd like to allow use of pynumdiff without convex optimization |
91 | 92 | if diff_method in [lineardiff, velocity]: |
92 | 93 | try: import cvxpy |
93 | 94 | except: warn(f"Cannot import cvxpy, skipping {diff_method} test."); return |
94 | 95 |
|
95 | | - plot = request.config.getoption("--plot") # Get the plot flag from pytest configuration |
96 | | - if plot: fig, axes = pyplot.subplots(len(test_funcs_and_derivs), 2, figsize=(12,7)) |
97 | | - |
98 | | - # loop over the test functions |
99 | | - for i,(latex,f,df) in enumerate(test_funcs_and_derivs): |
100 | | - x = f(t) # sample the function |
101 | | - x_noisy = x + noise # add a little noise |
102 | | - dxdt = df(t) # true values of the derivative at samples |
| 96 | + x = f(t) # sample the function |
| 97 | + x_noisy = x + noise # add a little noise |
| 98 | + dxdt = df(t) # true values of the derivative at samples |
103 | 99 |
|
104 | | - # differentiate without and with noise |
105 | | - x_hat, dxdt_hat = diff_method(x, dt, **params) if isinstance(params, dict) else diff_method(x, dt, params) \ |
106 | | - if isinstance(params, list) else diff_method(x, dt) |
107 | | - x_hat_noisy, dxdt_hat_noisy = diff_method(x_noisy, dt, **params) if isinstance(params, dict) \ |
108 | | - else diff_method(x_noisy, dt, params) if isinstance(params, list) else diff_method(x_noisy, dt) |
| 100 | + # differentiate without and with noise |
| 101 | + x_hat, dxdt_hat = diff_method(x, dt, **params) if isinstance(params, dict) else diff_method(x, dt, params) \ |
| 102 | + if isinstance(params, list) else diff_method(x, dt) |
| 103 | + x_hat_noisy, dxdt_hat_noisy = diff_method(x_noisy, dt, **params) if isinstance(params, dict) \ |
| 104 | + else diff_method(x_noisy, dt, params) if isinstance(params, list) else diff_method(x_noisy, dt) |
| 105 | + |
| 106 | + # check x_hat and x_hat_noisy are close to x and that dxdt_hat and dxdt_hat_noisy are close to dxdt |
| 107 | + #print("]\n[", end="") |
| 108 | + for j,(a,b) in enumerate([(x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy)]): |
| 109 | + l2_error = np.linalg.norm(a - b) |
| 110 | + linf_error = np.max(np.abs(a - b)) |
109 | 111 |
|
110 | | - # check x_hat and x_hat_noisy are close to x and that dxdt_hat and dxdt_hat_noisy are close to dxdt |
111 | | - #print("]\n[", end="") |
112 | | - for j,(a,b) in enumerate([(x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy)]): |
113 | | - l2_error = np.linalg.norm(a - b) |
114 | | - linf_error = np.max(np.abs(a - b)) |
115 | | - |
116 | | - #print(f"({int(np.ceil(np.log10(l2_error))) if l2_error> 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ") |
117 | | - log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j] |
118 | | - assert np.linalg.norm(a - b) < 10**log_l2_bound |
119 | | - assert np.max(np.abs(a - b)) < 10**log_linf_bound |
120 | | - if np.linalg.norm(a - b) < 10**(log_l2_bound - 1) or np.max(np.abs(a - b)) < 10**(log_linf_bound - 1): |
121 | | - print(f"Improvement detected for method {diff_method}") |
122 | | - |
123 | | - if plot: |
124 | | - axes[i, 0].plot(t, f(t), label=r"$x(t)$") |
125 | | - axes[i, 0].plot(t, x, 'C0+') |
126 | | - axes[i, 0].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$") |
127 | | - axes[i, 0].plot(t, dxdt_hat, 'C1+') |
128 | | - axes[i, 0].set_ylabel(latex, rotation=0, labelpad=50) |
129 | | - if i < len(test_funcs_and_derivs)-1: axes[i, 0].set_xticklabels([]) |
130 | | - else: axes[i, 0].set_xlabel('t') |
131 | | - if i == 0: axes[i, 0].set_title('noiseless') |
132 | | - axes[i, 1].plot(t, f(t), label=r"$x(t)$") |
133 | | - axes[i, 1].plot(t, x_noisy, 'C0+') |
134 | | - axes[i, 1].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$") |
135 | | - axes[i, 1].plot(t, dxdt_hat_noisy, 'C1+') |
136 | | - if i < len(test_funcs_and_derivs)-1: axes[i, 1].set_xticklabels([]) |
137 | | - else: axes[i, 1].set_xlabel('t') |
138 | | - axes[i, 1].set_yticklabels([]) |
139 | | - if i == 0: axes[i, 1].set_title('with noise') |
| 112 | + #print(f"({int(np.ceil(np.log10(l2_error))) if l2_error> 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ") |
| 113 | + log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j] |
| 114 | + assert np.linalg.norm(a - b) < 10**log_l2_bound |
| 115 | + assert np.max(np.abs(a - b)) < 10**log_linf_bound |
| 116 | + if 0 < np.linalg.norm(a - b) < 10**(log_l2_bound - 1) or 0 < np.max(np.abs(a - b)) < 10**(log_linf_bound - 1): |
| 117 | + print(f"Improvement detected for method {str(diff_method)}") |
140 | 118 |
|
141 | | - if plot: |
142 | | - axes[-1,-1].legend() |
143 | | - pyplot.tight_layout() |
144 | | - pyplot.show() |
| 119 | + if request.config.getoption("--plot"): # Get the plot flag from pytest configuration |
| 120 | + fig, axes = request.config.plots[diff_method] |
| 121 | + axes[i, 0].plot(t, f(t), label=r"$x(t)$") |
| 122 | + axes[i, 0].plot(t, x, 'C0+') |
| 123 | + axes[i, 0].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$") |
| 124 | + axes[i, 0].plot(t, dxdt_hat, 'C1+') |
| 125 | + axes[i, 0].set_ylabel(latex_name, rotation=0, labelpad=50) |
| 126 | + if i < len(test_funcs_and_derivs)-1: axes[i, 0].set_xticklabels([]) |
| 127 | + else: axes[i, 0].set_xlabel('t') |
| 128 | + if i == 0: axes[i, 0].set_title('noiseless') |
| 129 | + axes[i, 1].plot(t, f(t), label=r"$x(t)$") |
| 130 | + axes[i, 1].plot(t, x_noisy, 'C0+') |
| 131 | + axes[i, 1].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$") |
| 132 | + axes[i, 1].plot(t, dxdt_hat_noisy, 'C1+') |
| 133 | + if i < len(test_funcs_and_derivs)-1: axes[i, 1].set_xticklabels([]) |
| 134 | + else: axes[i, 1].set_xlabel('t') |
| 135 | + axes[i, 1].set_yticklabels([]) |
| 136 | + if i == 0: axes[i, 1].set_title('with noise') |
0 commit comments