From fa9d22d4f413cff83aa6bed8e8587a16aa3590f6 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Wed, 17 Jan 2024 15:32:43 -0500 Subject: [PATCH] improve log-factorial calculation for improved speed --- .../time_series_uncertainty_quantification.py | 129 +++++++++--------- 1 file changed, 67 insertions(+), 62 deletions(-) diff --git a/aerosandbox/tools/statistics/time_series_uncertainty_quantification.py b/aerosandbox/tools/statistics/time_series_uncertainty_quantification.py index 6f66d86a6..e5bf93d87 100644 --- a/aerosandbox/tools/statistics/time_series_uncertainty_quantification.py +++ b/aerosandbox/tools/statistics/time_series_uncertainty_quantification.py @@ -5,6 +5,8 @@ import aerosandbox.numpy as np from aerosandbox.tools.pretty_plots.utilities.natural_univariate_spline import NaturalUnivariateSpline as Spline +from scipy import signal +from aerosandbox.tools.code_benchmarking import Timer def estimate_noise_standard_deviation( @@ -29,9 +31,6 @@ def estimate_noise_standard_deviation( The algorithm used in this function is a highly-optimized version of the math described in this repository, part of an upcoming paper: https://github.com/peterdsharpe/aircraft-polar-reconstruction-from-flight-test - The repository is currently private, but will be public at some point; if you would like access to it, - please contact Peter Sharpe at pds@mit.edu. - Args: data: A 1D NumPy array of time-series data. @@ -46,28 +45,34 @@ def estimate_noise_standard_deviation( raise ValueError("Data must have at least 2 points.") if estimator_order is None: - estimator_order = min( - max( - 1, - int(len(data) ** 0.5) - ), - 1000 - ) + estimator_order = int(np.clip(len(data) ** 0.5, 1, 1000)) ##### Noise Variance Reconstruction ##### from scipy.special import gammaln ln_factorial = lambda x: gammaln(x + 1) ### For speed, pre-compute the log-factorial of integers from 1 to estimator_order - ln_f = ln_factorial(np.arange(estimator_order + 1)) + # ln_f = ln_factorial(np.arange(estimator_order + 1)) + ln_f = np.cumsum( + np.log( + np.concatenate([ + [1], + np.arange(1, estimator_order + 1) + ]) + ) + ) ### Create a convolutional kernel to vectorize the summation - coefficients = np.exp( - 2 * ln_f[estimator_order] - ln_f - ln_f[::-1] - 0.5 * ln_factorial(2 * estimator_order) - ) * (-1) ** np.arange(estimator_order + 1) - coefficients -= np.mean(coefficients) # Remove any bias introduced by floating-point error + log_coeffs = ( + 2 * ln_f[estimator_order] - ln_f - ln_f[::-1] - 0.5 * ln_factorial(2 * estimator_order) + ) + indices = np.nonzero(log_coeffs >= np.log(1e-20) + log_coeffs[estimator_order // 2])[0] + coefficients = np.exp(log_coeffs[indices[0]:indices[-1] + 1]) + coefficients[::2] *= -1 # Flip the sign on every other coefficient + coefficients -= np.mean(coefficients) # Remove any bias introduced by floating-point error - sample_stdev = np.convolve(data, coefficients[::-1], 'valid') + # sample_stdev = signal.convolve(data, coefficients[::-1], 'valid') + sample_stdev = signal.oaconvolve(data, coefficients[::-1], 'valid') return np.mean(sample_stdev ** 2) ** 0.5 @@ -244,56 +249,56 @@ def bootstrap_fits( if __name__ == '__main__': - # np.random.seed(1) - # N = 1000 - # f_sample_over_f_signal = 1000 - # - # t = np.arange(N) - # y = np.sin(2 * np.pi / f_sample_over_f_signal * t) + 0.1 * np.random.randn(len(t)) - # - # print(estimate_noise_standard_deviation(y)) + np.random.seed(1) + N = 1000 + f_sample_over_f_signal = 1000 - d = dict(np.load("raw_data.npz")) + t = np.arange(N) + y = np.sin(2 * np.pi / f_sample_over_f_signal * t) + 0.1 * np.random.randn(len(t)) - x = d["airspeed"] - y = d["voltage"] * d["current"] + print(estimate_noise_standard_deviation(y, 1)) - # estimate_noise_standard_deviation(x) + # d = dict(np.load("raw_data.npz")) + # + # x = d["airspeed"] + # y = d["voltage"] * d["current"] # - # x_fit, y_bootstrap_fits = bootstrap_fits( + # # estimate_noise_standard_deviation(x) + # # + # # x_fit, y_bootstrap_fits = bootstrap_fits( + # # x, y, + # # x_stdev=None, + # # y_stdev=None, + # # n_bootstraps=20, + # # spline_degree=5, + # # ) + # import matplotlib.pyplot as plt + # import aerosandbox.tools.pretty_plots as p + # + # fig, ax = plt.subplots(figsize=(7, 4)) + # + # p.plot_with_bootstrapped_uncertainty( # x, y, # x_stdev=None, - # y_stdev=None, - # n_bootstraps=20, - # spline_degree=5, + # y_stdev=estimate_noise_standard_deviation(y[np.argsort(x)]), + # ci=[0.75, 0.95], + # color="coral", + # n_bootstraps=100, + # n_fit_points=200, + # # ci_to_alpha_mapping=lambda ci: 0.4, + # normalize=False, + # spline_degree=3, + # ) + # plt.xlim(x.min(), x.max()) + # plt.ylim(-10, 800) + # p.set_ticks(1, 0.25, 100, 25) + # plt.legend( + # loc="lower right" + # ) + # p.show_plot( + # xlabel="Cruise Airspeed [m/s]", + # ylabel="Electrical Power Required [W]", + # title="Raw Data", + # legend=False, + # dpi=300 # ) - import matplotlib.pyplot as plt - import aerosandbox.tools.pretty_plots as p - - fig, ax = plt.subplots(figsize=(7, 4)) - - p.plot_with_bootstrapped_uncertainty( - x, y, - x_stdev=None, - y_stdev=estimate_noise_standard_deviation(y[np.argsort(x)]), - ci=[0.75, 0.95], - color="coral", - n_bootstraps=100, - n_fit_points=200, - # ci_to_alpha_mapping=lambda ci: 0.4, - normalize=False, - spline_degree=3, - ) - plt.xlim(x.min(), x.max()) - plt.ylim(-10, 800) - p.set_ticks(1, 0.25, 100, 25) - plt.legend( - loc="lower right" - ) - p.show_plot( - xlabel="Cruise Airspeed [m/s]", - ylabel="Electrical Power Required [W]", - title="Raw Data", - legend=False, - dpi=300 - )