Examples/nlopt

examples/nlopt/README.md

@@ -0,0 +1,18 @@
+# nlopt Examples
+
+Each sub-directory contains a self-contained example. The order in
+which the examples are to appear is specified in `order.json` (an
+array of directory names in the expected order).
+
+In each example directory you'll find:
+
+* `config.toml` - must conform to the specification outlined here:
+  https://docs.pyscript.net/latest/user-guide/configuration/ This is
+  parsed and ultimately turned into a JSON representation as part of
+  the package's API object.
+* `setup.py` - Python code for contextual and environmental setup,
+  NOT SEEN BY THE END USER, but is run before the `code.py` code is
+  evaluated. Allows us to create useful (IPython) shims, avoid
+  repeating boilerplate and whatnot.
+* `code.py` - the actual code added to the editor which forms the
+  practical example of using the package.

examples/nlopt/constrained_tutorial/code.py

@@ -0,0 +1,71 @@
+# ---------------------------------------------------------------------
+# The classic NLopt tutorial problem: minimize sqrt(x2)
+# subject to x2 >= (a*x1 + b)^3 for two pairs (a, b).
+# ---------------------------------------------------------------------
+import numpy as np
+import nlopt
+import matplotlib.pyplot as plt
+
+
+heading("Constrained minimization with two cubic inequalities")
+note(
+    "Minimize <code>f(x) = sqrt(x2)</code> subject to "
+    "<code>x2 &ge; (2*x1)^3</code> and "
+    "<code>x2 &ge; (-x1 + 1)^3</code>. The known optimum is "
+    "<code>(1/3, 8/27)</code> with value "
+    "<code>sqrt(8/27) &approx; 0.5443</code>."
+)
+
+
+def objective(x, grad):
+    if grad.size > 0:
+        grad[0] = 0.0
+        grad[1] = 0.5 / np.sqrt(x[1])
+    return np.sqrt(x[1])
+
+
+def cubic_constraint(x, grad, a, b):
+    """NLopt expects constraints in the form g(x) <= 0."""
+    if grad.size > 0:
+        grad[0] = 3 * a * (a * x[0] + b) ** 2
+        grad[1] = -1.0
+    return (a * x[0] + b) ** 3 - x[1]
+
+
+opt = nlopt.opt(nlopt.LD_MMA, 2)
+opt.set_min_objective(objective)
+opt.set_lower_bounds([-float("inf"), 0.0])
+# Use lambdas to inject the (a, b) parameters into the constraint.
+opt.add_inequality_constraint(lambda x, g: cubic_constraint(x, g, 2, 0), 1e-8)
+opt.add_inequality_constraint(lambda x, g: cubic_constraint(x, g, -1, 1), 1e-8)
+opt.set_xtol_rel(1e-6)
+
+x_star = opt.optimize([1.234, 5.678])
+f_star = opt.last_optimum_value()
+
+note(
+    f"Optimum at <code>x = ({x_star[0]:.5f}, {x_star[1]:.5f})</code>, "
+    f"<code>f(x) = {f_star:.6f}</code>."
+)
+
+# Visualize the feasible region and the optimum.
+x1 = np.linspace(-0.5, 1.0, 400)
+upper_a = (2 * x1) ** 3
+upper_b = (-x1 + 1) ** 3
+feasible_lower = np.maximum(upper_a, upper_b)
+
+fig, ax = plt.subplots(figsize=(7, 4.5))
+ax.plot(x1, upper_a, "--", color="gray", label="$x_2 = (2 x_1)^3$")
+ax.plot(x1, upper_b, ":", color="gray", label="$x_2 = (-x_1+1)^3$")
+ax.fill_between(x1, feasible_lower, 1.5, color="lightblue",
+                alpha=0.5, label="Feasible region")
+ax.plot(x_star[0], x_star[1], "o", color="crimson",
+        markersize=9, label=f"Optimum ({x_star[0]:.3f}, {x_star[1]:.3f})")
+ax.set_xlim(-0.5, 1.0)
+ax.set_ylim(0, 1.2)
+ax.set_xlabel("$x_1$")
+ax.set_ylabel("$x_2$")
+ax.set_title("Feasible region and constrained optimum")
+ax.legend(loc="upper right", fontsize=9)
+fig.tight_layout()
+display(fig, append=True)

examples/nlopt/constrained_tutorial/config.toml

@@ -0,0 +1 @@
+packages = ["nlopt", "numpy", "matplotlib"]

examples/nlopt/constrained_tutorial/setup.py

@@ -0,0 +1,18 @@
+"""Lightweight setup for later cells. No IPython shim here."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(*args, **kwargs, target=__pyscript_display_target__)
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+

examples/nlopt/global_vs_local/code.py

@@ -0,0 +1,79 @@
+# ---------------------------------------------------------------------
+# When the landscape has many local minima, derivative-free global
+# algorithms shine. Here we use DIRECT-L on the Rastrigin function,
+# a classic test problem with a forest of local minima.
+# ---------------------------------------------------------------------
+
+import numpy as np
+import nlopt
+import matplotlib.pyplot as plt
+
+
+heading("Global optimization on the Rastrigin function")
+note(
+    "Rastrigin in 2D: "
+    "<code>f(x) = 20 + sum(x_i^2 - 10 cos(2&pi; x_i))</code>. "
+    "It has hundreds of local minima inside [-5.12, 5.12]^2, with "
+    "the global minimum at the origin. A local gradient method "
+    "started far from zero will get stuck; a global method finds it."
+)
+
+
+def rastrigin(x, grad):
+    # Derivative-free algorithms call us with grad.size == 0,
+    # so we don't need to fill it in.
+    return 20.0 + np.sum(x * x - 10.0 * np.cos(2.0 * np.pi * x))
+
+
+bounds = 5.12
+
+# A naive local search from a poor starting point.
+local_opt = nlopt.opt(nlopt.LN_NELDERMEAD, 2)
+local_opt.set_min_objective(rastrigin)
+local_opt.set_lower_bounds([-bounds, -bounds])
+local_opt.set_upper_bounds([bounds, bounds])
+local_opt.set_xtol_rel(1e-6)
+x_local = local_opt.optimize([4.0, -3.5])
+f_local = local_opt.last_optimum_value()
+
+# DIRECT-L: a Lipschitzian global search that does not need
+# gradients and respects bound constraints.
+global_opt = nlopt.opt(nlopt.GN_DIRECT_L, 2)
+global_opt.set_min_objective(rastrigin)
+global_opt.set_lower_bounds([-bounds, -bounds])
+global_opt.set_upper_bounds([bounds, bounds])
+# Global methods need a stopping budget; cap evaluations.
+global_opt.set_maxeval(2000)
+x_global = global_opt.optimize([4.0, -3.5])
+f_global = global_opt.last_optimum_value()
+
+note(
+    f"Local Nelder-Mead landed at "
+    f"<code>({x_local[0]:.3f}, {x_local[1]:.3f})</code> "
+    f"with <code>f = {f_local:.4f}</code>."
+)
+note(
+    f"Global DIRECT-L landed at "
+    f"<code>({x_global[0]:.3f}, {x_global[1]:.3f})</code> "
+    f"with <code>f = {f_global:.4f}</code> (true minimum is 0)."
+)
+
+# Plot the landscape with both solutions.
+grid = np.linspace(-bounds, bounds, 300)
+gx, gy = np.meshgrid(grid, grid)
+gz = 20.0 + (gx ** 2 - 10 * np.cos(2 * np.pi * gx)) \
+    + (gy ** 2 - 10 * np.cos(2 * np.pi * gy))
+
+fig, ax = plt.subplots(figsize=(7, 5.5))
+mesh = ax.pcolormesh(gx, gy, gz, shading="auto", cmap="viridis")
+ax.plot(x_local[0], x_local[1], "o", color="white",
+        markeredgecolor="black", markersize=10, label="Local (Nelder-Mead)")
+ax.plot(x_global[0], x_global[1], "*", color="red",
+        markeredgecolor="black", markersize=16, label="Global (DIRECT-L)")
+ax.set_title("Rastrigin landscape: local vs. global solver")
+ax.set_xlabel("$x_1$")
+ax.set_ylabel("$x_2$")
+ax.legend(loc="upper right")
+fig.colorbar(mesh, ax=ax, label="f(x)")
+fig.tight_layout()
+display(fig, append=True)

examples/nlopt/global_vs_local/config.toml

@@ -0,0 +1 @@
+packages = ["nlopt", "numpy", "matplotlib"]

examples/nlopt/global_vs_local/setup.py

@@ -0,0 +1,17 @@
+"""Lightweight setup for later cells. No IPython shim here."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(*args, **kwargs, target=__pyscript_display_target__)
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)

examples/nlopt/minimize_quadratic/code.py

@@ -0,0 +1,58 @@
+"""
+A first taste of NLopt: minimize a smooth, gradient-based objective.
+
+NLopt offers dozens of algorithms for nonlinear optimization,
+both local and global, with or without constraints. The Python
+interface revolves around a single object: `nlopt.opt(algorithm,
+n_dims)`. You set bounds, an objective, stopping criteria, and
+then call `optimize(x0)`.
+
+Docs: https://nlopt.readthedocs.io/en/latest/NLopt_Python_Reference/
+"""
+from IPython.core.display import display, HTML
+
+import numpy as np
+import nlopt
+
+
+heading("Finding the bottom of a tilted bowl")
+note(
+    "We minimize f(x, y) = (x - 3)^2 + 2 * (y + 1)^2. "
+    "The minimum is at (3, -1) with value 0. We use MMA, a "
+    "gradient-based local algorithm, so the objective fills in "
+    "<code>grad</code> in-place when NLopt asks for it."
+)
+
+# Track how many times the objective is called.
+call_count = {"n": 0}
+
+
+def bowl(x, grad):
+    """Objective function with analytic gradient."""
+    call_count["n"] += 1
+    if grad.size > 0:
+        # IMPORTANT: write into grad in place, do not rebind it.
+        grad[0] = 2.0 * (x[0] - 3.0)
+        grad[1] = 4.0 * (x[1] + 1.0)
+    return (x[0] - 3.0) ** 2 + 2.0 * (x[1] + 1.0) ** 2
+
+
+# LD_MMA = Local, Derivative-based, Method of Moving Asymptotes.
+opt = nlopt.opt(nlopt.LD_MMA, 2)
+opt.set_min_objective(bowl)
+opt.set_lower_bounds([-10.0, -10.0])
+opt.set_upper_bounds([10.0, 10.0])
+opt.set_xtol_rel(1e-6)
+
+x_star = opt.optimize([0.0, 0.0])
+f_star = opt.last_optimum_value()
+
+note(
+    f"Found minimum at "
+    f"<code>x = ({x_star[0]:.6f}, {x_star[1]:.6f})</code>, "
+    f"with <code>f(x) = {f_star:.3e}</code>, "
+    f"after <strong>{call_count['n']}</strong> objective evaluations."
+)
+
+# NLopt result codes are small positive integers on success.
+note(f"Result code: {opt.last_optimize_result()} (positive means success).")

examples/nlopt/minimize_quadratic/config.toml

@@ -0,0 +1 @@
+packages = ["nlopt", "numpy"]

examples/nlopt/minimize_quadratic/setup.py

@@ -0,0 +1,36 @@
+"""Shim setup for the first example. Includes the full IPython shim."""
+import sys
+import types
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+ipython = types.ModuleType("IPython")
+core = types.ModuleType("IPython.core")
+core_display = types.ModuleType("IPython.core.display")
+core_display.display = display
+core_display.HTML = HTML
+ipython.core = core
+core.display = core_display
+ipython.get_ipython = lambda: None
+ipython.display = core_display
+sys.modules["IPython"] = ipython
+sys.modules["IPython.core"] = core
+sys.modules["IPython.core.display"] = core_display
+sys.modules["IPython.display"] = core_display
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)

examples/nlopt/order.json

@@ -0,0 +1,5 @@
+[
+    "minimize_quadratic",
+    "constrained_tutorial",
+    "global_vs_local"
+]