Resolve PyNumero incompatibilities with NumPy2

.github/workflows/test_pr_and_main.yml

@@ -141,6 +141,15 @@ jobs:
           TARGET: linux
           PYENV: pip
 
+        - os: ubuntu-latest
+          python: 3.12
+          other: /numpy2
+          slim: 1
+          skip_doctest: 1
+          TARGET: linux
+          PYENV: pip
+          PACKAGES: "gurobipy dill numpy>2.0 scipy networkx"
+
         - os: ubuntu-latest
           python: 3.9
           other: /pyutilib

pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_block.py

@@ -15,10 +15,10 @@
 from pyomo.core.expr.visitor import identify_variables
 import pyomo.environ as pyo
 
+from pyomo.common.dependencies import networkx_available as nx_available
 from pyomo.contrib.pynumero.dependencies import (
     numpy as np,
     numpy_available,
-    scipy,
     scipy_available,
 )
 
@@ -151,6 +151,7 @@ def flow_out_eqn(m, t):
 
 
 class TestExternalGreyBoxBlock(unittest.TestCase):
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_construct_scalar(self):
         m = pyo.ConcreteModel()
         m.ex_block = ExternalGreyBoxBlock(concrete=True)
@@ -171,6 +172,7 @@ def test_construct_scalar(self):
         self.assertEqual(len(block.outputs), 0)
         self.assertEqual(len(block._equality_constraint_names), 2)
 
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_construct_indexed(self):
         block = ExternalGreyBoxBlock([0, 1, 2], concrete=True)
         self.assertIs(type(block), IndexedExternalGreyBoxBlock)
@@ -192,6 +194,7 @@ def test_construct_indexed(self):
             self.assertEqual(len(b._equality_constraint_names), 2)
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_solve_square(self):
         m = pyo.ConcreteModel()
         m.ex_block = ExternalGreyBoxBlock(concrete=True)
@@ -234,6 +237,7 @@ def test_solve_square(self):
         self.assertAlmostEqual(m_ex.y.value, y.value, delta=1e-8)
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_optimize(self):
         m = pyo.ConcreteModel()
         m.ex_block = ExternalGreyBoxBlock(concrete=True)
@@ -292,6 +296,7 @@ def test_optimize(self):
         self.assertAlmostEqual(m_ex.y.value, y.value, delta=1e-8)
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_optimize_with_cyipopt_for_inner_problem(self):
         # Use CyIpopt, rather than the default SciPy solvers,
         # for the inner problem
@@ -427,6 +432,7 @@ def test_optimize_no_decomposition(self):
         self.assertAlmostEqual(m_ex.x.value, x.value, delta=1e-8)
         self.assertAlmostEqual(m_ex.y.value, y.value, delta=1e-8)
 
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_construct_dynamic(self):
         m = make_dynamic_model()
         time = m.time
@@ -504,6 +510,7 @@ def test_construct_dynamic(self):
         )
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_solve_square_dynamic(self):
         # Create the "external model"
         m = make_dynamic_model()
@@ -571,6 +578,7 @@ def linking_constraint_rule(m, i, t):
             self.assertStructuredAlmostEqual(values, target_values, delta=1e-5)
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_optimize_dynamic(self):
         # Create the "external model"
         m = make_dynamic_model()
@@ -653,6 +661,7 @@ def linking_constraint_rule(m, i, t):
             self.assertStructuredAlmostEqual(values, target_values, delta=1e-5)
 
     @unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+    @unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
     def test_optimize_dynamic_references(self):
         """
         When when pre-existing variables are attached to the EGBB
@@ -717,7 +726,8 @@ def test_optimize_dynamic_references(self):
             self.assertStructuredAlmostEqual(values, target_values, delta=1e-5)
 
 
-class TestPyomoNLPWithGreyBoxBLocks(unittest.TestCase):
+@unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
+class TestPyomoNLPWithGreyBoxBlocks(unittest.TestCase):
     def test_set_and_evaluate(self):
         m = pyo.ConcreteModel()
         m.ex_block = ExternalGreyBoxBlock(concrete=True)

pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_model.py

@@ -13,6 +13,7 @@
 import pyomo.common.unittest as unittest
 import pyomo.environ as pyo
 
+from pyomo.common.dependencies import networkx_available as nx_available
 from pyomo.contrib.pynumero.dependencies import (
     numpy as np,
     numpy_available,
@@ -513,6 +514,7 @@ def test_explicit_zeros(self):
         np.testing.assert_allclose(hess.data, data, rtol=1e-8)
 
 
+@unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
 class TestExternalPyomoModel(unittest.TestCase):
     def test_evaluate_SimpleModel1(self):
         model = SimpleModel1()
@@ -838,6 +840,7 @@ def test_evaluate_hessian_lagrangian_SimpleModel2x2_1(self):
             np.testing.assert_allclose(hess_lag, expected_hess_lag, rtol=1e-8)
 
 
+@unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
 class TestUpdatedHessianCalculationMethods(unittest.TestCase):
     """
     These tests exercise the methods for fast Hessian-of-Lagrangian
@@ -1021,6 +1024,7 @@ def test_evaluate_hessian_equality_constraints_order(self):
             )
 
 
+@unittest.skipUnless(nx_available, "SCCImplicitFunctionSolver requires networkx")
 class TestScaling(unittest.TestCase):
     def con_3_body(self, x, y, u, v):
         return 1e5 * x**2 + 1e4 * y**2 + 1e1 * u**2 + 1e0 * v**2

pyomo/contrib/pynumero/sparse/base_block.py

@@ -11,6 +11,8 @@
 
 # These classes are for checking types consistently and raising errors
 
+from ..dependencies import numpy as np
+
 
 class BaseBlockVector(object):
     """Base class for block vectors"""
@@ -177,3 +179,128 @@ def transpose(self, *axes):
     def tostring(self, order='C'):
         msg = "tostring not implemented for {}".format(self.__class__.__name__)
         raise NotImplementedError(msg)
+
+
+#: NumPy ufuncs that take one vector and are compatible with pyNumero vectors
+vec_unary_ufuncs = {
+    ## MATH ufuncs
+    np.negative,
+    np.positive,
+    np.absolute,
+    np.fabs,
+    np.rint,
+    np.sign,
+    np.conj,
+    np.conjugate,
+    np.exp,
+    np.exp2,
+    np.log,
+    np.log2,
+    np.log10,
+    np.expm1,
+    np.log1p,
+    np.sqrt,
+    np.square,
+    np.cbrt,
+    np.reciprocal,
+    ## TRIG ufuncs
+    np.sin,
+    np.cos,
+    np.tan,
+    np.arcsin,
+    np.arccos,
+    np.arctan,
+    np.sinh,
+    np.cosh,
+    np.tanh,
+    np.arcsinh,
+    np.arccosh,
+    np.arctanh,
+    np.degrees,
+    np.radians,
+    np.deg2rad,
+    np.rad2deg,
+    ## COMPARISON ufuncs
+    np.logical_not,
+    ## BIT-TWIDDLING ufuncs
+    np.invert,
+    ## FLOATING ufuncs
+    np.isfinite,
+    np.isinf,
+    np.isnan,
+    # np.isnat,  # only defined for datetime
+    np.fabs,  # numpy docs list here and in MATH
+    np.signbit,
+    np.spacing,
+    # np.modf, # disabled because shape is not preserved
+    # np.frexp, # disabled because shape is not preserved
+    np.floor,
+    np.ceil,
+    np.trunc,
+    # OTHER (not listed in ufuncs docs)
+    np.abs,
+}
+
+#: NumPy ufuncs that take two vectors and are compatible with pyNumero vectors
+vec_binary_ufuncs = {
+    ## MATH ufuncs
+    np.add,
+    np.subtract,
+    np.multiply,
+    # np.matmult, # disabled because shape is not preserved
+    np.divide,
+    np.logaddexp,
+    np.logaddexp2,
+    np.true_divide,
+    np.floor_divide,
+    np.power,
+    np.float_power,
+    np.remainder,
+    np.mod,
+    np.fmod,
+    # np.divmod,  # dieabled because shape is not preserved
+    np.heaviside,
+    np.gcd,
+    np.lcm,
+    ## TRIG ufuncs
+    np.arctan2,
+    np.hypot,
+    ## BIT-TWIDDLING ufuncs
+    np.bitwise_and,
+    np.bitwise_or,
+    np.bitwise_xor,
+    np.left_shift,
+    np.right_shift,
+    ## COMPARISON ufuncs
+    np.greater,
+    np.greater_equal,
+    np.less,
+    np.less_equal,
+    np.not_equal,
+    np.equal,
+    np.logical_and,
+    np.logical_or,
+    np.logical_xor,
+    np.maximum,
+    np.minimum,
+    np.fmax,
+    np.fmin,
+    ## FLOATING ufincs
+    np.copysign,
+    np.nextafter,
+    np.ldexp,
+    np.fmod,  # numpy docs list here and in MATH
+}
+
+#: NumPy ufuncs can be used as reductions for pyNumero vectors
+vec_associative_reductions = {
+    np.add,
+    np.multiply,
+    np.bitwise_and,
+    np.bitwise_or,
+    np.bitwise_xor,
+    np.maximum,
+    np.minimum,
+    np.fmax,
+    np.fmin,
+}