diff --git a/pyomo/contrib/piecewise/piecewise_linear_function.py b/pyomo/contrib/piecewise/piecewise_linear_function.py index 666beeddbd6..f0a92ef139f 100644 --- a/pyomo/contrib/piecewise/piecewise_linear_function.py +++ b/pyomo/contrib/piecewise/piecewise_linear_function.py @@ -9,6 +9,9 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +import logging + +from pyomo.common import DeveloperError from pyomo.common.autoslots import AutoSlots from pyomo.common.collections import ComponentMap from pyomo.common.dependencies import numpy as np @@ -30,6 +33,8 @@ # be zero. ZERO_TOLERANCE = 1e-8 +logger = logging.getLogger(__name__) + class PiecewiseLinearFunctionData(_BlockData): _Block_reserved_words = Any @@ -282,11 +287,34 @@ def _construct_from_function_and_points(self, obj, parent, nonlinear_function): logger.error("Unable to triangulate the set of input points.") raise - obj._points = [pt for pt in points] - obj._simplices = [simplex for simplex in map(tuple, triangulation.simplices)] + # Get the points for the triangulation because they might not all be + # there if any were coplanar. + obj._points = [pt for pt in map(tuple, triangulation.points)] + obj._simplices = [] + for simplex in triangulation.simplices: + # Check whether or not the simplex is degenerate. If it is, we will + # just drop it. + points = triangulation.points[simplex].transpose() + if ( + np.linalg.det( + points[:, 1:] + - np.append(points[:, : dimension - 1], points[:, [0]], axis=1) + ) + != 0 + ): + obj._simplices.append(tuple(sorted(simplex))) + + # It's possible that qhull dropped some points if there were numerical + # issues with them (e.g., if they were redundant). We'll be polite and + # tell the user: + for pt in triangulation.coplanar: + logger.info( + "The Delaunay triangulation dropped the point with index " + "%s from the triangulation." % pt[0] + ) return self._construct_from_function_and_simplices( - obj, parent, nonlinear_function + obj, parent, nonlinear_function, simplices_are_user_defined=False ) def _construct_from_univariate_function_and_segments(self, obj, func): @@ -301,7 +329,9 @@ def _construct_from_univariate_function_and_segments(self, obj, func): return obj @_define_handler(_handlers, True, False, True, False) - def _construct_from_function_and_simplices(self, obj, parent, nonlinear_function): + def _construct_from_function_and_simplices( + self, obj, parent, nonlinear_function, simplices_are_user_defined=True + ): if obj._simplices is None: obj._get_simplices_from_arg(self._simplices_rule(parent, obj._index)) simplices = obj._simplices @@ -335,12 +365,31 @@ def _construct_from_function_and_simplices(self, obj, parent, nonlinear_function A[i, j + 1] = nonlinear_function(*pt) A[i + 1, :] = 0 A[i + 1, dimension] = -1 - # This system has a solution unless there's a bug--we know there is - # a hyperplane that passes through dimension + 1 points (and the - # last equation scales it so that the coefficient for the output - # of the nonlinear function dimension is -1, so we can just read - # off the linear equation in the x space). - normal = np.linalg.solve(A, b) + # This system has a solution unless there's a bug--we filtered the + # simplices to make sure they are full-dimensional, so we know there + # is a hyperplane that passes through these dimension + 1 points (and the + # last equation scales it so that the coefficient for the output of + # the nonlinear function dimension is -1, so we can just read off + # the linear equation in the x space). + try: + normal = np.linalg.solve(A, b) + except np.linalg.LinAlgError as e: + logger.warning('LinAlgError: %s' % e) + msg = ( + "When calculating the hyperplane approximation over the simplex " + "with index %s, the matrix was unexpectedly singular. This " + "likely means that this simplex is degenerate" % num_piece + ) + + if simplices_are_user_defined: + raise ValueError(msg) + # otherwise it's our fault, and I was hoping this is unreachable + # code... + raise DeveloperError( + msg + + "and that it should have been filtered out of the triangulation" + ) + obj._linear_functions.append(_multivariate_linear_functor(normal)) return obj diff --git a/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py b/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py index abb5bec0c9f..5c91381d113 100644 --- a/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py +++ b/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py @@ -10,9 +10,11 @@ # ___________________________________________________________________________ from io import StringIO +import logging import pickle from pyomo.common.dependencies import attempt_import +from pyomo.common.log import LoggingIntercept import pyomo.common.unittest as unittest from pyomo.contrib.piecewise import PiecewiseLinearFunction from pyomo.core.expr.compare import ( @@ -314,7 +316,7 @@ def g2(x1, x2): self.assertEqual(str(m.c.body.expr), "pw(x1, x2)") self.assertIs(m.c.body.expr.pw_linear_function, m.pw) - @unittest.skipUnless(numpy_available, "numpy are not available") + @unittest.skipUnless(numpy_available, "numpy is not available") def test_evaluate_pw_linear_function(self): # NOTE: This test requires numpy because it is used to check which # simplex a point is in @@ -337,3 +339,124 @@ def g2(x1, x2): self.assertAlmostEqual(m.pw(1, 3), g1(1, 3)) self.assertAlmostEqual(m.pw(2.5, 6), g2(2.5, 6)) self.assertAlmostEqual(m.pw(0.2, 4.3), g2(0.2, 4.3)) + + +class TestTriangulationProducesDegenerateSimplices(unittest.TestCase): + cube_extreme_pt_indices = [ + {10, 11, 13, 14, 19, 20, 22, 23}, # right bottom back + {9, 10, 12, 13, 18, 19, 21, 22}, # right bottom front + {0, 1, 3, 4, 9, 10, 12, 13}, # left bottom front + {1, 2, 4, 5, 10, 11, 13, 14}, # left bottom back + {3, 4, 6, 7, 12, 13, 15, 16}, # left top front + {4, 5, 7, 8, 13, 14, 16, 17}, # left top back + {12, 13, 15, 16, 21, 22, 24, 25}, # right top front + {13, 14, 16, 17, 22, 23, 25, 26}, # right top back + ] + + def make_model(self): + m = ConcreteModel() + + m.f = lambda x1, x2, y: x1 * x2 + y + # This is a 2x2 stack of cubes, so there are 8 total cubes, each of which + # will get divided into 6 simplices. + m.points = [ + (-2.0, 0.0, 1.0), + (-2.0, 0.0, 4.0), + (-2.0, 0.0, 7.0), + (-2.0, 1.5, 1.0), + (-2.0, 1.5, 4.0), + (-2.0, 1.5, 7.0), + (-2.0, 3.0, 1.0), + (-2.0, 3.0, 4.0), + (-2.0, 3.0, 7.0), + (-1.5, 0.0, 1.0), + (-1.5, 0.0, 4.0), + (-1.5, 0.0, 7.0), + (-1.5, 1.5, 1.0), + (-1.5, 1.5, 4.0), + (-1.5, 1.5, 7.0), + (-1.5, 3.0, 1.0), + (-1.5, 3.0, 4.0), + (-1.5, 3.0, 7.0), + (-1.0, 0.0, 1.0), + (-1.0, 0.0, 4.0), + (-1.0, 0.0, 7.0), + (-1.0, 1.5, 1.0), + (-1.0, 1.5, 4.0), + (-1.0, 1.5, 7.0), + (-1.0, 3.0, 1.0), + (-1.0, 3.0, 4.0), + (-1.0, 3.0, 7.0), + ] + return m + + @unittest.skipUnless( + scipy_available and numpy_available, "scipy and/or numpy are not available" + ) + def test_degenerate_simplices_filtered(self): + m = self.make_model() + pw = m.approx = PiecewiseLinearFunction(points=m.points, function=m.f) + + # check that all the points got used + self.assertEqual(len(pw._points), 27) + for p_model, p_pw in zip(m.points, pw._points): + self.assertEqual(p_model, p_pw) + + # Started with a 2x2 grid of cubes, and each is divided into 6 + # simplices. It's crazy degenerate in terms of *how* this is done, but + # that's the point of this test. + self.assertEqual(len(pw._simplices), 48) + simplex_in_cube = {idx: 0 for idx in range(8)} + for simplex in pw._simplices: + for i, vertex_set in enumerate(self.cube_extreme_pt_indices): + if set(simplex).issubset(vertex_set): + simplex_in_cube[i] += 1 + # verify the simplex is full-dimensional + pts = np.array([pw._points[j] for j in simplex]).transpose() + A = pts[:, 1:] - np.append(pts[:, :2], pts[:, [0]], axis=1) + self.assertNotEqual(np.linalg.det(A), 0) + + # Check that they are 6 to a cube, as expected + for num in simplex_in_cube.values(): + self.assertEqual(num, 6) + + @unittest.skipUnless( + scipy_available and numpy_available, "scipy and/or numpy are not available" + ) + def test_redundant_points_logged(self): + m = self.make_model() + # add a redundant point + m.points.append((-2, 0, 1)) + + out = StringIO() + with LoggingIntercept( + out, 'pyomo.contrib.piecewise.piecewise_linear_function', level=logging.INFO + ): + m.approx = PiecewiseLinearFunction(points=m.points, function=m.f) + + self.assertIn( + "The Delaunay triangulation dropped the point with index 27 " + "from the triangulation", + out.getvalue(), + ) + + @unittest.skipUnless(numpy_available, "numpy is not available") + def test_user_given_degenerate_simplex_error(self): + m = self.make_model() + with self.assertRaisesRegex( + ValueError, + "When calculating the hyperplane approximation over the simplex " + "with index 0, the matrix was unexpectedly singular. This " + "likely means that this simplex is degenerate", + ): + m.pw = PiecewiseLinearFunction( + simplices=[ + ( + (-2.0, 0.0, 1.0), + (-2.0, 0.0, 4.0), + (-2.0, 1.5, 1.0), + (-2.0, 1.5, 4.0), + ) + ], + function=m.f, + )