From 2ec495e5f8086a304abe7a5f71571651705dd427 Mon Sep 17 00:00:00 2001
From: Connor Olding <cloningdonor@gmail.com>
Date: Sat, 20 May 2023 19:03:46 -0700
Subject: [PATCH] add optimizers from cobyqa

---
 thursday/candidates/__init__.py |  1 +
 thursday/candidates/cobyqa.py   | 82 +++++++++++++++++++++++++++++++++
 thursday/parties.py             |  1 +
 3 files changed, 84 insertions(+)
 create mode 100644 thursday/candidates/cobyqa.py

diff --git a/thursday/candidates/__init__.py b/thursday/candidates/__init__.py
index 8fa86d4..4be675b 100644
--- a/thursday/candidates/__init__.py
+++ b/thursday/candidates/__init__.py
@@ -1,3 +1,4 @@
+from .cobyqa import COBYQA_OPTIMIZERS
 from .dlib import dlib_cube
 from .evolopy import make_evolopy
 from .fcmaes import (
diff --git a/thursday/candidates/cobyqa.py b/thursday/candidates/cobyqa.py
new file mode 100644
index 0000000..1c3ac9c
--- /dev/null
+++ b/thursday/candidates/cobyqa.py
@@ -0,0 +1,82 @@
+from cobyqa import minimize
+from cobyqa.optimize import EXIT_RHOEND_SUCCESS, EXIT_MAXFEV_WARNING
+import numpy as np
+import sys
+
+
+def _solve(
+    objective, x0, bounds, maxfun, *, rhobeg=None, rhoend=None, npt=None, tweaks=None
+):
+    prng = np.random.default_rng()
+
+    options = dict()
+    if npt is not None:
+        options["npt"] = npt
+    options["rhobeg"] = 1.0 if rhobeg is None else rhobeg
+    options["rhoend"] = 1e-6 if rhoend is None else rhoend
+    options["maxfev"] = maxfun
+    options["maxiter"] = 10_000
+
+    narrowing = False
+
+    advanced = dict()
+    if tweaks is None:
+        pass
+
+    elif tweaks == "bobby":
+        # use settings similar to pybobyqa. this seems to help a lot.
+        if rhobeg is None:
+            options["rhobeg"] = 0.1
+        if rhoend is None:
+            options["rhoend"] = 1e-8
+
+    else:
+        assert False, tweaks
+
+    nf = 0
+    while nf < maxfun:
+        options["maxfev"] = maxfun - nf
+        res = minimize(
+            objective, x0, xl=bounds[0], xu=bounds[1], options=options, **advanced
+        )
+
+        if not res.success and res.status == EXIT_MAXFEV_WARNING:
+            pass  # assert _objective(check) == budget, "nfev mismatch"
+        elif res.success and res.status == EXIT_RHOEND_SUCCESS:
+            pass
+        else:
+            print(f"{tweaks=}", res.success, res.message, file=sys.stderr)
+
+        x0 = prng.uniform(size=len(x0)) if np.allclose(x0, res.x) else res.x
+        nf += res.nfev
+
+        if narrowing:
+            options["rhobeg"] /= 2
+
+    return res
+
+
+def cobyqa_default_cube(objective, size, budget):
+    x0, bounds = np.full(size, 0.5), (np.full(size, 0.0), np.full(size, 1.0))
+    res = _solve(objective, x0, bounds, budget)
+    return res.fun, res.x, res.nfev
+
+
+def cobyqa_bobby_cube(objective, size, budget):
+    x0, bounds = np.full(size, 0.5), (np.full(size, 0.0), np.full(size, 1.0))
+    res = _solve(objective, x0, bounds, budget, tweaks="bobby")
+    return res.fun, res.x, res.nfev
+
+
+def cobyqa_large_bobby_cube(objective, size, budget):
+    x0, bounds = np.full(size, 0.5), (np.full(size, 0.0), np.full(size, 1.0))
+    npt = (size + 1) * (size + 2) // 2
+    res = _solve(objective, x0, bounds, budget, tweaks="bobby", npt=npt)
+    return res.fun, res.x, res.nfev
+
+
+COBYQA_OPTIMIZERS = [
+    cobyqa_default_cube,
+    cobyqa_bobby_cube,
+    cobyqa_large_bobby_cube,
+]
diff --git a/thursday/parties.py b/thursday/parties.py
index 144e187..8d90df5 100644
--- a/thursday/parties.py
+++ b/thursday/parties.py
@@ -115,6 +115,7 @@ else:
 
 parties = dict(
     baseline=BASELINE_OPTIMIZERS,
+    cobyqa=COBYQA_OPTIMIZERS,
     everything=FUCKING_EVERYTHING,
     evolopy=EVOLOPY_OPTIMIZERS,
     negative=PREVIOUSLY_NEGATIVE,