MAINT: optimize.quadratic_assignment: SPEC-007 transition to rng (scipy#21848)

* MAINT: optimize.quadratic_assignment: SPEC-007 transition to rng 

---------

Co-authored-by: Matt Haberland <[email protected]>

Author: andyfaff

scipy/optimize/_qap.py

@@ -1,5 +1,7 @@
 import numpy as np
 import operator
+import warnings
+import numbers
 from . import (linear_sum_assignment, OptimizeResult)
 from ._optimize import _check_unknown_options
 
@@ -8,6 +10,7 @@
 
 QUADRATIC_ASSIGNMENT_METHODS = ['faq', '2opt']
 
+
 def quadratic_assignment(A, B, method="faq", options=None):
     r"""
     Approximates solution to the quadratic assignment problem and
@@ -60,13 +63,24 @@ def quadratic_assignment(A, B, method="faq", options=None):
             ``partial_match[i, 1]`` of `B`. The array has shape ``(m, 2)``,
             where ``m`` is not greater than the number of nodes, :math:`n`.
 
-        rng : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}
-            If `seed` is None (or `np.random`), the `numpy.random.RandomState`
-            singleton is used.
-            If `seed` is an int, a new ``RandomState`` instance is used,
-            seeded with `seed`.
-            If `seed` is already a ``Generator`` or ``RandomState`` instance then
-            that instance is used.
+        rng : `numpy.random.Generator`, optional
+            Pseudorandom number generator state. When `rng` is None, a new
+            `numpy.random.Generator` is created using entropy from the
+            operating system. Types other than `numpy.random.Generator` are
+            passed to `numpy.random.default_rng` to instantiate a ``Generator``.
+
+            .. versionchanged:: 1.15.0
+                As part of the `SPEC-007 <https://scientific-python.org/specs/spec-0007/>`_
+                transition from use of `numpy.random.RandomState` to
+                `numpy.random.Generator` is occurring. Supplying
+                `np.random.RandomState` to this function will now emit a
+                `DeprecationWarning`. In SciPy 1.17 its use will raise an exception.
+                In addition relying on global state using `np.random.seed`
+                will emit a `FutureWarning`. In SciPy 1.17 the global random number
+                generator will no longer be used.
+                Use of an int-like seed will raise a `FutureWarning`, in SciPy 1.17 it
+                will be normalized via `np.random.default_rng` rather than
+                `np.random.RandomState`.
 
         For method-specific options, see
         :func:`show_options('quadratic_assignment') <show_options>`.
@@ -112,11 +126,12 @@ def quadratic_assignment(A, B, method="faq", options=None):
     --------
     >>> import numpy as np
     >>> from scipy.optimize import quadratic_assignment
+    >>> rng = np.random.default_rng()
     >>> A = np.array([[0, 80, 150, 170], [80, 0, 130, 100],
     ...               [150, 130, 0, 120], [170, 100, 120, 0]])
     >>> B = np.array([[0, 5, 2, 7], [0, 0, 3, 8],
     ...               [0, 0, 0, 3], [0, 0, 0, 0]])
-    >>> res = quadratic_assignment(A, B)
+    >>> res = quadratic_assignment(A, B, options={'rng': rng})
     >>> print(res)
          fun: 3260
      col_ind: [0 3 2 1]
@@ -159,7 +174,7 @@ def quadratic_assignment(A, B, method="faq", options=None):
     ...               [8, 5, 0, 2], [6, 1, 2, 0]])
     >>> B = np.array([[0, 1, 8, 4], [1, 0, 5, 2],
     ...               [8, 5, 0, 5], [4, 2, 5, 0]])
-    >>> res = quadratic_assignment(A, B)
+    >>> res = quadratic_assignment(A, B, options={'rng': rng})
     >>> print(res)
          fun: 178
      col_ind: [1 0 3 2]
@@ -170,7 +185,7 @@ def quadratic_assignment(A, B, method="faq", options=None):
 
     >>> guess = np.array([np.arange(len(A)), res.col_ind]).T
     >>> res = quadratic_assignment(A, B, method="2opt",
-    ...                            options = {'partial_guess': guess})
+    ...     options = {'rng': rng, 'partial_guess': guess})
     >>> print(res)
          fun: 176
      col_ind: [1 2 3 0]
@@ -186,10 +201,38 @@ def quadratic_assignment(A, B, method="faq", options=None):
                "2opt": _quadratic_assignment_2opt}
     if method not in methods:
         raise ValueError(f"method {method} must be in {methods}.")
+
+    _spec007_transition(options.get("rng", None))
     res = methods[method](A, B, **options)
     return res
 
 
+def _spec007_transition(rng):
+    if isinstance(rng, np.random.RandomState):
+        warnings.warn(
+            "Use of `RandomState` with `quadratic_assignment` is deprecated"
+            " and will result in an exception in SciPy 1.17",
+            DeprecationWarning,
+            stacklevel=2
+        )
+    if ((rng is None or rng is np.random) and
+            np.random.mtrand._rand._bit_generator._seed_seq is None):
+        warnings.warn(
+            "The NumPy global RNG was seeded by calling `np.random.seed`."
+            " From SciPy 1.17, this function will no longer use the global RNG.",
+            FutureWarning,
+            stacklevel=2
+        )
+    if isinstance(rng, numbers.Integral | np.integer):
+        warnings.warn(
+            "The behavior when the rng option is an integer is changing: the value"
+            " will be normalized using np.random.default_rng beginning in SciPy 1.17,"
+            " and the resulting Generator will be used to generate random numbers.",
+            FutureWarning,
+            stacklevel=2
+        )
+
+
 def _calc_score(A, B, perm):
     # equivalent to objective function but avoids matmul
     return np.sum(A * B[perm][:, perm])
@@ -284,13 +327,8 @@ def _quadratic_assignment_faq(A, B,
         ``partial_match[i, 1]`` of `B`. The array has shape ``(m, 2)``, where
         ``m`` is not greater than the number of nodes, :math:`n`.
 
-    rng : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
-        If `seed` is None (or `np.random`), the `numpy.random.RandomState`
-        singleton is used.
-        If `seed` is an int, a new ``RandomState`` instance is used,
-        seeded with `seed`.
-        If `seed` is already a ``Generator`` or ``RandomState`` instance then
-        that instance is used.
+    rng : {None, int, `numpy.random.Generator`}, optional
+        Pseudorandom number generator state. See `quadratic_assignment` for details.
     P0 : 2-D array, "barycenter", or "randomized" (default: "barycenter")
         Initial position. Must be a doubly-stochastic matrix [3]_.
 
@@ -312,7 +350,7 @@ def _quadratic_assignment_faq(A, B,
         Integer specifying the max number of Frank-Wolfe iterations performed.
     tol : float (default: 0.03)
         Tolerance for termination. Frank-Wolfe iteration terminates when
-        :math:`\frac{||P_{i}-P_{i+1}||_F}{\sqrt{m')}} \leq tol`,
+        :math:`\frac{||P_{i}-P_{i+1}||_F}{\sqrt{m'}} \leq tol`,
         where :math:`i` is the iteration number.
 
     Returns
@@ -343,34 +381,36 @@ def _quadratic_assignment_faq(A, B,
     As mentioned above, a barycenter initialization often results in a better
     solution than a single random initialization.
 
-    >>> from numpy.random import default_rng
-    >>> rng = default_rng()
+    >>> from scipy.optimize import quadratic_assignment
+    >>> import numpy as np
+    >>> rng = np.random.default_rng()
     >>> n = 15
     >>> A = rng.random((n, n))
     >>> B = rng.random((n, n))
-    >>> res = quadratic_assignment(A, B)  # FAQ is default method
+    >>> options = {"rng": rng}
+    >>> res = quadratic_assignment(A, B, options=options)  # FAQ is default method
     >>> print(res.fun)
-    46.871483385480545  # may vary
+    47.797048706380636  # may vary
 
-    >>> options = {"P0": "randomized"}  # use randomized initialization
+    >>> options = {"rng": rng, "P0": "randomized"}  # use randomized initialization
     >>> res = quadratic_assignment(A, B, options=options)
     >>> print(res.fun)
-    47.224831071310625 # may vary
+    47.37287069769966 # may vary
 
     However, consider running from several randomized initializations and
     keeping the best result.
 
     >>> res = min([quadratic_assignment(A, B, options=options)
     ...            for i in range(30)], key=lambda x: x.fun)
     >>> print(res.fun)
-    46.671852533681516 # may vary
+    46.55974835248574 # may vary
 
-    The '2-opt' method can be used to further refine the results.
+    The '2-opt' method can be used to attempt to refine the results.
 
-    >>> options = {"partial_guess": np.array([np.arange(n), res.col_ind]).T}
+    >>> options = {"partial_guess": np.array([np.arange(n), res.col_ind]).T, "rng": rng}
     >>> res = quadratic_assignment(A, B, method="2opt", options=options)
     >>> print(res.fun)
-    46.47160735721583 # may vary
+    46.55974835248574 # may vary
 
     References
     ----------
@@ -578,13 +618,8 @@ def _quadratic_assignment_2opt(A, B, maximize=False, rng=None,
     -------
     maximize : bool (default: False)
         Maximizes the objective function if ``True``.
-    rng : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
-        If `seed` is None (or `np.random`), the `numpy.random.RandomState`
-        singleton is used.
-        If `seed` is an int, a new ``RandomState`` instance is used,
-        seeded with `seed`.
-        If `seed` is already a ``Generator`` or ``RandomState`` instance then
-        that instance is used.
+    rng : {None, int, `numpy.random.Generator`}, optional
+        Pseudorandom number generator state. See `quadratic_assignment` for details.
     partial_match : 2-D array of integers, optional (default: None)
         Fixes part of the matching. Also known as a "seed" [2]_.