gh-38250: provide compatibility with numpy 2.0

    
This will fix #38242

Basically, most of the problems are with numpy switching to a different
way to display its data,
showing type. We either use a compatibility version, or add ellipses in
doctests.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x ] The title is concise and informative.
- [x ] The description explains in detail what this PR is about.
- [x ] I have linked a relevant issue or discussion.
- [x ] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #38250
Reported by: Dima Pasechnik
Reviewer(s):

src/doc/en/faq/faq-usage.rst

@@ -324,7 +324,7 @@ ints. For example::
     sage: RealNumber = float; Integer = int
     sage: from scipy import stats
     sage: stats.ttest_ind([1,2,3,4,5], [2,3,4,5,.6])
-    Ttest...Result(statistic=0.0767529..., pvalue=0.940704...)
+    Ttest...Result(statistic=...0.0767529..., pvalue=...0.940704...)
     sage: stats.uniform(0,15).ppf([0.5,0.7])
     array([  7.5,  10.5])
 

src/doc/en/prep/Quickstarts/Statistics-and-Distributions.rst

@@ -24,6 +24,8 @@ NumPy provides, for example, functions to compute the arithmetic mean and
 the standard deviation::
 
     sage: import numpy as np
+    sage: if int(np.version.short_version[0]) > 1:
+    ....:     np.set_printoptions(legacy="1.25")
     sage: np.mean([1, 2, 3, 5])
     2.75
 

src/doc/en/thematic_tutorials/numerical_sage/numpy.rst

@@ -7,6 +7,8 @@ import it.
 ::
 
     sage: import numpy
+    sage: if int(numpy.version.short_version[0]) > 1:
+    ....:     numpy.set_printoptions(legacy="1.25")  # to ensure numpy 2.0 compatibility
 
 The basic object of computation in NumPy is an array. It is simple to
 create an array.

src/doc/it/faq/faq-usage.rst

@@ -305,7 +305,7 @@ anziché Integer di Sage. Ad esempio::
     sage: RealNumber = float; Integer = int
     sage: from scipy import stats
     sage: stats.ttest_ind([1,2,3,4,5], [2,3,4,5,.6])
-    Ttest...Result(statistic=0.0767529..., pvalue=0.940704...)
+    Ttest...Result(statistic=...0.0767529..., pvalue=...0.940704...)
     sage: stats.uniform(0,15).ppf([0.5,0.7])
     array([  7.5,  10.5])
 

src/sage/arith/misc.py

@@ -2277,7 +2277,7 @@ def power_mod(a, n, m):
 
         sage: from numpy import int32                                                   # needs numpy
         sage: power_mod(int32(2), int32(390), int32(391))                               # needs numpy
-        285
+        ...285...
         sage: from gmpy2 import mpz
         sage: power_mod(mpz(2), mpz(390), mpz(391))
         mpz(285)

src/sage/calculus/interpolators.pyx

@@ -27,6 +27,9 @@ Development supported by NSF award No. 0702939.
 import numpy as np
 cimport numpy as np
 
+if int(np.version.short_version[0]) > 1:
+    np.set_printoptions(legacy="1.25")
+
 from math import pi
 cdef double TWOPI = 2*pi
 

src/sage/calculus/riemann.pyx

@@ -44,6 +44,9 @@ from sage.calculus.integration import numerical_integral
 import numpy as np
 cimport numpy as np
 
+if int(np.version.short_version[0]) > 1:
+    np.set_printoptions(legacy="1.25")
+
 from math import pi
 from math import sin
 from math import cos

src/sage/combinat/fully_packed_loop.py

@@ -66,6 +66,9 @@ def _make_color_list(n, colors=None, color_map=None, randomize=False):
     r"""
     TESTS::
 
+        sage: import numpy as np
+        sage: if int(np.version.short_version[0]) > 1:
+        ....:     np.set_printoptions(legacy="1.25")
         sage: from sage.combinat.fully_packed_loop import _make_color_list
         sage: _make_color_list(5)
         sage: _make_color_list(5, ['blue', 'red'])

src/sage/env.py

@@ -296,7 +296,7 @@ def sage_include_directories(use_sources=False):
         sage: import sage.env
         sage: sage.env.sage_include_directories()
         ['...',
-         '.../numpy/core/include',
+         '.../numpy/...core/include',
          '.../include/python...']
 
     To check that C/C++ files are correctly found, we verify that we can

src/sage/functions/special.py

@@ -217,6 +217,9 @@ class SphericalHarmonic(BuiltinFunction):
         sage: spherical_harmonic(1, 1, pi/2, pi).n()  # abs tol 1e-14                   # needs sage.symbolic
         0.345494149471335
         sage: from scipy.special import sph_harm  # NB: arguments x and y are swapped   # needs scipy
+        sage: import numpy as np                                                        # needs scipy
+        sage: if int(np.version.short_version[0]) > 1:                                  # needs scipy
+        ....:     np.set_printoptions(legacy="1.25")                                    # needs scipy
         sage: sph_harm(1, 1, pi.n(), (pi/2).n())  # abs tol 1e-14                       # needs scipy sage.symbolic
         (0.3454941494713355-4.231083042742082e-17j)
 

src/sage/graphs/generators/classical_geometries.py

@@ -1315,7 +1315,7 @@ def CossidentePenttilaGraph(q):
 
     from sage.libs.gap.libgap import libgap
     adj_list = libgap.function_factory("""function(q)
-        local z, e, so, G, nu, G1, G0, B, T, s, O1, O2, x;
+        local z, e, so, G, nu, G1, G0, B, T, s, O1, O2, x, sqo;
         LoadPackage("grape");
         G0:=SO(3,q^2);
         so:=GeneratorsOfGroup(G0);

src/sage/graphs/generic_graph.py

@@ -18420,6 +18420,9 @@ def shortest_path_all_pairs(self, by_weight=False, algorithm=None,
                 M = self.adjacency_matrix(vertices=int_to_vertex)
 
             # We call the Floyd-Warshall method from SciPy
+            import numpy  # to ensure numpy 2.0 compatibility
+            if int(numpy.version.short_version[0]) > 1:
+                numpy.set_printoptions(legacy="1.25")
             from numpy import array as np_array
             from scipy.sparse.csgraph import floyd_warshall
             dd, pp = floyd_warshall(np_array(M), directed=self.is_directed(),

src/sage/graphs/graph_generators_pyx.pyx

@@ -57,7 +57,7 @@ def RandomGNP(n, p, bint directed=False, bint loops=False, seed=None):
 
         sage: from numpy import mean                                                    # needs numpy
         sage: abs(mean([RandomGNP(200, .2).density() for i in range(30)]) - .2) < .001  # needs numpy
-        True
+        ...True...
         sage: RandomGNP(150, .2, loops=True)
         Traceback (most recent call last):
         ...

src/sage/matrix/matrix1.pyx

@@ -705,12 +705,15 @@ cdef class Matrix(Matrix0):
         Type ``numpy.typecodes`` for a list of the possible
         typecodes::
 
-            sage: import numpy                                                          # needs numpy
-            sage: sorted(numpy.typecodes.items())                                       # needs numpy
+            sage: import numpy                            # needs numpy
+            sage: numpy.typecodes.items()                 # needs numpy # random
             [('All', '?bhilqpBHILQPefdgFDGSUVOMm'), ('AllFloat', 'efdgFDG'),
-             ('AllInteger', 'bBhHiIlLqQpP'), ('Character', 'c'), ('Complex', 'FDG'),
-             ('Datetime', 'Mm'), ('Float', 'efdg'), ('Integer', 'bhilqp'),
-             ('UnsignedInteger', 'BHILQP')]
+            ...
+
+        For instance, you can see possibilities for real floating point numbers::
+
+            sage: numpy.typecodes['Float']                # needs numpy
+            'efdg'
 
         Alternatively, numpy automatically calls this function (via
         the magic :meth:`__array__` method) to convert Sage matrices

src/sage/modules/free_module_element.pyx

@@ -555,7 +555,11 @@ def vector(arg0, arg1=None, arg2=None, sparse=None, immutable=False):
         R = None
 
     try:
+        import numpy
         from numpy import ndarray
+        if int(numpy.version.short_version[0]) > 1:
+            numpy.set_printoptions(legacy="1.25")
+
     except ImportError:
         pass
     else:
@@ -1188,7 +1192,11 @@ cdef class FreeModuleElement(Vector):   # abstract base class
             over Rational Field to numpy array of type <... 'float'>:
             setting an array element with a sequence.
         """
+        import numpy
         from numpy import array
+        if int(numpy.version.short_version[0]) > 1:
+            numpy.set_printoptions(legacy="1.25")
+
         try:
             return array(self, dtype=dtype)
         except ValueError as e:

src/sage/numerical/optimize.py

@@ -155,6 +155,10 @@ def find_root(f, a, b, xtol=10e-13, rtol=2.0**-50, maxiter=100, full_output=Fals
         b = max(s_1, s_2)
 
     import scipy.optimize
+    import numpy
+    if int(numpy.version.short_version[0]) > 1:
+        numpy.set_printoptions(legacy="1.25")
+
     g = lambda x: float(f(x))
     brentqRes = scipy.optimize.brentq(g, a, b,
                                  full_output=full_output, xtol=xtol, rtol=rtol, maxiter=maxiter)
@@ -288,6 +292,10 @@ def find_local_minimum(f, a, b, tol=1.48e-08, maxfun=500):
     a = float(a)
     b = float(b)
     import scipy.optimize
+    import numpy
+    if int(numpy.version.short_version[0]) > 1:
+        numpy.set_printoptions(legacy="1.25")
+
     xmin, fval, iter, funcalls = scipy.optimize.fminbound(f, a, b, full_output=1, xtol=tol, maxfun=maxfun)
     return fval, xmin
 
@@ -376,6 +384,8 @@ def minimize(func, x0, gradient=None, hessian=None, algorithm="default",
         sage: def rosen(x):  # The Rosenbrock function
         ....:    return sum(100.0r*(x[1r:]-x[:-1r]**2.0r)**2.0r + (1r-x[:-1r])**2.0r)
         sage: import numpy
+        sage: if int(numpy.version.short_version[0]) > 1:
+        ....:     numpy.set_printoptions(legacy="1.25")
         sage: from numpy import zeros
         sage: def rosen_der(x):
         ....:    xm = x[1r:-1r]
@@ -393,6 +403,9 @@ def minimize(func, x0, gradient=None, hessian=None, algorithm="default",
     from sage.structure.element import Expression
     from sage.ext.fast_callable import fast_callable
     import numpy
+    if int(numpy.version.short_version[0]) > 1:
+        numpy.set_printoptions(legacy="1.25")
+
     from scipy import optimize
     if isinstance(func, Expression):
         var_list = func.variables()
@@ -530,6 +543,8 @@ def minimize_constrained(func,cons,x0,gradient=None,algorithm='default', **args)
     from sage.structure.element import Expression
     from sage.ext.fast_callable import fast_callable
     import numpy
+    if int(numpy.version.short_version[0]) > 1:
+        numpy.set_printoptions(legacy="1.25")
     from scipy import optimize
     function_type = type(lambda x,y: x+y)
 
@@ -652,6 +667,8 @@ def find_fit(data, model, initial_guess=None, parameters=None, variables=None, s
     ``lmdif`` and ``lmder`` algorithms.
     """
     import numpy
+    if int(numpy.version.short_version[0]) > 1:
+        numpy.set_printoptions(legacy="1.25")
 
     if not isinstance(data, numpy.ndarray):
         try:

src/sage/plot/arrow.py

@@ -53,6 +53,9 @@ def get_minmax_data(self):
 
         EXAMPLES::
 
+            sage: import numpy  # to ensure numpy 2.0 compatibility
+            sage: if int(numpy.version.short_version[0]) > 1:
+            ....:     numpy.set_printoptions(legacy="1.25")
             sage: from sage.plot.arrow import CurveArrow
             sage: b = CurveArrow(path=[[(0,0),(.5,.5),(1,0)],[(.5,1),(0,0)]],
             ....:                options={})

src/sage/plot/complex_plot.pyx

@@ -1282,15 +1282,13 @@ def rgb_to_hls(rgb):
         raise ValueError("Last dimension of input array must be 3; "
                          "shape {} was found.".format(rgb.shape))
     in_shape = rgb.shape
-    rgb = np.array(
-        rgb, copy=False, dtype=np.dtype(float), ndmin=2
-    )
+    rgb = np.asarray(rgb, dtype=np.dtype(float))
     rgb_max = rgb.max(-1)
     rgb_min = rgb.min(-1)
     l = (rgb_max + rgb_min)/2.0  # lightness
 
     hls = np.zeros_like(rgb)
-    delta = rgb.ptp(-1)
+    delta = np.ptp(rgb, -1)
     s = np.zeros_like(delta)
 
     ipos = delta > 0

src/sage/plot/graphics.py

@@ -2764,7 +2764,7 @@ def matplotlib(self, filename=None,
             sage: xmin, xmax = sub.get_xlim()
             sage: ymin, ymax = sub.get_ylim()
             sage: xmin > xmax, ymin > ymax
-            (True, True)
+            (...True..., ...True...)
         """
         if not isinstance(ticks, (list, tuple)):
             ticks = (ticks, None)

src/sage/plot/histogram.py

@@ -92,6 +92,8 @@ def get_minmax_data(self):
             {'xmax': 10.0, 'xmin': 3.0, 'ymax': 0.476190476190..., 'ymin': 0}
         """
         import numpy
+        if int(numpy.version.short_version[0]) > 1:
+            numpy.set_printoptions(legacy="1.25")
 
         # Extract these options (if they are not None) and pass them to
         # histogram()

src/sage/plot/multigraphics.py

@@ -784,6 +784,9 @@ def _add_subplot(self, figure, index, **options):
             True
             sage: G.position(1)
             (0.2, 0.3, 0.4, 0.1)
+            sage: import numpy  # to ensure numpy 2.0 compatibility
+            sage: if int(numpy.version.short_version[0]) > 1:
+            ....:     numpy.set_printoptions(legacy="1.25")
             sage: ax1.get_position().bounds  # tol 1.0e-13
             (0.2, 0.3, 0.4000000000000001, 0.10000000000000003)
 
@@ -1295,6 +1298,9 @@ def position(self, index):
             sage: g1 = plot(sin(x), (x, -pi, pi))
             sage: g2 = circle((0,1), 1.)
             sage: G = graphics_array([g1, g2])
+            sage: import numpy  # to ensure numpy 2.0 compatibility
+            sage: if int(numpy.version.short_version[0]) > 1:
+            ....:     numpy.set_printoptions(legacy="1.25")
             sage: G.position(0)  # tol 5.0e-3
             (0.025045451349937315,
              0.03415488992713045,

src/sage/plot/plot3d/implicit_surface.pyx

@@ -468,6 +468,7 @@ cdef class MarchingCubesTriangles(MarchingCubes):
             sage: cube_marcher.y_vertices.tolist()
             [[[<1.0, 0.5, 0.0>, None]], [[None, None]]]
             sage: cube_marcher.x_vertices.any() # This shouldn't affect the X vertices.
+            ...
         """
         (self.y_vertices, self.y_vertices_swapped) = \
             (self.y_vertices_swapped, self.y_vertices)
@@ -574,6 +575,7 @@ cdef class MarchingCubesTriangles(MarchingCubes):
             sage: cube_marcher.x_vertices.tolist()
             [[None, None], [None, <1.5, 1.0, 1.0>]]
             sage: cube_marcher.y_vertices.any() or cube_marcher.z_vertices.any() # This shouldn't affect the Y or Z vertices.
+            ...
         """
         cdef bint has_prev = (_prev is not None)
         cdef bint has_next = (_next is not None)

src/sage/plot/scatter_plot.py

@@ -63,9 +63,9 @@ def get_minmax_data(self):
             sage: s = scatter_plot([[0,1],[2,4],[3.2,6]])
             sage: d = s.get_minmax_data()
             sage: d['xmin']
-            0.0
+            ...0.0...
             sage: d['ymin']
-            1.0
+            ...1.0...
         """
         return {'xmin': self.xdata.min(),
                 'xmax': self.xdata.max(),

src/sage/plot/streamline_plot.py

@@ -71,6 +71,9 @@ def get_minmax_data(self):
         EXAMPLES::
 
             sage: x, y = var('x y')
+            sage: import numpy  # to ensure numpy 2.0 compatibility
+            sage: if int(numpy.version.short_version[0]) > 1:
+            ....:     numpy.set_printoptions(legacy="1.25")
             sage: d = streamline_plot((.01*x, x+y), (x,10,20), (y,10,20))[0].get_minmax_data()
             sage: d['xmin']
             10.0

src/sage/probability/probability_distribution.pyx

@@ -140,6 +140,9 @@ cdef class ProbabilityDistribution:
              1.8,
              2.0]
         """
+        import numpy as np
+        if int(np.version.short_version[0]) > 1:
+            np.set_printoptions(legacy="1.25")
         import pylab
         ell = [float(self.get_random_element()) for _ in range(num_samples)]
         S = pylab.hist(ell, bins, density=True)

src/sage/rings/complex_double.pyx

@@ -2558,10 +2558,10 @@ cdef class ComplexToCDF(Morphism):
 
         sage: # needs numpy
         sage: import numpy
-        sage: f = CDF.coerce_map_from(numpy.complex_)
-        sage: f(numpy.complex_(I))
+        sage: f = CDF.coerce_map_from(numpy.complex128)
+        sage: f(numpy.complex128(I))
         1.0*I
-        sage: f(numpy.complex_(I)).parent()
+        sage: f(numpy.complex128(I)).parent()
         Complex Double Field
     """
     def __init__(self, R):
@@ -2578,7 +2578,7 @@ cdef class ComplexToCDF(Morphism):
         EXAMPLES::
 
             sage: import numpy                                                          # needs numpy
-            sage: CDF(numpy.complex_(I))    # indirect doctest                          # needs numpy
+            sage: CDF(numpy.complex128(I))    # indirect doctest                          # needs numpy
             1.0*I
         """
         cdef ComplexDoubleElement z = <ComplexDoubleElement>ComplexDoubleElement.__new__(ComplexDoubleElement)
@@ -2592,7 +2592,7 @@ cdef class ComplexToCDF(Morphism):
         EXAMPLES::
 
             sage: import numpy                                                          # needs numpy
-            sage: f = sage.rings.complex_double.ComplexToCDF(numpy.complex_)            # needs numpy
+            sage: f = sage.rings.complex_double.ComplexToCDF(numpy.complex128)          # needs numpy
             sage: f._repr_type()                                                        # needs numpy
             'Native'
         """

src/sage/rings/integer.pyx

@@ -593,6 +593,8 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
         Test comparisons with numpy types (see :issue:`13386` and :issue:`18076`)::
 
             sage: import numpy                                                          # needs numpy
+            sage: if int(numpy.version.short_version[0]) > 1:                           # needs numpy
+            ....:     numpy.set_printoptions(legacy="1.25")                             # needs numpy
             sage: numpy.int8('12') == 12                                                # needs numpy
             True
             sage: 12 == numpy.int8('12')                                                # needs numpy

src/sage/rings/polynomial/polynomial_element.pyx

@@ -8851,7 +8851,7 @@ cdef class Polynomial(CommutativePolynomial):
                          'you expect.')
 
                 import numpy
-                from numpy.linalg.linalg import LinAlgError
+                from numpy.linalg import LinAlgError
                 from sage.rings.complex_double import CDF
 
                 numpy_dtype = ('complex' if input_complex else 'double')

src/sage/rings/real_mpfi.pyx

@@ -231,6 +231,8 @@ specified if given a non-interval and an interval::
 TESTS::
 
     sage: import numpy                                                                  # needs numpy
+    sage: if int(numpy.version.short_version[0]) > 1:                                   # needs numpy
+    ....:     numpy.set_printoptions(legacy="1.25")                                     # needs numpy
     sage: RIF(2) == numpy.int8('2')                                                     # needs numpy
     True
     sage: numpy.int8('2') == RIF(2)                                                     # needs numpy