Merge pull request #360 from jcapriot/random_generator

Update use of `numpy`'s random number generators.

.ci/environment_test.yml

@@ -4,21 +4,23 @@ channels:
 dependencies:
   - numpy>=1.22.4
   - scipy>=1.8
+
   # optionals
   - vtk>=6
   - pyvista
   - omf
   - matplotlib
+
   # documentation
   - sphinx
   - pydata-sphinx-theme==0.13.3
   - sphinx-gallery>=0.1.13
   - numpydoc>=1.5
   - jupyter
   - graphviz
-  - pymatsolver>=0.1.2
   - pillow
   - pooch
+
   # testing
   - sympy
   - pytest

discretize/base/base_mesh.py

@@ -2419,15 +2419,15 @@ def get_face_inner_product_deriv(
         >>> import numpy as np
         >>> import matplotlib as mpl
         >>> mpl.rcParams.update({'font.size': 14})
-        >>> np.random.seed(45)
+        >>> rng = np.random.default_rng(45)
         >>> mesh = TensorMesh([[(1, 4)], [(1, 4)]])
 
         Define a model, and a random vector to multiply the derivative with,
         then we grab the respective derivative function and calculate the
         sparse matrix,
 
-        >>> m = np.random.rand(mesh.nC)  # physical property parameters
-        >>> u = np.random.rand(mesh.nF)  # vector of shape (n_faces)
+        >>> m = rng.random(mesh.nC)  # physical property parameters
+        >>> u = rng.random(mesh.nF)  # vector of shape (n_faces)
         >>> Mf = mesh.get_face_inner_product(m)
         >>> F = mesh.get_face_inner_product_deriv(m)  # Function handle
         >>> dFdm_u = F(u)
@@ -2458,8 +2458,8 @@ def get_face_inner_product_deriv(
         function handle :math:`\mathbf{F}(\mathbf{u})` and plot the evaluation
         of this function on a spy plot.
 
-        >>> m = np.random.rand(mesh.nC, 3)  # anisotropic physical property parameters
-        >>> u = np.random.rand(mesh.nF)     # vector of shape (n_faces)
+        >>> m = rng.random((mesh.nC, 3))  # anisotropic physical property parameters
+        >>> u = rng.random(mesh.nF)     # vector of shape (n_faces)
         >>> Mf = mesh.get_face_inner_product(m)
         >>> F = mesh.get_face_inner_product_deriv(m)  # Function handle
         >>> dFdm_u = F(u)
@@ -2602,14 +2602,14 @@ def get_edge_inner_product_deriv(
         >>> import numpy as np
         >>> import matplotlib as mpl
         >>> mpl.rcParams.update({'font.size': 14})
-        >>> np.random.seed(45)
+        >>> rng = np.random.default_rng(45)
         >>> mesh = TensorMesh([[(1, 4)], [(1, 4)]])
 
         Next we create a random isotropic model vector, and a random vector to multiply
         the derivative with (for illustration purposes).
 
-        >>> m = np.random.rand(mesh.nC)  # physical property parameters
-        >>> u = np.random.rand(mesh.nF)  # vector of shape (n_edges)
+        >>> m = rng.random(mesh.nC)  # physical property parameters
+        >>> u = rng.random(mesh.nF)  # vector of shape (n_edges)
         >>> Me = mesh.get_edge_inner_product(m)
         >>> F = mesh.get_edge_inner_product_deriv(m)  # Function handle
         >>> dFdm_u = F(u)
@@ -2640,8 +2640,8 @@ def get_edge_inner_product_deriv(
         function handle :math:`\mathbf{F}(\mathbf{u})` and plot the evaluation
         of this function on a spy plot.
 
-        >>> m = np.random.rand(mesh.nC, 3)  # physical property parameters
-        >>> u = np.random.rand(mesh.nF)     # vector of shape (n_edges)
+        >>> m = rng.random((mesh.nC, 3))  # physical property parameters
+        >>> u = rng.random(mesh.nF)     # vector of shape (n_edges)
         >>> Me = mesh.get_edge_inner_product(m)
         >>> F = mesh.get_edge_inner_product_deriv(m)  # Function handle
         >>> dFdm_u = F(u)
@@ -4128,9 +4128,9 @@ def get_interpolation_matrix(
         >>> from discretize import TensorMesh
         >>> import numpy as np
         >>> import matplotlib.pyplot as plt
-        >>> np.random.seed(14)
+        >>> rng = np.random.default_rng(14)
 
-        >>> locs = np.random.rand(50)*0.8+0.1
+        >>> locs = rng.random(50)*0.8+0.1
         >>> dense = np.linspace(0, 1, 200)
         >>> fun = lambda x: np.cos(2*np.pi*x)
 

discretize/mixins/mpl_mod.py

@@ -470,7 +470,7 @@ def plot_slice(
 
         >>> from matplotlib import pyplot as plt
         >>> import discretize
-        >>> from pymatsolver import Solver
+        >>> from scipy.sparse.linalg import spsolve
         >>> hx = [(5, 2, -1.3), (2, 4), (5, 2, 1.3)]
         >>> hy = [(2, 2, -1.3), (2, 6), (2, 2, 1.3)]
         >>> hz = [(2, 2, -1.3), (2, 6), (2, 2, 1.3)]
@@ -482,7 +482,7 @@ def plot_slice(
         >>> q[[4, 4], [4, 4], [2, 6]]=[-1, 1]
         >>> q = discretize.utils.mkvc(q)
         >>> A = M.face_divergence * M.cell_gradient
-        >>> b = Solver(A) * (q)
+        >>> b = spsolve(A, q)
 
         and finaly, plot the vector values of the result, which are defined on faces
 

discretize/tests.py

@@ -78,11 +78,10 @@
     "You break it, you fix it.",
 ]
 
-# Initiate random number generator
-rng = np.random.default_rng()
+_happiness_rng = np.random.default_rng()
 
 
-def setup_mesh(mesh_type, nC, nDim):
+def setup_mesh(mesh_type, nC, nDim, random_seed=None):
     """Generate arbitrary mesh for testing.
 
     For the mesh type, number of cells along each axis and dimension specified,
@@ -100,19 +99,25 @@ def setup_mesh(mesh_type, nC, nDim):
         number of base mesh cells and must be a power of 2.
     nDim : int
         The dimension of the mesh. Must be 1, 2 or 3.
+    random_seed : numpy.random.Generator, int, optional
+        If ``random`` is in `mesh_type`, this is the random number generator to use for
+        creating that random mesh. If an integer or None it is used to seed a new
+        `numpy.random.default_rng`.
 
     Returns
     -------
     discretize.base.BaseMesh
         A discretize mesh of class specified by the input argument *mesh_type*
     """
+    if "random" in mesh_type:
+        rng = np.random.default_rng(random_seed)
     if "TensorMesh" in mesh_type:
         if "uniform" in mesh_type:
             h = [nC, nC, nC]
         elif "random" in mesh_type:
-            h1 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
-            h2 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
-            h3 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
+            h1 = rng.random(nC) * nC * 0.5 + nC * 0.5
+            h2 = rng.random(nC) * nC * 0.5 + nC * 0.5
+            h3 = rng.random(nC) * nC * 0.5 + nC * 0.5
             h = [hi / np.sum(hi) for hi in [h1, h2, h3]]  # normalize
         else:
             raise Exception("Unexpected mesh_type")
@@ -127,14 +132,14 @@ def setup_mesh(mesh_type, nC, nDim):
             else:
                 h = [nC, nC, nC]
         elif "random" in mesh_type:
-            h1 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
+            h1 = rng.random(nC) * nC * 0.5 + nC * 0.5
             if "symmetric" in mesh_type:
                 h2 = [
                     2 * np.pi,
                 ]
             else:
-                h2 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
-            h3 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
+                h2 = rng.random(nC) * nC * 0.5 + nC * 0.5
+            h3 = rng.random(nC) * nC * 0.5 + nC * 0.5
             h = [hi / np.sum(hi) for hi in [h1, h2, h3]]  # normalize
             h[1] = h[1] * 2 * np.pi
         else:
@@ -179,9 +184,9 @@ def setup_mesh(mesh_type, nC, nDim):
         if "uniform" in mesh_type or "notatree" in mesh_type:
             h = [nC, nC, nC]
         elif "random" in mesh_type:
-            h1 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
-            h2 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
-            h3 = np.random.rand(nC) * nC * 0.5 + nC * 0.5
+            h1 = rng.random(nC) * nC * 0.5 + nC * 0.5
+            h2 = rng.random(nC) * nC * 0.5 + nC * 0.5
+            h3 = rng.random(nC) * nC * 0.5 + nC * 0.5
             h = [hi / np.sum(hi) for hi in [h1, h2, h3]]  # normalize
         else:
             raise Exception("Unexpected mesh_type")
@@ -216,13 +221,13 @@ class for the given operator. Within the test class, the user sets the parameter
 
     OrderTest inherits from :class:`unittest.TestCase`.
 
-    Parameters
+    Attributes
     ----------
     name : str
         Name the convergence test
     meshTypes : list of str
         List denoting the mesh types on which the convergence will be tested.
-        List entries must of the list {'uniformTensorMesh', 'randomTensorMesh',
+        List entries must be of the list {'uniformTensorMesh', 'randomTensorMesh',
         'uniformCylindricalMesh', 'randomCylindricalMesh', 'uniformTree', 'randomTree',
         'uniformCurv', 'rotateCurv', 'sphereCurv'}
     expectedOrders : float or list of float (default = 2.0)
@@ -236,6 +241,10 @@ class for the given operator. Within the test class, the user sets the parameter
         for the meshes used in the convergence test; e.g. [4, 8, 16, 32]
     meshDimension : int
         Mesh dimension. Must be 1, 2 or 3
+    random_seed : numpy.random.Generator, int, optional
+        If ``random`` is in `mesh_type`, this is the random number generator
+        used generate the random meshes, if an ``int`` or ``None``, it used to seed
+        a new `numpy.random.default_rng`.
 
     Notes
     -----
@@ -302,6 +311,7 @@ class for the given operator. Within the test class, the user sets the parameter
     meshTypes = ["uniformTensorMesh"]
     _meshType = meshTypes[0]
     meshDimension = 3
+    random_seed = None
 
     def setupMesh(self, nC):
         """Generate mesh and set as current mesh for testing.
@@ -316,7 +326,9 @@ def setupMesh(self, nC):
         Float
             Maximum cell width for the mesh
         """
-        mesh, max_h = setup_mesh(self._meshType, nC, self.meshDimension)
+        mesh, max_h = setup_mesh(
+            self._meshType, nC, self.meshDimension, random_seed=self.random_seed
+        )
         self.M = mesh
         return max_h
 
@@ -335,7 +347,7 @@ def getError(self):
         """
         return 1.0
 
-    def orderTest(self):
+    def orderTest(self, random_seed=None):
         """Perform an order test.
 
         For number of cells specified in meshSizes setup mesh, call getError
@@ -360,6 +372,9 @@ def orderTest(self):
                 "expectedOrders must have the same length as the meshTypes"
             )
 
+        if random_seed is not None:
+            self.random_seed = random_seed
+
         def test_func(n_cells):
             max_h = self.setupMesh(n_cells)
             err = self.getError()
@@ -496,9 +511,9 @@ class in older versions of `discretize`.
             np.testing.assert_allclose(orders[-1], expected_order, rtol=rtol)
         elif test_type == "all":
             np.testing.assert_allclose(orders, expected_order, rtol=rtol)
-        print(np.random.choice(happiness))
+        print(_happiness_rng.choice(happiness))
     except AssertionError as err:
-        print(np.random.choice(sadness))
+        print(_happiness_rng.choice(sadness))
         raise err
 
     return orders
@@ -552,6 +567,7 @@ def check_derivative(
     tolerance=0.85,
     eps=1e-10,
     ax=None,
+    random_seed=None,
 ):
     """Perform a basic derivative check.
 
@@ -560,8 +576,8 @@ def check_derivative(
 
     Parameters
     ----------
-    fctn : function
-        Function handle
+    fctn : callable
+        The function to test.
     x0 : numpy.ndarray
         Point at which to check derivative
     num : int, optional
@@ -570,7 +586,7 @@ def check_derivative(
         If *True*, plot the convergence of the approximation of the derivative
     dx : numpy.ndarray, optional
         Step direction. By default, this parameter is set to *None* and a random
-        step direction is chosen.
+        step direction is chosen using `rng`.
     expectedOrder : int, optional
         The expected order of convergence for the numerical derivative
     tolerance : float, optional
@@ -580,6 +596,10 @@ def check_derivative(
     ax : matplotlib.pyplot.Axes, optional
         An axis object for the convergence plot if *plotIt = True*.
         Otherwise, the function will create a new axis.
+    random_seed : numpy.random.Generator, int, optional
+        If `dx` is ``None``, this is the random number generator to use for
+        generating a step direction. If an integer or None, it is used to seed
+        a new `numpy.random.default_rng`.
 
     Returns
     -------
@@ -591,10 +611,11 @@ def check_derivative(
     >>> from discretize import tests, utils
     >>> import numpy as np
     >>> import matplotlib.pyplot as plt
+    >>> rng = np.random.default_rng(786412)
 
     >>> def simplePass(x):
     ...     return np.sin(x), utils.sdiag(np.cos(x))
-    >>> passed = tests.check_derivative(simplePass, np.random.randn(5))
+    >>> passed = tests.check_derivative(simplePass, rng.standard_normal(5), random_seed=rng)
     ==================== check_derivative ====================
     iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
     ---------------------------------------------------------
@@ -611,6 +632,7 @@ def check_derivative(
     __tracebackhide__ = True
     # matplotlib is a soft dependencies for discretize,
     # lazy-loaded to decrease load time of discretize.
+
     try:
         import matplotlib
         import matplotlib.pyplot as plt
@@ -627,7 +649,8 @@ def check_derivative(
     x0 = mkvc(x0)
 
     if dx is None:
-        dx = np.random.randn(len(x0))
+        rng = np.random.default_rng(random_seed)
+        dx = rng.standard_normal(len(x0))
 
     h = np.logspace(-1, -num, num)
     E0 = np.ones(h.shape)
@@ -693,14 +716,17 @@ def _plot_it(axes, passed):
             # Thus it has no higher order derivatives.
             pass
         else:
-            test = np.mean(order1) > tolerance * expectedOrder
+            order_mean = np.mean(order1)
+            expected = tolerance * expectedOrder
+            test = order_mean > expected
             if not test:
                 raise AssertionError(
-                    f"\n Order mean {np.mean(order1)} is not greater than"
-                    f" {tolerance} of the expected order {expectedOrder}."
+                    f"\n Order mean {order_mean} is not greater than"
+                    f" {expected} = tolerance: {tolerance} "
+                    f"* expected order: {expectedOrder}."
                 )
         print("{0!s} PASS! {1!s}".format("=" * 25, "=" * 25))
-        print(np.random.choice(happiness) + "\n")
+        print(_happiness_rng.choice(happiness) + "\n")
         if plotIt:
             _plot_it(ax, True)
     except AssertionError as err:
@@ -709,7 +735,7 @@ def _plot_it(axes, passed):
                 "*" * 57, "<" * 25, ">" * 25, "*" * 57
             )
         )
-        print(np.random.choice(sadness) + "\n")
+        print(_happiness_rng.choice(sadness) + "\n")
         if plotIt:
             _plot_it(ax, False)
         raise err
@@ -773,6 +799,7 @@ def assert_isadjoint(
     rtol=1e-6,
     atol=0.0,
     assert_error=True,
+    random_seed=None,
 ):
     r"""Do a dot product test for the forward operator and its adjoint operator.
 
@@ -823,6 +850,9 @@ def assert_isadjoint(
         assertion error if failed). If set to False, the result of the test is
         returned as boolean and a message is printed.
 
+    random_seed : numpy.random.Generator, int, optional
+        The random number generator to use for the adjoint test. If an integer or None
+        it is used to seed a new `numpy.random.default_rng`.
 
     Returns
     -------
@@ -837,6 +867,8 @@ def assert_isadjoint(
     """
     __tracebackhide__ = True
 
+    rng = np.random.default_rng(random_seed)
+
     def random(size, iscomplex):
         """Create random data of size and dtype of <size>."""
         out = rng.standard_normal(size)