f(x) = sqrt(x2) subject to "
+ "x2 ≥ (2*x1)^3 and "
+ "x2 ≥ (-x1 + 1)^3. The known optimum is "
+ "(1/3, 8/27) with value "
+ "sqrt(8/27) ≈ 0.5443."
+)
+
+
+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 x = ({x_star[0]:.5f}, {x_star[1]:.5f}), "
+ f"f(x) = {f_star:.6f}."
+)
+
+# 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)
diff --git a/examples/nlopt/constrained_tutorial/config.toml b/examples/nlopt/constrained_tutorial/config.toml
new file mode 100644
index 0000000..6eb7e61
--- /dev/null
+++ b/examples/nlopt/constrained_tutorial/config.toml
@@ -0,0 +1 @@
+packages = ["nlopt", "numpy", "matplotlib"]
diff --git a/examples/nlopt/constrained_tutorial/setup.py b/examples/nlopt/constrained_tutorial/setup.py
new file mode 100644
index 0000000..bd3abbe
--- /dev/null
+++ b/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"{text}
"), append=True) + diff --git a/examples/nlopt/global_vs_local/code.py b/examples/nlopt/global_vs_local/code.py new file mode 100644 index 0000000..85c5041 --- /dev/null +++ b/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: " + "f(x) = 20 + sum(x_i^2 - 10 cos(2π x_i)). "
+ "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"({x_local[0]:.3f}, {x_local[1]:.3f}) "
+ f"with f = {f_local:.4f}."
+)
+note(
+ f"Global DIRECT-L landed at "
+ f"({x_global[0]:.3f}, {x_global[1]:.3f}) "
+ f"with f = {f_global:.4f} (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)
diff --git a/examples/nlopt/global_vs_local/config.toml b/examples/nlopt/global_vs_local/config.toml
new file mode 100644
index 0000000..6eb7e61
--- /dev/null
+++ b/examples/nlopt/global_vs_local/config.toml
@@ -0,0 +1 @@
+packages = ["nlopt", "numpy", "matplotlib"]
diff --git a/examples/nlopt/global_vs_local/setup.py b/examples/nlopt/global_vs_local/setup.py
new file mode 100644
index 0000000..e6ff81d
--- /dev/null
+++ b/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"{text}
"), append=True) diff --git a/examples/nlopt/minimize_quadratic/code.py b/examples/nlopt/minimize_quadratic/code.py new file mode 100644 index 0000000..644ddc7 --- /dev/null +++ b/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 " + "grad 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"x = ({x_star[0]:.6f}, {x_star[1]:.6f}), "
+ f"with f(x) = {f_star:.3e}, "
+ f"after {call_count['n']} objective evaluations."
+)
+
+# NLopt result codes are small positive integers on success.
+note(f"Result code: {opt.last_optimize_result()} (positive means success).")
diff --git a/examples/nlopt/minimize_quadratic/config.toml b/examples/nlopt/minimize_quadratic/config.toml
new file mode 100644
index 0000000..a6da9b5
--- /dev/null
+++ b/examples/nlopt/minimize_quadratic/config.toml
@@ -0,0 +1 @@
+packages = ["nlopt", "numpy"]
diff --git a/examples/nlopt/minimize_quadratic/setup.py b/examples/nlopt/minimize_quadratic/setup.py
new file mode 100644
index 0000000..5986b1d
--- /dev/null
+++ b/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"{text}
"), append=True) diff --git a/examples/nlopt/order.json b/examples/nlopt/order.json new file mode 100644 index 0000000..275f5bd --- /dev/null +++ b/examples/nlopt/order.json @@ -0,0 +1,5 @@ +[ + "minimize_quadratic", + "constrained_tutorial", + "global_vs_local" +]