Move benchmarks from studies to tutorial; add photo

tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/README.md

@@ -0,0 +1,14 @@
+# Optimization Benchmark Problems
+
+This folder compares optimization performance with AeroSandbox to various other common optimization paradigms.
+
+## AeroSandbox vs. Black-Box Optimization Methods
+
+This chart shows optimization performance on the [N-dimensional Rosenbrock problem](https://en.wikipedia.org/wiki/Rosenbrock_function#Multidimensional_generalizations). Here, $N$ is the number of design variables, which is a convenient knob to dial up or down the difficulty of the problem. The problem is defined as:
+
+* minimize $\sum_{i=1}^{N-1} [ 100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2]$
+
+
+![benchmark_nd_rosenbrock](./nd_rosenbrock/benchmark_nd_rosenbrock.png)
+
+

tutorial/01 - Optimization and Math/01 - Optimization Benchmark Problems/nd_rosenbrock/run_times.py

@@ -0,0 +1,328 @@
+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 matplotlib.patheffects as path_effects
+
+
+# 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))
+
+    try:
+        sol = opti.solve(verbose=False, max_iter=100000000)  # solve
+    except RuntimeError:
+        raise ValueError(f"N={N} failed!")
+
+    if not np.allclose(sol(x), 1, atol=1e-4):
+        raise ValueError(f"N={N} failed!")
+
+    return sol.stats()['n_call_nlp_f']
+
+
+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,
+        )
+    )
+
+    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,
+        )
+    )
+
+    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,
+        )
+    )
+
+    if not np.allclose(res.x, 1, atol=1e-4):
+        raise ValueError(f"N={N} failed!")
+
+    return res.nfev
+
+
+def solve_scipy_genetic(N=10):
+    res = optimize.differential_evolution(
+        func=objective,
+        bounds=[(-10, 10)] * N,
+        maxiter=1000000000,
+        x0=get_initial_guess(N),
+    )
+
+    if not np.allclose(res.x, 1, atol=1e-4):
+        raise ValueError(f"N={N} failed!")
+
+    return res.nfev
+
+
+if __name__ == '__main__':
+
+    solvers = {
+        "AeroSandbox": solve_aerosandbox,
+        "BFGS" : solve_scipy_bfgs,
+        "SLSQP": solve_scipy_slsqp,
+        "Nelder-Mead": solve_scipy_nm,
+        "Genetic"    : solve_scipy_genetic,
+    }
+
+    if False:
+        for solver_name, solver in solvers.items():
+            print(f"Running {solver_name}...")
+            solver(N=2)
+
+            N_ideal = 2.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
+                except KeyboardInterrupt:
+                    break
+                print(f"{solver_name}: N={N}, t={t}, nfev={nfev}")
+                with open(f"{solver_name.lower()}_times.csv", "a") as f:
+                    f.write(f"{N},{t},{nfev}\n")
+
+                if N > 10e3:
+                    break
+                if t > 120:
+                    break
+                if nfev > 1e6:
+                    break
+
+    import matplotlib.pyplot as plt
+    import aerosandbox.tools.pretty_plots as p
+    import pandas as pd
+
+    fig, ax = plt.subplots(figsize=(5.2, 4))
+
+    # Define a list of distinguishable colors
+    import copy
+    colors = p.sns.husl_palette(
+        n_colors=len(solvers) - 1,
+        h=0,
+        s=0.25,
+        l=0.6,
+    )
+    fallback_colors = itertools.cycle(colors)
+
+    name_remaps = {
+        # "aerosandbox": "AeroSandbox",
+        "Nelder-Mead": "Nelder\nMead",
+    }
+
+    color_remaps = {
+        "AeroSandbox": "royalblue",
+    }
+
+    notables = ["AeroSandbox"]
+
+    for i, solver_name in enumerate(solvers.keys()):
+
+        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')
+
+        x = df["N"].values
+        y = df["nfev"].values
+
+        label = solver_name
+
+        if label in color_remaps:
+            color = color_remaps[label]
+        else:
+            color = next(fallback_colors)
+
+        line, = plt.plot(
+            x, y, ".",
+            alpha=0.2,
+            color=color
+        )
+
+
+        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"])
+            )
+
+            # 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,
+                "a2": 1,
+                "b2": 3,
+                "c": 0,
+            },
+            parameter_bounds={
+                "a1": [0, np.inf],
+                "b1": [0, 10],
+                "a2": [0, np.inf],
+                "b2": [0, 10],
+                "c": [0, np.min(y)],
+            },
+            residual_norm_type="L1",
+            put_residuals_in_logspace=True,
+            verbose=False
+        )
+        x_plot = np.geomspace(x.min(), x.max(), 500)
+        p.plot_smooth(
+            x_plot, fit(x_plot), "-",
+            function_of="x",
+            color=color,
+            alpha=0.8,
+            resample_resolution=10000
+        )
+
+        if label in name_remaps:
+            label_to_write = name_remaps[label]
+        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,
+                xy=(x[-1], fit(x[-1])),
+                xytext=(-5, -45),
+                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,
+                                            ),
+                ],
+            )
+        else:
+            txt = ax.annotate(
+                label_to_write,
+                xy=(x[-1], fit(x[-1])),
+                xytext=(4, 0),
+                textcoords="offset points",
+                fontsize=7,
+                zorder=4,
+                alpha=0.7,
+                color=color,
+                horizontalalignment='left',
+                verticalalignment='center',
+                path_effects=[
+                    path_effects.withStroke(linewidth=2, foreground=ax.get_facecolor(),
+                                            alpha=0.3,
+                                            ),
+                ]
+            )
+
+    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)))
+
+    p.show_plot(
+        "AeroSandbox vs.\nBlack-Box Optimization Methods",
+        # "\nfor the N-Dimensional Rosenbrock Problem",
+        "\nProblem Size\n(# of Design Variables)",
+        "Computational\nCost\n\n(# of Function\nEvaluations)",
+        set_ticks=False,
+        legend=False,
+        dpi=600,
+        savefig=["benchmark_nd_rosenbrock.pdf", "benchmark_nd_rosenbrock.png"]
+    )