[DOC] make some doctests in user,reference pass pytest

1. Add `import numpy as np` in rst files
2. Update NumPy repr for array (whitespace)
3. Update bytearray representation
4. Fix tiny output formatting (`<class ...>`, etc)
5. Format tracebacks
6. Skip random number tests or some platform-dependent outputs
7. Add `<matplotlib. ... at 0x...>` or similar output lines where
   missing
8. Set seed

12 files changed, 69 insertions, 36 deletions

diff --git a/doc/source/reference/arrays.datetime.rst b/doc/source/reference/arrays.datetime.rst
index 63c93821b..37e0a5463 100644
--- a/doc/source/reference/arrays.datetime.rst
+++ b/doc/source/reference/arrays.datetime.rst
@@ -2,6 +2,9 @@
 
 .. _arrays.datetime:
 
+.. for doctest:
+    >>> import numpy as np
+
 ************************
 Datetimes and Timedeltas
 ************************
diff --git a/doc/source/reference/arrays.dtypes.rst b/doc/source/reference/arrays.dtypes.rst
index 8606bc8f1..aec5539e8 100644
--- a/doc/source/reference/arrays.dtypes.rst
+++ b/doc/source/reference/arrays.dtypes.rst
@@ -1,5 +1,8 @@
 .. currentmodule:: numpy
 
+.. for doctest:
+    >>> import numpy as np
+
 .. _arrays.dtypes:
 
 **********************************
diff --git a/doc/source/reference/arrays.ndarray.rst b/doc/source/reference/arrays.ndarray.rst
index 0f703b475..889fb103c 100644
--- a/doc/source/reference/arrays.ndarray.rst
+++ b/doc/source/reference/arrays.ndarray.rst
@@ -1,5 +1,8 @@
 .. currentmodule:: numpy
 
+.. for doctest:
+    >>> import numpy as np
+
 .. _arrays.ndarray:
 
 ******************************************
@@ -54,13 +57,13 @@ objects implementing the :class:`buffer` or :ref:`array
 
    >>> y = x[:,1]
    >>> y
-   array([2, 5])
+   array([2, 5], dtype=int32)
    >>> y[0] = 9 # this also changes the corresponding element in x
    >>> y
-   array([9, 5])
+   array([9, 5], dtype=int32)
    >>> x
    array([[1, 9, 3],
-          [4, 5, 6]])
+          [4, 5, 6]], dtype=int32)
 
 
 Constructing arrays
diff --git a/doc/source/reference/arrays.nditer.cython.rst b/doc/source/reference/arrays.nditer.cython.rst
index 43aad9927..9e51162f2 100644
--- a/doc/source/reference/arrays.nditer.cython.rst
+++ b/doc/source/reference/arrays.nditer.cython.rst
@@ -1,3 +1,6 @@
+.. for doctest:
+    >>> import numpy as np
+
 Putting the Inner Loop in Cython
 ================================
 
@@ -49,7 +52,7 @@ Here's how this looks.
     ...
     >>> a = np.arange(6).reshape(2,3)
     >>> sum_squares_py(a)
-    array(55.0)
+    array(55.)
     >>> sum_squares_py(a, axis=-1)
     array([  5.,  50.])
 
@@ -117,11 +120,11 @@ as our native Python/NumPy code did.
 
 .. admonition:: Example
 
-    >>> from sum_squares import sum_squares_cy
+    >>> from sum_squares import sum_squares_cy #doctest: +SKIP
     >>> a = np.arange(6).reshape(2,3)
-    >>> sum_squares_cy(a)
+    >>> sum_squares_cy(a) #doctest: +SKIP
     array(55.0)
-    >>> sum_squares_cy(a, axis=-1)
+    >>> sum_squares_cy(a, axis=-1) #doctest: +SKIP
     array([  5.,  50.])
 
 Doing a little timing in IPython shows that the reduced overhead and
diff --git a/doc/source/reference/arrays.nditer.rst b/doc/source/reference/arrays.nditer.rst
index 72a04f73e..89074d517 100644
--- a/doc/source/reference/arrays.nditer.rst
+++ b/doc/source/reference/arrays.nditer.rst
@@ -3,6 +3,7 @@
 .. for doctests
    The last section on Cython is 'included' at the end of this file. The tests
    for that section are disabled.
+   >>> import numpy as np
 
 .. _arrays.nditer:
 
@@ -489,9 +490,9 @@ reasons.
 
     >>> b = np.zeros((3,))
     >>> square([1,2,3], out=b)
-    array([ 1.,  4.,  9.])
+    array([1.,  4.,  9.])
     >>> b
-    array([ 1.,  4.,  9.])
+    array([1.,  4.,  9.])
 
     >>> square(np.arange(6).reshape(2,3), out=b)
     Traceback (most recent call last):
diff --git a/doc/source/reference/arrays.scalars.rst b/doc/source/reference/arrays.scalars.rst
index c691e802f..42fdd18bb 100644
--- a/doc/source/reference/arrays.scalars.rst
+++ b/doc/source/reference/arrays.scalars.rst
@@ -1,3 +1,6 @@
+.. for doctest:
+    >>> import numpy as np
+
 .. _arrays.scalars:
 
 *******
diff --git a/doc/source/reference/c-api/coremath.rst b/doc/source/reference/c-api/coremath.rst
index e129fdd77..2a851e937 100644
--- a/doc/source/reference/c-api/coremath.rst
+++ b/doc/source/reference/c-api/coremath.rst
@@ -1,3 +1,6 @@
+.. for doctest:
+   >>> import numpy as np
+
 NumPy core libraries
 ====================
 
diff --git a/doc/source/reference/distutils.rst b/doc/source/reference/distutils.rst
index f201ba668..33ebeb62c 100644
--- a/doc/source/reference/distutils.rst
+++ b/doc/source/reference/distutils.rst
@@ -1,3 +1,6 @@
+.. for doctest:
+    >>> import numpy as np
+
 **********************************
 Packaging (:mod:`numpy.distutils`)
 **********************************
@@ -188,6 +191,8 @@ Info are easily retrieved from the `get_info` function in
 
   >>> info = np.distutils.misc_util.get_info('npymath')
   >>> config.add_extension('foo', sources=['foo.c'], extra_info=info)
+  <numpy.distutils.extension.Extension('foo') at 0x...>
+
 
 An additional list of paths to look for .ini files can be given to `get_info`.
 
diff --git a/doc/source/reference/maskedarray.generic.rst b/doc/source/reference/maskedarray.generic.rst
index d3849c50d..4d8d9750a 100644
--- a/doc/source/reference/maskedarray.generic.rst
+++ b/doc/source/reference/maskedarray.generic.rst
@@ -433,7 +433,7 @@ and entries of the output masked array are masked wherever the corresponding
 input fall outside the validity domain::
 
    >>> x = ma.array([-1, 1, 0, 2, 3], mask=[0, 0, 0, 0, 1])
-   >>> np.log(x)
+   >>> np.log(x) # doctest: +SKIP
    masked_array(data=[--, 0.0, --, 0.6931471805599453, --],
                 mask=[ True, False,  True, False,  True],
           fill_value=1e+20)
@@ -467,7 +467,7 @@ Suppose now that we wish to print that same data, but with the missing values
 replaced by the average value.
 
    >>> print(mx.filled(mx.mean()))
-   [ 0.  1.  2.  3.  4.]
+   [0.  1.  2.  3.  4.]
 
 
 Numerical operations
diff --git a/doc/source/reference/random/generator.rst b/doc/source/reference/random/generator.rst
index 7934be98a..4a863ebf3 100644
--- a/doc/source/reference/random/generator.rst
+++ b/doc/source/reference/random/generator.rst
@@ -12,6 +12,9 @@ random values from useful distributions. The default BitGenerator used by
 can be changed by passing an instantized BitGenerator to ``Generator``.
 
 
+.. for doctest:
+    >>> import numpy as np
+
 .. autofunction:: default_rng
 
 .. autoclass:: Generator
@@ -71,7 +74,7 @@ By default, `Generator.permuted` returns a copy.  To operate in-place with
 `Generator.permuted`, pass the same array as the first argument *and* as
 the value of the ``out`` parameter.  For example,
 
-    >>> rng = np.random.default_rng()
+    >>> rng = np.random.default_rng(12345)
     >>> x = np.arange(0, 15).reshape(3, 5)
     >>> x
     array([[ 0,  1,  2,  3,  4],
@@ -79,9 +82,9 @@ the value of the ``out`` parameter.  For example,
            [10, 11, 12, 13, 14]])
     >>> y = rng.permuted(x, axis=1, out=x)
     >>> x
-    array([[ 1,  0,  2,  4,  3],  # random
-           [ 6,  7,  8,  9,  5],
-           [10, 14, 11, 13, 12]])
+    array([[ 4,  3,  0,  2,  1],
+           [ 9,  7,  6,  8,  5],
+           [10, 12, 13, 11, 14]])
 
 Note that when ``out`` is given, the return value is ``out``:
 
@@ -97,16 +100,16 @@ which dimension of the input array to use as the sequence. In the case of a
 two-dimensional array, ``axis=0`` will, in effect, rearrange the rows of the
 array, and  ``axis=1`` will rearrange the columns.  For example
 
-    >>> rng = np.random.default_rng()
+    >>> rng = np.random.default_rng(2345)
     >>> x = np.arange(0, 15).reshape(3, 5)
     >>> x
     array([[ 0,  1,  2,  3,  4],
            [ 5,  6,  7,  8,  9],
            [10, 11, 12, 13, 14]])
     >>> rng.permutation(x, axis=1)
-    array([[ 1,  3,  2,  0,  4],  # random
-           [ 6,  8,  7,  5,  9],
-           [11, 13, 12, 10, 14]])
+    array([[ 4,  2,  1,  3,  0],
+           [ 9,  7,  6,  8,  5],
+           [14, 12, 11, 13, 10]])
 
 Note that the columns have been rearranged "in bulk": the values within
 each column have not changed.
@@ -117,9 +120,9 @@ independently of the others.  Compare the following example of the use of
 `Generator.permuted` to the above example of `Generator.permutation`:
 
     >>> rng.permuted(x, axis=1)
-    array([[ 1,  0,  2,  4,  3],  # random
-           [ 5,  7,  6,  9,  8],
-           [10, 14, 12, 13, 11]])
+    array([[ 1,  2,  0,  3,  4],
+           [ 7,  9,  8,  6,  5],
+           [13, 11, 10, 14, 12]])
 
 In this example, the values within each row (i.e. the values along
 ``axis=1``) have been shuffled independently.  This is not a "bulk"
@@ -131,11 +134,11 @@ Shuffling non-NumPy sequences
 a sequence that is not a NumPy array, it shuffles that sequence in-place.
 For example,
 
-    >>> rng = np.random.default_rng()
+    >>> rng = np.random.default_rng(3456)
     >>> a = ['A', 'B', 'C', 'D', 'E']
     >>> rng.shuffle(a)  # shuffle the list in-place
     >>> a
-    ['B', 'D', 'A', 'E', 'C']  # random
+    ['B', 'E', 'A', 'D', 'C']
 
 Distributions
 =============
diff --git a/doc/source/reference/routines.polynomials.classes.rst b/doc/source/reference/routines.polynomials.classes.rst
index 5f575bed1..e36ef6e97 100644
--- a/doc/source/reference/routines.polynomials.classes.rst
+++ b/doc/source/reference/routines.polynomials.classes.rst
@@ -1,3 +1,6 @@
+.. for doctest:
+    >>> import numpy as np
+
 Using the Convenience Classes
 =============================
 
@@ -59,11 +62,11 @@ first is the coefficients, the second is the domain, and the third is the
 window::
 
    >>> p.coef
-   array([ 1.,  2.,  3.])
+   array([1., 2., 3.])
    >>> p.domain
-   array([-1.,  1.])
+   array([-1,  1])
    >>> p.window
-   array([-1.,  1.])
+   array([-1,  1])
 
 Printing a polynomial yields the polynomial expression in a more familiar
 format::
@@ -77,7 +80,7 @@ representation is also available (default on Windows). The polynomial string
 format can be toggled at the package-level with the 
 `~numpy.polynomial.set_default_printstyle` function::
 
-   >>> numpy.polynomial.set_default_printstyle('ascii')
+   >>> np.polynomial.set_default_printstyle('ascii')
    >>> print(p)
    1.0 + 2.0 x**1 + 3.0 x**2
 
@@ -137,9 +140,9 @@ Evaluation::
    array([  1.,   6.,  17.,  34.,  57.])
    >>> x = np.arange(6).reshape(3,2)
    >>> p(x)
-   array([[  1.,   6.],
-          [ 17.,  34.],
-          [ 57.,  86.]])
+   array([[ 1.,   6.],
+          [17.,  34.],
+          [57.,  86.]])
 
 Substitution:
 
@@ -294,7 +297,7 @@ polynomials up to degree 5 are plotted below.
     ...     ax = plt.plot(x, T.basis(i)(x), lw=2, label=f"$T_{i}$")
     ...
     >>> plt.legend(loc="upper left")
-    <matplotlib.legend.Legend object at 0x3b3ee10>
+    <matplotlib.legend.Legend object at 0x...>
     >>> plt.show()
 
 In the range -1 <= `x` <= 1 they are nice, equiripple functions lying between +/- 1.
@@ -309,7 +312,7 @@ The same plots over the range -2 <= `x` <= 2 look very different:
     ...     ax = plt.plot(x, T.basis(i)(x), lw=2, label=f"$T_{i}$")
     ...
     >>> plt.legend(loc="lower right")
-    <matplotlib.legend.Legend object at 0x3b3ee10>
+    <matplotlib.legend.Legend object at 0x...>
     >>> plt.show()
 
 As can be seen, the "good" parts have shrunk to insignificance. In using
@@ -335,12 +338,12 @@ illustrated below for a fit to a noisy sine curve.
     >>> y = np.sin(x) + np.random.normal(scale=.1, size=x.shape)
     >>> p = T.fit(x, y, 5)
     >>> plt.plot(x, y, 'o')
-    [<matplotlib.lines.Line2D object at 0x2136c10>]
+    [<matplotlib.lines.Line2D object at 0x...>]
     >>> xx, yy = p.linspace()
     >>> plt.plot(xx, yy, lw=2)
-    [<matplotlib.lines.Line2D object at 0x1cf2890>]
+    [<matplotlib.lines.Line2D object at 0x...>]
     >>> p.domain
-    array([ 0.        ,  6.28318531])
+    array([0.        ,  6.28318531])
     >>> p.window
     array([-1.,  1.])
     >>> plt.show()
diff --git a/doc/source/reference/routines.polynomials.rst b/doc/source/reference/routines.polynomials.rst
index 4aea963c0..75bed4694 100644
--- a/doc/source/reference/routines.polynomials.rst
+++ b/doc/source/reference/routines.polynomials.rst
@@ -1,5 +1,8 @@
 .. _routines.polynomial:
 
+.. for doctest:
+   >>> import numpy as np
+
 Polynomials
 ***********