Add seed options for stochastic NLopt algorithms

docs/source/algorithms.md

@@ -3375,6 +3375,14 @@ algorithms. Please add the
 addition to optimagic when using an NLOPT algorithm. To install nlopt run
 `conda install nlopt`.
 
+PRAXIS, ESCH, ISRES, CRS2-LM, and randomized DIRECT accept a `seed` option. It must be
+an integer between 0 and `2**32 - 1`, inclusive. When set, it initializes NLopt's random
+number generator before the optimization. When `None`, optimagic leaves NLopt's existing
+random state unchanged. The random state is shared by NLopt optimizers on the same
+thread in builds with thread-local storage, and may be shared across threads in other
+builds. Nested optimizations can therefore change the random stream of an outer
+optimization. The deterministic DIRECT variants do not use the seed.
+
 ```{eval-rst}
 .. dropdown:: nlopt_bobyqa
 
@@ -3451,8 +3459,7 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     fail under its certain rank-preserving transformations (e.g., will lead to
     a non-quadratic shape of the objective function).
 
-    The algorithm is not determenistic and it is not possible to achieve
-    detereminancy via seed setting.
+    The algorithm can take random steps. Use ``seed`` to control those steps.
 
     The algorithm failed on a simple benchmark function with finite parameter bounds.
     Passing arguments `lower_bounds` and `upper_bounds` has been disabled for this
@@ -3476,6 +3483,7 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     - **stopping.maxfun** (int): If the maximum number of function
       evaluation is reached, the optimization stops but we do not count this
       as convergence.
+    - **seed** (int or None): Random seed; see the NLopt notes above.
 
 ```
 
@@ -3807,6 +3815,13 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     - "DIRECT_L_RAND_NOSCAL"
     - "DIRECT_RAND"
 
+    .. note::
+
+       NLopt 2.10.1 selects the deterministic unscaled DIRECT-L implementation for
+       all unscaled variants, including the randomized variant. Its seed therefore
+       has no effect on the search. This has been corrected in NLopt's development
+       branch.
+
     - **convergence.xtol_rel** (float):  Stop when the relative movement
       between parameter vectors is smaller than this.
     - **convergence.xtol_abs** (float): Stop when the absolute movement
@@ -3821,6 +3836,8 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     - **locally_biased** (bool): Whether the "L" version of the algorithm is selected.
     - **random_search** (bool): Whether the randomized version of the algorithm is selected.
     - **unscaled_bounds** (bool): Whether the "NOSCAL" version of the algorithm is selected.
+    - **seed** (int or None): Random seed; see the NLopt notes above.
+
 ```
 
 ```{eval-rst}
@@ -3849,6 +3866,8 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     - **stopping.maxfun** (int): If the maximum number of function
       evaluation is reached, the optimization stops but we do not count this
       as convergence.
+    - **seed** (int or None): Random seed; see the NLopt notes above.
+
 ```
 
 ```{eval-rst}
@@ -3880,6 +3899,8 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
     - **stopping.maxfun** (int): If the maximum number of
       function evaluation is reached, the optimization stops but we do not count
       this as convergence.
+    - **seed** (int or None): Random seed; see the NLopt notes above.
+
 ```
 
 ```{eval-rst}
@@ -3914,6 +3935,8 @@ addition to optimagic when using an NLOPT algorithm. To install nlopt run
       convergence.
     - **population_size** (int): Size of the population. If None, it's set to be
       10 * (number of parameters + 1).
+    - **seed** (int or None): Random seed; see the NLopt notes above.
+
 ```
 
 ## Optimizers from iminuit

src/optimagic/optimizers/nlopt_optimizers.py

@@ -5,8 +5,10 @@
 """
 
 from dataclasses import dataclass
+from typing import Annotated
 
 import numpy as np
+from annotated_types import Ge, Le
 from numpy.typing import NDArray
 
 from optimagic import mark
@@ -36,6 +38,9 @@
     import nlopt
 
 
+_NloptSeed = Annotated[int, Ge(0), Le(2**32 - 1)]
+
+
 @mark.minimizer(
     name="nlopt_bobyqa",
     solver_type=AggregationLevel.SCALAR,
@@ -138,6 +143,7 @@ class NloptPRAXIS(Algorithm):
     convergence_ftol_rel: NonNegativeFloat = CONVERGENCE_FTOL_REL
     convergence_ftol_abs: NonNegativeFloat = CONVERGENCE_FTOL_ABS
     stopping_maxfun: PositiveInt = STOPPING_MAXFUN
+    seed: _NloptSeed | None = None
 
     def _solve_internal_problem(
         self, problem: InternalOptimizationProblem, x0: NDArray[np.float64]
@@ -152,6 +158,7 @@ def _solve_internal_problem(
             convergence_ftol_abs=self.convergence_ftol_abs,
             stopping_max_eval=self.stopping_maxfun,
             algorithm=nlopt.LN_PRAXIS,
+            seed=self.seed,
         )
 
         return res
@@ -569,6 +576,7 @@ class NloptDirect(Algorithm):
     locally_biased: bool = False
     random_search: bool = False
     unscaled_bounds: bool = False
+    seed: _NloptSeed | None = None
 
     def _solve_internal_problem(
         self, problem: InternalOptimizationProblem, x0: NDArray[np.float64]
@@ -603,6 +611,7 @@ def _solve_internal_problem(
             convergence_ftol_abs=self.convergence_ftol_abs,
             stopping_max_eval=self.stopping_maxfun,
             algorithm=algo,
+            seed=self.seed if self.random_search else None,
         )
         return res
 
@@ -629,6 +638,7 @@ class NloptESCH(Algorithm):
     convergence_ftol_rel: NonNegativeFloat = CONVERGENCE_FTOL_REL
     convergence_ftol_abs: NonNegativeFloat = CONVERGENCE_FTOL_ABS
     stopping_maxfun: PositiveInt = STOPPING_MAXFUN_GLOBAL
+    seed: _NloptSeed | None = None
 
     def _solve_internal_problem(
         self, problem: InternalOptimizationProblem, x0: NDArray[np.float64]
@@ -643,6 +653,7 @@ def _solve_internal_problem(
             convergence_ftol_abs=self.convergence_ftol_abs,
             stopping_max_eval=self.stopping_maxfun,
             algorithm=nlopt.GN_ESCH,
+            seed=self.seed,
         )
 
         return res
@@ -670,6 +681,7 @@ class NloptISRES(Algorithm):
     convergence_ftol_rel: NonNegativeFloat = CONVERGENCE_FTOL_REL
     convergence_ftol_abs: NonNegativeFloat = CONVERGENCE_FTOL_ABS
     stopping_maxfun: PositiveInt = STOPPING_MAXFUN_GLOBAL
+    seed: _NloptSeed | None = None
 
     def _solve_internal_problem(
         self, problem: InternalOptimizationProblem, x0: NDArray[np.float64]
@@ -684,6 +696,7 @@ def _solve_internal_problem(
             convergence_ftol_abs=self.convergence_ftol_abs,
             stopping_max_eval=self.stopping_maxfun,
             algorithm=nlopt.GN_ISRES,
+            seed=self.seed,
         )
 
         return res
@@ -712,6 +725,7 @@ class NloptCRS2LM(Algorithm):
     convergence_ftol_abs: NonNegativeFloat = CONVERGENCE_FTOL_ABS
     stopping_maxfun: PositiveInt = STOPPING_MAXFUN_GLOBAL
     population_size: PositiveInt | None = None
+    seed: _NloptSeed | None = None
 
     def _solve_internal_problem(
         self, problem: InternalOptimizationProblem, x0: NDArray[np.float64]
@@ -731,6 +745,7 @@ def _solve_internal_problem(
             stopping_max_eval=self.stopping_maxfun,
             algorithm=nlopt.GN_CRS2_LM,
             population_size=population_size,
+            seed=self.seed,
         )
 
         return res
@@ -749,6 +764,7 @@ def _minimize_nlopt(
     stopping_max_eval=None,
     population_size=None,
     nonlinear_constraints=None,
+    seed=None,
 ):
     """Run actual nlopt optimization argument, set relevant attributes."""
 
@@ -785,6 +801,8 @@ def func(x, grad):
         for constr in _get_nlopt_constraints(nonlinear_constraints, filter_type="ineq"):
             opt.add_inequality_mconstraint(constr["fun"], constr["tol"])
     opt.set_min_objective(func)
+    if seed is not None:
+        nlopt.srand(seed)
     solution_x = opt.optimize(x0)
     return _process_nlopt_results(opt, solution_x, is_global)
 

tests/optimagic/optimizers/test_nlopt_optimizers.py

@@ -0,0 +1,95 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from optimagic import mark
+from optimagic.config import IS_NLOPT_INSTALLED
+from optimagic.exceptions import InvalidAlgoOptionError
+from optimagic.optimization.optimize import minimize
+from optimagic.optimizers.nlopt_optimizers import NloptESCH
+from optimagic.parameters.bounds import Bounds
+
+if IS_NLOPT_INSTALLED:
+    import nlopt
+
+
+pytestmark = pytest.mark.skipif(not IS_NLOPT_INSTALLED, reason="nlopt not installed")
+
+
+def _run(algorithm, options, *, flat=False):
+    evaluations = []
+
+    @mark.scalar
+    def objective(x):
+        evaluations.append(x.copy())
+        return 1.0 if flat else float(x @ x)
+
+    bounds = (
+        None
+        if algorithm == "nlopt_praxis"
+        else Bounds(lower=np.full(3, -1.0), upper=np.full(3, 1.0))
+    )
+    minimize(
+        fun=objective,
+        params=np.array([0.5, -0.4, 0.3]),
+        algorithm=algorithm,
+        algo_options={"stopping_maxfun": 80, **options},
+        collect_history=False,
+        bounds=bounds,
+    )
+    return np.array(evaluations)
+
+
+@pytest.mark.parametrize(
+    ("algorithm", "options"),
+    [
+        ("nlopt_praxis", {}),
+        ("nlopt_esch", {}),
+        ("nlopt_isres", {}),
+        ("nlopt_crs2_lm", {}),
+        ("nlopt_direct", {"locally_biased": True, "random_search": True}),
+    ],
+)
+def test_seed_controls_evaluation_sequence(algorithm, options):
+    # A flat objective makes PRAXIS take random steps.
+    flat = algorithm == "nlopt_praxis"
+    first = _run(algorithm, {**options, "seed": 0}, flat=flat)
+    repeated = _run(algorithm, {**options, "seed": 0}, flat=flat)
+    different = _run(algorithm, {**options, "seed": 1}, flat=flat)
+
+    assert_allclose(first, repeated, rtol=0, atol=0)
+    assert first.shape != different.shape or not np.array_equal(first, different)
+
+
+def test_none_seed_preserves_external_seed():
+    nlopt.srand(123)
+    first = _run("nlopt_esch", {"seed": None})
+    nlopt.srand(123)
+    second = _run("nlopt_esch", {"seed": None})
+    nlopt.srand(456)
+    third = _run("nlopt_esch", {"seed": None})
+
+    assert_allclose(first, second, rtol=0, atol=0)
+    assert first.shape != third.shape or not np.array_equal(first, third)
+
+
+def test_deterministic_direct_does_not_reset_random_state():
+    nlopt.srand(456)
+    expected = _run("nlopt_esch", {"seed": None})
+    nlopt.srand(456)
+    _run("nlopt_direct", {"seed": 123})
+    actual = _run("nlopt_esch", {"seed": None})
+
+    assert_allclose(actual, expected, rtol=0, atol=0)
+
+
+@pytest.mark.parametrize("seed", [0, 2**32 - 1])
+def test_seed_limits_are_valid(seed):
+    assert NloptESCH(seed=seed).seed == seed
+    _run("nlopt_esch", {"seed": seed})
+
+
+@pytest.mark.parametrize("seed", [-1, 2**32])
+def test_seed_outside_portable_range_is_rejected(seed):
+    with pytest.raises(InvalidAlgoOptionError):
+        NloptESCH(seed=seed)