diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/README.md b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/README.md index 31396774..4e16419b 100644 --- a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/README.md +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/README.md @@ -16,4 +16,8 @@ The performance of AeroSandbox (with CasADi backend) is compared against existin ![benchmark_nd_rosenbrock](./nd_rosenbrock/benchmark_nd_rosenbrock.png) +## AeroSandbox vs. Disciplined Optimization Methods +In this chart, runtime is used instead of function evaluations, because the GPkit API doesn't easily expose this information from the underlying solver. + +![benchmark_gp_beam](./gp_beam/benchmark_gp_beam.png) \ No newline at end of file diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/aerosandbox_times.csv b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/aerosandbox_times.csv new file mode 100644 index 00000000..8861196b --- /dev/null +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/aerosandbox_times.csv @@ -0,0 +1,78 @@ +6,0.0092133,2 +7,0.0093947,2 +8,0.0108742,2 +9,0.0119871,2 +10,0.0100001,2 +11,0.0099996,2 +12,0.0089999,2 +13,0.01,2 +14,0.011,2 +16,0.0106868,2 +17,0.0110006,2 +19,0.0104631,2 +21,0.0111357,2 +23,0.0102841,2 +25,0.0109994,2 +28,0.0111787,2 +31,0.0109992,2 +34,0.0115629,2 +37,0.0116452,2 +41,0.0115181,2 +45,0.0133307,2 +49,0.0128663,2 +54,0.0117651,2 +60,0.0130029,2 +66,0.0129992,2 +72,0.0127002,2 +79,0.0137136,2 +87,0.0140218,2 +96,0.0139998,2 +106,0.0160009,2 +116,0.0169992,2 +128,0.0170001,2 +141,0.0269826,2 +155,0.0183702,2 +170,0.0183993,2 +187,0.0173063,2 +206,0.0177846,2 +226,0.0175911,2 +249,0.0195237,2 +274,0.0210766,2 +301,0.0249957,2 +331,0.0253988,2 +364,0.0265567,2 +401,0.0253783,2 +441,0.0267146,2 +485,0.0288208,2 +534,0.0324002,2 +587,0.0349842,2 +646,0.0595186,2 +710,0.0605462,2 +781,0.0554812,2 +859,0.0807664,2 +945,0.0593648,2 +1040,0.0568037,2 +1144,0.0677671,2 +1258,0.0693652,2 +1384,0.0706875,2 +1522,0.0707606,2 +1675,0.0824281,2 +1842,0.0858479,2 +2026,0.0885053,2 +2229,0.0971889,2 +2452,0.104217,2 +2697,0.1217041,2 +2967,0.1324828,2 +3263,0.1398847,2 +3590,0.1668577,2 +3949,0.1693555,2 +4344,0.169189,2 +4778,0.1819305,2 +5256,0.2108755,2 +5781,0.2239938,2 +6359,0.251578,2 +6995,0.2634903,2 +7695,0.2908258,2 +8464,0.3135229,2 +9311,0.3742784,2 +10242,0.3932527,2 diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.pdf b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.pdf new file mode 100644 index 00000000..8a4ddd93 Binary files /dev/null and b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.pdf differ diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.png b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.png new file mode 100644 index 00000000..93762b14 Binary files /dev/null and b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/benchmark_gp_beam.png differ diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_cvxopt_times.csv b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_cvxopt_times.csv new file mode 100644 index 00000000..feede223 --- /dev/null +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_cvxopt_times.csv @@ -0,0 +1,43 @@ +6,0.089875,nan +7,0.1037407,nan +8,0.1252925,nan +9,0.1469695,nan +10,0.1640054,nan +11,0.1865629,nan +12,0.2174818,nan +13,0.2407184,nan +14,0.2663685,nan +16,0.3855672,nan +17,0.3457179,nan +19,0.3871485,nan +21,0.4178118,nan +23,0.4459063,nan +25,0.5363831,nan +28,0.6282798,nan +31,0.7110482,nan +34,0.8038039,nan +37,0.9210395,nan +41,1.0729629,nan +45,1.3092143,nan +49,1.4039818,nan +54,1.7042651,nan +60,1.9846408,nan +66,2.2841193,nan +72,2.609569,nan +79,3.0817283,nan +87,3.4560014,nan +96,5.9335542,nan +106,6.8537624,nan +116,8.0498565,nan +128,10.0156526,nan +141,11.4241825,nan +155,13.0482585,nan +170,14.9169516,nan +187,16.8947642,nan +206,21.6764217,nan +226,22.8835737,nan +249,27.7233109,nan +274,34.0943054,nan +301,45.1328843,nan +331,50.1436696,nan +364,56.0820475,nan diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_mosek_times.csv b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_mosek_times.csv new file mode 100644 index 00000000..0e9911d6 --- /dev/null +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/gpkit_mosek_times.csv @@ -0,0 +1,59 @@ +6,0.0248533,nan +7,0.0343567,nan +8,0.0268454,nan +9,0.0316529,nan +10,0.0321701,nan +11,0.0309571,nan +12,0.0367912,nan +13,0.0369759,nan +14,0.0943213,nan +16,0.0397125,nan +17,0.0447057,nan +19,0.0490063,nan +21,0.0618139,nan +23,0.0559357,nan +25,0.0652641,nan +28,0.0682731,nan +31,0.0749005,nan +34,0.0811166,nan +37,0.1059255,nan +41,0.0971243,nan +45,0.2158601,nan +49,0.1333879,nan +54,0.1493238,nan +60,0.1518241,nan +66,0.169105,nan +72,0.2810773,nan +79,0.203531,nan +87,0.305913,nan +96,0.3626245,nan +106,0.4750413,nan +116,0.3205944,nan +128,0.3754863,nan +141,0.6446782,nan +155,0.5003476,nan +170,0.7573683,nan +187,0.7188326,nan +206,0.9415509,nan +226,1.2407348,nan +249,1.4527316,nan +274,1.2807559,nan +301,2.1001261,nan +331,2.2158967,nan +364,3.7489755,nan +401,2.5389554,nan +441,3.0720139,nan +485,4.2228518,nan +534,7.043779,nan +587,7.4332837,nan +646,10.5903846,nan +710,5.2605252,nan +781,18.6644266,nan +859,14.9577836,nan +945,20.682793,nan +1040,23.749001,nan +1144,19.4155647,nan +1258,14.5918238,nan +1384,22.4536076,nan +1522,31.3027442,nan +1675,37.0666398,nan diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/run_times.py b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/run_times.py index fd98891b..01376325 100644 --- a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/run_times.py +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/gp_beam/run_times.py @@ -1,143 +1,230 @@ from aerosandbox.tools.code_benchmarking import time_function import aerosandbox as asb import aerosandbox.numpy as np -from scipy import optimize -import random, itertools +import itertools import matplotlib.patheffects as path_effects +import pytest -# Problem is unimodal for N=2, N=3, and N>=8. Bimodal for 4<=N<=7. Global min is always a vector of ones. - -def get_initial_guess(N): - rng = np.random.default_rng(0) - return rng.uniform(-10, 10, N) - -def objective(x): - return np.mean( - 100 * (x[1:] - x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2 - ) - def solve_aerosandbox(N=10): - opti = asb.Opti() # set up an optimization environment - - x = opti.variable(init_guess=get_initial_guess(N)) - opti.minimize(objective(x)) + import aerosandbox as asb + import aerosandbox.numpy as np - try: - sol = opti.solve(verbose=False, max_iter=100000000) # solve - except RuntimeError: - raise ValueError(f"N={N} failed!") + L = 6 # m, overall beam length + EI = 1.1e4 # N*m^2, bending stiffness + q = 110 * np.ones(N) # N/m, distributed load - if not np.allclose(sol(x), 1, atol=1e-4): - raise ValueError(f"N={N} failed!") + x = np.linspace(0, L, N) # m, node locations - return sol.stats()['n_call_nlp_f'] + opti = asb.Opti() + w = opti.variable(init_guess=np.zeros(N)) # m, displacement -def solve_scipy_bfgs(N=10): - res = optimize.minimize( - fun=objective, - x0=get_initial_guess(N), - method="BFGS", - tol=1e-8, - options=dict( - maxiter=np.Inf, - ) + th = opti.derivative_of( # rad, slope + w, with_respect_to=x, + derivative_init_guess=np.zeros(N), ) - if not np.allclose(res.x, 1, atol=1e-4): - raise ValueError(f"N={N} failed!") - - return res.nfev - - -def solve_scipy_slsqp(N=10): - res = optimize.minimize( - fun=objective, - x0=get_initial_guess(N), - method="SLSQP", - tol=1e-8, - options=dict( - maxiter=1000000000, - ) + M = opti.derivative_of( # N*m, moment + th * EI, with_respect_to=x, + derivative_init_guess=np.zeros(N), ) - if not np.allclose(res.x, 1, atol=1e-4): - raise ValueError(f"N={N} failed!") - - return res.nfev - - -def solve_scipy_nm(N=10): - res = optimize.minimize( - fun=objective, - x0=get_initial_guess(N), - method="Nelder-Mead", - options=dict( - maxiter=np.Inf, - maxfev=np.Inf, - xatol=1e-8, - adaptive=True, - ) + V = opti.derivative_of( # N, shear + M, with_respect_to=x, + derivative_init_guess=np.zeros(N), ) - if not np.allclose(res.x, 1, atol=1e-4): - raise ValueError(f"N={N} failed!") - - return res.nfev + opti.constrain_derivative( + variable=V, with_respect_to=x, + derivative=q, + ) + opti.subject_to([ + w[0] == 0, + th[0] == 0, + M[-1] == 0, + V[-1] == 0, + ]) + + sol = opti.solve(verbose=False) + + print(sol(w[-1])) + assert sol(w[-1]) == pytest.approx(1.62, abs=0.01) + + return sol.stats()['n_call_nlp_f'] # return number of function evaluations + + +eps = 1e-20 # has to be quite large for consistent cvxopt printouts; + + +def solve_gpkit_cvxopt(N=10): + import numpy as np + from gpkit import parse_variables, Model, ureg + from gpkit.small_scripts import mag + + class Beam(Model): + """Discretization of the Euler beam equations for a distributed load. + + Variables + --------- + EI [N*m^2] Bending stiffness + dx [m] Length of an element + L 5 [m] Overall beam length + + Boundary Condition Variables + ---------------------------- + V_tip eps [N] Tip loading + M_tip eps [N*m] Tip moment + th_base eps [-] Base angle + w_base eps [m] Base deflection + + Node Variables of length N + -------------------------- + q 100*np.ones(N) [N/m] Distributed load + V [N] Internal shear + M [N*m] Internal moment + th [-] Slope + w [m] Displacement + + Upper Unbounded + --------------- + w_tip + + """ + + @parse_variables(__doc__, globals()) + def setup(self, N=4): + # minimize tip displacement (the last w) + self.cost = self.w_tip = w[-1] + return { + "definition of dx" : L == (N - 1) * dx, + "boundary_conditions" : [ + V[-1] >= V_tip, + M[-1] >= M_tip, + th[0] >= th_base, + w[0] >= w_base + ], + # below: trapezoidal integration to form a piecewise-linear + # approximation of loading, shear, and so on + # shear and moment increase from tip to base (left > right) + "shear integration" : + V[:-1] >= V[1:] + 0.5 * dx * (q[:-1] + q[1:]), + "moment integration" : + M[:-1] >= M[1:] + 0.5 * dx * (V[:-1] + V[1:]), + # slope and displacement increase from base to tip (right > left) + "theta integration" : + th[1:] >= th[:-1] + 0.5 * dx * (M[1:] + M[:-1]) / EI, + "displacement integration": + w[1:] >= w[:-1] + 0.5 * dx * (th[1:] + th[:-1]) + } + + b = Beam(N=N, substitutions={"L": 6, "EI": 1.1e4, "q": 110 * np.ones(N)}) + sol = b.solve('cvxopt', verbosity=0) + + assert sol("w")[-1].to('m').magnitude == pytest.approx(1.62, abs=0.01) + + return np.nan + + +def solve_gpkit_mosek(N=10): + import numpy as np + from gpkit import parse_variables, Model, ureg + from gpkit.small_scripts import mag + + class Beam(Model): + """Discretization of the Euler beam equations for a distributed load. + + Variables + --------- + EI [N*m^2] Bending stiffness + dx [m] Length of an element + L 5 [m] Overall beam length + + Boundary Condition Variables + ---------------------------- + V_tip eps [N] Tip loading + M_tip eps [N*m] Tip moment + th_base eps [-] Base angle + w_base eps [m] Base deflection + + Node Variables of length N + -------------------------- + q 100*np.ones(N) [N/m] Distributed load + V [N] Internal shear + M [N*m] Internal moment + th [-] Slope + w [m] Displacement + + Upper Unbounded + --------------- + w_tip + + """ + + @parse_variables(__doc__, globals()) + def setup(self, N=4): + # minimize tip displacement (the last w) + self.cost = self.w_tip = w[-1] + return { + "definition of dx" : L == (N - 1) * dx, + "boundary_conditions" : [ + V[-1] >= V_tip, + M[-1] >= M_tip, + th[0] >= th_base, + w[0] >= w_base + ], + # below: trapezoidal integration to form a piecewise-linear + # approximation of loading, shear, and so on + # shear and moment increase from tip to base (left > right) + "shear integration" : + V[:-1] >= V[1:] + 0.5 * dx * (q[:-1] + q[1:]), + "moment integration" : + M[:-1] >= M[1:] + 0.5 * dx * (V[:-1] + V[1:]), + # slope and displacement increase from base to tip (right > left) + "theta integration" : + th[1:] >= th[:-1] + 0.5 * dx * (M[1:] + M[:-1]) / EI, + "displacement integration": + w[1:] >= w[:-1] + 0.5 * dx * (th[1:] + th[:-1]) + } -def solve_scipy_genetic(N=10): - res = optimize.differential_evolution( - func=objective, - bounds=[(-10, 10)] * N, - maxiter=1000000000, - x0=get_initial_guess(N), - ) + b = Beam(N=N, substitutions={"L": 6, "EI": 1.1e4, "q": 110 * np.ones(N)}) + sol = b.solve('mosek_conif', verbosity=0) - if not np.allclose(res.x, 1, atol=1e-4): - raise ValueError(f"N={N} failed!") + assert sol("w")[-1].to('m').magnitude == pytest.approx(1.62, abs=0.01) - return res.nfev + return np.nan if __name__ == '__main__': solvers = { - "AeroSandbox": solve_aerosandbox, - "BFGS" : solve_scipy_bfgs, - "SLSQP": solve_scipy_slsqp, - "Nelder-Mead": solve_scipy_nm, - "Genetic" : solve_scipy_genetic, + "AeroSandbox" : solve_aerosandbox, + "GPkit_cvxopt": solve_gpkit_cvxopt, + "GPkit_mosek" : solve_gpkit_mosek, } - if False: + if False: # If True, runs the benchmark and appends data to respective *.csv files for solver_name, solver in solvers.items(): print(f"Running {solver_name}...") - solver(N=2) + solver(N=5) - N_ideal = 2.0 + N_ideal = 5.0 Ns_attempted = [] while True: N_ideal *= 1.1 - # print(f"Trying N_ideal={N_ideal}...") N = np.round(N_ideal).astype(int) if N in Ns_attempted: continue - # if 4 <= N <= 7: - # continue - print(f"Trying N={N}...") Ns_attempted.append(N) try: t, nfev = time_function( lambda: solver(N=N), - # desired_runtime=0.25, - # runtime_reduction=lambda x: np.percentile(x, 5) ) except ValueError: continue @@ -160,19 +247,14 @@ def solve_scipy_genetic(N=10): fig, ax = plt.subplots(figsize=(5.2, 4)) - # Define a list of distinguishable colors - import copy - colors = p.sns.husl_palette( + fallback_colors = itertools.cycle(p.sns.husl_palette( n_colors=len(solvers) - 1, - h=0, - s=0.25, - l=0.6, - ) - fallback_colors = itertools.cycle(colors) + h=0, s=0.25, l=0.6, + )) name_remaps = { - # "aerosandbox": "AeroSandbox", - "Nelder-Mead": "Nelder\nMead", + "GPkit_cvxopt": "GPkit\n(cvxopt)", + "GPkit_mosek" : "GPkit\n(mosek)", } color_remaps = { @@ -181,23 +263,25 @@ def solve_scipy_genetic(N=10): notables = ["AeroSandbox"] - for i, solver_name in enumerate(solvers.keys()): + for i, solver_name in enumerate(solvers.keys()): # For each solver... + # Reads the data from file df = pd.read_csv(f"{solver_name.lower()}_times.csv", header=None, names=["N", "t", "nfev"]) aggregate_cols = [col for col in df.columns if col != 'N'] df = df.groupby('N', as_index=False)[aggregate_cols].mean() df = df.sort_values('N') + # Determines which columns to plot x = df["N"].values - y = df["nfev"].values + y = df["t"].values - label = solver_name - - if label in color_remaps: - color = color_remaps[label] + # Figures out which color to use + if solver_name in color_remaps: + color = color_remaps[solver_name] else: color = next(fallback_colors) + # Plots the raw data line, = plt.plot( x, y, ".", alpha=0.2, @@ -205,33 +289,30 @@ def solve_scipy_genetic(N=10): ) + # Makes a curve fit and plots that def model(x, p): return ( - p["c"] - + np.exp(p["b1"] * np.log(x) + p["a1"]) - + np.exp(p["b2"] * np.log(x) + p["a2"]) + p["c"] + + np.exp(p["b1"] * np.log(x) + p["a1"]) + + np.exp(p["b2"] * np.log(x) + p["a2"]) ) - # return p["a"] * x ** p["b"] + p["c"] - fit = asb.FittedModel( model=model, x_data=x, y_data=y, parameter_guesses={ - "a1": 0, - "b1": 2, + "a1": 1, + "b1": 1, "a2": 1, "b2": 3, - "c": 0, + "c" : 0, }, parameter_bounds={ - "a1": [0, np.inf], "b1": [0, 10], - "a2": [0, np.inf], "b2": [0, 10], - "c": [0, np.min(y)], + "c" : [0, np.min(y)], }, residual_norm_type="L1", put_residuals_in_logspace=True, @@ -246,31 +327,14 @@ def model(x, p): resample_resolution=10000 ) - if label in name_remaps: - label_to_write = name_remaps[label] + # Writes the label for each plot + if solver_name in name_remaps: + label_to_write = name_remaps[solver_name] else: - label_to_write = label - - if label in notables: - # txt = ax.annotate( - # label, - # xy=(x[-1], fit(x[-1])), - # xytext=(0, -8), - # textcoords="offset points", - # fontsize=10, - # zorder=5, - # alpha=0.9, - # color=color, - # horizontalalignment='right', - # verticalalignment='top', - # path_effects=[ - # path_effects.withStroke(linewidth=2, foreground=ax.get_facecolor(), - # alpha=0.8, - # ), - # ], - # rotation=33 - # ) - txt = ax.annotate( + label_to_write = solver_name + + if solver_name in notables: + ax.annotate( label_to_write, xy=(x[-1], fit(x[-1])), xytext=(-5, -45), @@ -288,7 +352,7 @@ def model(x, p): ], ) else: - txt = ax.annotate( + ax.annotate( label_to_write, xy=(x[-1], fit(x[-1])), xytext=(4, 0), @@ -308,21 +372,19 @@ def model(x, p): plt.xscale("log") plt.yscale("log") - plt.xlim(left=1, right=1e4) - plt.ylim(bottom=10) from aerosandbox.tools.string_formatting import eng_string ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, pos: eng_string(x))) - ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, pos: eng_string(x))) + ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, pos: f"{x:.4g}")) p.show_plot( - "AeroSandbox vs.\nBlack-Box Optimization Methods", - # "\nfor the N-Dimensional Rosenbrock Problem", + "AeroSandbox vs. Disciplined Methods" + "\nfor the GP-Compatible Beam Problem", "\nProblem Size\n(# of Design Variables)", - "Computational\nCost\n\n(# of Function\nEvaluations)", + "Computational\nCost\n\n(Wall-clock\nruntime,\nin seconds)", set_ticks=False, legend=False, dpi=600, - savefig="benchmark_nd_rosenbrock.pdf" + savefig=["benchmark_gp_beam.pdf", "benchmark_gp_beam.png"] ) diff --git a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/nd_rosenbrock/run_times.py b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/nd_rosenbrock/run_times.py index d3930c94..5720c6a6 100644 --- a/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/nd_rosenbrock/run_times.py +++ b/tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/nd_rosenbrock/run_times.py @@ -2,7 +2,7 @@ import aerosandbox as asb import aerosandbox.numpy as np from scipy import optimize -import random, itertools +import itertools import matplotlib.patheffects as path_effects @@ -111,33 +111,27 @@ def solve_scipy_genetic(N=10): "Genetic" : solve_scipy_genetic, } - if False: + if False: # If True, runs the benchmark and appends data to respective *.csv files for solver_name, solver in solvers.items(): print(f"Running {solver_name}...") - solver(N=2) + solver(N=5) - N_ideal = 2.0 + N_ideal = 5.0 Ns_attempted = [] while True: N_ideal *= 1.1 - # print(f"Trying N_ideal={N_ideal}...") N = np.round(N_ideal).astype(int) if N in Ns_attempted: continue - # if 4 <= N <= 7: - # continue - print(f"Trying N={N}...") Ns_attempted.append(N) try: t, nfev = time_function( lambda: solver(N=N), - # desired_runtime=0.25, - # runtime_reduction=lambda x: np.percentile(x, 5) ) except ValueError: continue @@ -160,18 +154,12 @@ def solve_scipy_genetic(N=10): fig, ax = plt.subplots(figsize=(5.2, 4)) - # Define a list of distinguishable colors - import copy - colors = p.sns.husl_palette( + fallback_colors = itertools.cycle(p.sns.husl_palette( n_colors=len(solvers) - 1, - h=0, - s=0.25, - l=0.6, - ) - fallback_colors = itertools.cycle(colors) + h=0, s=0.25, l=0.6, + )) name_remaps = { - # "aerosandbox": "AeroSandbox", "Nelder-Mead": "Nelder\nMead", } @@ -181,23 +169,25 @@ def solve_scipy_genetic(N=10): notables = ["AeroSandbox"] - for i, solver_name in enumerate(solvers.keys()): + for i, solver_name in enumerate(solvers.keys()): # For each solver... + # Reads the data from file df = pd.read_csv(f"{solver_name.lower()}_times.csv", header=None, names=["N", "t", "nfev"]) aggregate_cols = [col for col in df.columns if col != 'N'] df = df.groupby('N', as_index=False)[aggregate_cols].mean() df = df.sort_values('N') + # Determines which columns to plot x = df["N"].values y = df["nfev"].values - label = solver_name - - if label in color_remaps: - color = color_remaps[label] + # Figures out which color to use + if solver_name in color_remaps: + color = color_remaps[solver_name] else: color = next(fallback_colors) + # Plots the raw data line, = plt.plot( x, y, ".", alpha=0.2, @@ -205,6 +195,7 @@ def solve_scipy_genetic(N=10): ) + # Makes a curve fit and plots that def model(x, p): return ( p["c"] @@ -246,31 +237,14 @@ def model(x, p): resample_resolution=10000 ) - if label in name_remaps: - label_to_write = name_remaps[label] + # Writes the label for each plot + if solver_name in name_remaps: + label_to_write = name_remaps[solver_name] else: - label_to_write = label - - if label in notables: - # txt = ax.annotate( - # label, - # xy=(x[-1], fit(x[-1])), - # xytext=(0, -8), - # textcoords="offset points", - # fontsize=10, - # zorder=5, - # alpha=0.9, - # color=color, - # horizontalalignment='right', - # verticalalignment='top', - # path_effects=[ - # path_effects.withStroke(linewidth=2, foreground=ax.get_facecolor(), - # alpha=0.8, - # ), - # ], - # rotation=33 - # ) - txt = ax.annotate( + label_to_write = solver_name + + if solver_name in notables: + ax.annotate( label_to_write, xy=(x[-1], fit(x[-1])), xytext=(-5, -45), @@ -288,7 +262,7 @@ def model(x, p): ], ) else: - txt = ax.annotate( + ax.annotate( label_to_write, xy=(x[-1], fit(x[-1])), xytext=(4, 0),