From 888d6580ec63074ca7dd12c8e3650ad4ef4d3ee2 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Tue, 30 Jan 2024 12:39:36 -0500 Subject: [PATCH] tune accuracy of integrate_discrete; improve plotting --- aerosandbox/numpy/integrate_discrete.py | 50 ++++--------------- ..._curvature_display_spectral_convergence.py | 43 +++++++--------- 2 files changed, 28 insertions(+), 65 deletions(-) diff --git a/aerosandbox/numpy/integrate_discrete.py b/aerosandbox/numpy/integrate_discrete.py index 2b030ce7..2ada1c6b 100644 --- a/aerosandbox/numpy/integrate_discrete.py +++ b/aerosandbox/numpy/integrate_discrete.py @@ -13,8 +13,8 @@ def integrate_discrete_intervals( method_endpoints: str = "lower_order", ): """ - Given a set of points (x_i, f_i) from a function, computes the integral of that function over each set of - adjacent points ("intervals"). + Given a set of sampled points (x_i, f_i) from a function, computes the integral of that function over each set of + adjacent points ("intervals"). Does this via a reconstruction approach, with several methods available. In general, N points will yield N-1 integrals (one for each "interval" between points). @@ -264,13 +264,13 @@ def integrate_discrete_squared_curvature( method: str = "hybrid_simpson_cubic", ): """ - Given a set of points (x_i, f_i) from a function f(x), computes the following quantity: + Given a set of sampled points (x_i, f_i) from a function f(x), computes the following quantity: int_{x[0]}^{x[-1]} (f''(x))^2 dx This is useful for regularization of smooth curves (i.e., encouraging smooth functions as optimization results). - Performs this through one of several methods, specified by `method`: + Performs this through one of several reconstruction-based methods, specified by `method`: * "cubic": On each interval, reconstructs a piecewise cubic polynomial. This cubic is the unique polynomial that passes through the two points at the endpoints of the interval, plus the next point beyond each endpoint @@ -290,7 +290,7 @@ def integrate_discrete_squared_curvature( These two quadratics are then analytically differentiated twice, squared, and integrated over the interval. This requires much less calculation, since the quadratics have uniform curvature over the - interval, causing a lot of things to simplify. The result is then the arithmetic mean of the results of this + interval, causing a lot of things to simplify. The result is then computed by combining the results of this process for the two quadratic reconstructions. This is similar to a Simpson's rule integration, balanced between the two sides of the interval. In @@ -568,7 +568,9 @@ def integrate_discrete_squared_curvature( a = res_backward_simpson[slice(None, -1)] b = res_forward_simpson[slice(1, None)] - middle_intervals = (a + b) / 2 + + # middle_intervals = (a + b) / 2 + middle_intervals = ((a ** 2 + b ** 2) / 2 + 1e-100) ** 0.5 # This is more accurate across all frequencies last_interval = res_backward_simpson[slice(-1, None)] @@ -580,36 +582,6 @@ def integrate_discrete_squared_curvature( return res - elif method in ["subgradient"]: # TODO: Is this a duplicate of simpson? - x1 = x[:-3] - x2 = x[1:-2] - x3 = x[2:-1] - x4 = x[3:] - - f1 = f[:-3] - f2 = f[1:-2] - f3 = f[2:-1] - f4 = f[3:] - - h = x3 - x2 - hm = x2 - x1 - hp = x4 - x3 - - dfm = f2 - f1 - df = f3 - f2 - dfp = f4 - f3 - - slope = df / h - slopep = dfp / hp - slopem = dfm / hm - - ddfp = slopep - slope - ddfm = slope - slopem - return ( - ddfp ** 2 / h + - ddfm ** 2 / h - ) / 2 - elif method in ["hybrid_simpson_cubic"]: from aerosandbox.numpy.calculus import gradient dfdx = gradient( @@ -619,13 +591,10 @@ def integrate_discrete_squared_curvature( ) h = x[1:] - x[:-1] - f1 = f[:-1] - f2 = f[1:] + df = f[1:] - f[:-1] dfdx1 = dfdx[:-1] dfdx2 = dfdx[1:] - df = f2 - f1 - res = ( 4 * (dfdx1 ** 2 + dfdx1 * dfdx2 + dfdx2 ** 2) / h + 12 * df / h ** 2 * (df / h - dfdx1 - dfdx2) @@ -712,7 +681,6 @@ def f(x): f=f_vals, x=x_vals, # method="simpson" - # method="subgradient", method="hybrid_simpson_cubic", ) integral = np.sum(approx) diff --git a/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_display_spectral_convergence.py b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_display_spectral_convergence.py index 4cc08659..6d300415 100644 --- a/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_display_spectral_convergence.py +++ b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_display_spectral_convergence.py @@ -5,14 +5,14 @@ n_samples = 10000 -x = s.symbols("x") -k = s.symbols("k", integer=True, positive=True) +x = s.symbols("x", real=True) +k = s.symbols("k", positive=True, real=True) f = s.cos(k * x * 2 * s.pi) / k ** 2 d2fdx2 = f.diff(x, 2) exact = s.simplify(s.integrate(d2fdx2 ** 2, (x, 0, 1))) -print(float(exact)) +@np.vectorize def get_approx(period=10, method="cubic"): x_vals = np.linspace(0, 1, n_samples).astype(float) f_vals = s.lambdify( @@ -20,12 +20,6 @@ def get_approx(period=10, method="cubic"): f.subs(k, (n_samples - 1) / period), )(x_vals) - # import matplotlib.pyplot as plt - # import aerosandbox.tools.pretty_plots as p - # fig, ax = plt.subplots() - # ax.plot(x_vals, f_vals, ".-") - # p.show_plot() - approx = np.sum(integrate_discrete_squared_curvature( f=f_vals, x=x_vals, @@ -34,30 +28,31 @@ def get_approx(period=10, method="cubic"): return approx -periods = np.geomspace(1, 10000, 201) -# approxes = get_approx(periods) -# rel_errors = np.abs(approxes - float(exact)) / float(exact) + +periods = np.geomspace(2, n_samples, 1001) +exacts = s.lambdify(k, exact)((n_samples - 1) / periods) import matplotlib.pyplot as plt import aerosandbox.tools.pretty_plots as p + fig, ax = plt.subplots(2, 1, figsize=(6, 8)) for method in ["cubic", "simpson", "hybrid_simpson_cubic"]: - approxes = np.vectorize( - lambda period: get_approx(period, method=method) - )(periods) - rel_errors = np.abs(approxes - float(exact)) / float(exact) - ax[0].loglog(periods, rel_errors, ".-", label=method, alpha=0.8) - ax[1].semilogx(periods, approxes, ".-", label=method, alpha=0.8) + approxes = get_approx(periods, method) + ratio = approxes / exacts + rel_errors = np.abs(approxes - exacts) / exacts + ax[0].loglog(periods, rel_errors, ".-", label=method, alpha=0.8, markersize=3) + ax[1].semilogx(periods, ratio, ".-", label=method, alpha=0.8, markersize=3) plt.xlabel("Period [samples]") ax[0].set_ylabel("Relative Error") -ax[1].set_ylabel("Approximation") +ax[1].set_ylabel("Approximation / Exact") plt.sca(ax[1]) -p.hline( - float(exact), +ax[1].set_ylim(bottom=0) +ax[1].plot( + periods, + np.ones_like(periods), + label="Exact", color="k", linestyle="--", - label="Exact", alpha=0.5 ) - -p.show_plot() \ No newline at end of file +p.show_plot()