Fixed case of degeneracy in converting a ConvexHullDomain to LinQuadD…

…omain

psdr/domains/convexhull.py

@@ -22,6 +22,49 @@
 from ..geometry import unique_points
 from .euclidean import TOL
 
+
+
+
+def _hull_to_linineq(X):
+	r"""
+	"""
+	if X.shape[0] == 1:
+		A = None
+		b = None
+		vertices = X
+		A_eq = np.eye(X.shape[1])
+		b_eq = X.flatten()
+		return A, b, A_eq, b_eq, vertices
+
+	# Find a low-dimensional subspace on which this data lives
+	Xc = np.mean(X, axis = 0)
+	Xdiff = (X.T - Xc.reshape(-1,1)).T
+	U, s, VT = scipy.linalg.svd(Xdiff, full_matrices = False)
+	I = np.isclose(s, 0)
+	Q, _ = scipy.linalg.qr(VT[~I].T, mode = 'full')
+	dim = np.sum(~I)
+	# component where coordinates change
+	Qpar = Q[:,0:dim]
+	# orthogonal complement
+	Qperp = Q[:, dim:]
+
+	Y = X @ Qpar
+	if Y.shape[1] > 1:
+		hull = ConvexHull(Y)
+		A = hull.equations[:,:-1] @ Qpar.T
+		b = -hull.equations[:,-1]
+		vertices = X[hull.vertices]
+	else:
+		# this is an interval which is not handled by Qhull
+		A = np.array([1, -1]).reshape(2,1) @ Qpar.T
+		b = np.array([np.max(Y), -np.min(Y)])
+		vertices = np.array([X[np.argmax(Y)], X[np.argmin(Y)]])
+
+	A_eq = Qperp.T
+	b_eq = np.mean( X @ Qperp, axis = 0)
+	return A, b, A_eq, b_eq, vertices
+
+
 class ConvexHullDomain(LinQuadDomain):
 	r"""Define a domain that is the interior of a convex hull of points.
 
@@ -114,21 +157,30 @@ def to_linineq(self, **kwargs):
 		r""" Convert the domain into a LinIneqDomain
 
 		"""
-		if len(self) > 1:
-			hull = ConvexHull(self._X) 
-			A = hull.equations[:,:-1]
-			b = -hull.equations[:,-1]
-			dom_hull = LinQuadDomain(A = A, b = b, names = self.names, **kwargs)
-			dom_hull.vertices = np.copy(self._X[hull.vertices])
-			dom = dom_hull.add_constraints(A = self.A, b = self.b, A_eq = self.A_eq, b_eq = self.b_eq,
+		A, b, A_eq, b_eq, vertices = _hull_to_linineq(self._X)
+		dom_hull = LinQuadDomain(A = A, b = b, A_eq = A_eq, b_eq = b_eq, names = self.names, **kwargs)
+		dom_hull.vertices = np.copy(vertices)
+		# Add back in non-hull constraints
+		dom = dom_hull.add_constraints(A = self.A, b = self.b, A_eq = self.A_eq, b_eq = self.b_eq,
 				Ls = self.Ls, ys = self.ys, rhos = self.rhos)
-		else:
-			lb = self.corner([-1])
-			ub = self.corner([1])
-			dom = BoxDomain(lb, ub, names = self.names, **kwargs)
-
 		return dom
 
+#		if len(self) > 1:
+#			hull = ConvexHull(self._X) 
+#			vertices = self._X[hull.vertices]
+#			A = hull.equations[:,:-1]
+#			b = -hull.equations[:,-1]
+#			dom_hull = LinQuadDomain(A = A, b = b, names = self.names, **kwargs)
+#			dom_hull.vertices = np.copy(self._X[hull.vertices])
+#			dom = dom_hull.add_constraints(A = self.A, b = self.b, A_eq = self.A_eq, b_eq = self.b_eq,
+#				Ls = self.Ls, ys = self.ys, rhos = self.rhos)
+#		else:
+#			lb = self.corner([-1])
+#			ub = self.corner([1])
+#			dom = BoxDomain(lb, ub, names = self.names, **kwargs)
+#
+#		return dom
+
 	def coefficients(self, x, **kwargs):
 		r""" Find the coefficients of the convex combination of elements in the space yielding x
 

psdr/domains/euclidean.py

@@ -86,12 +86,17 @@ def is_empty(self):
 			return False
 		except EmptyDomainException:
 			return True
+		except UnboundedDomainException:
+			return False
 
 
 	@cached_property
 	def is_point(self):
 		try:
-			U = ortho_group.rvs(len(self))
+			if len(self) > 1:
+				U = ortho_group.rvs(len(self))
+			else:
+				U = np.ones((1,1))
 
 			for u in U:
 				x1 = self.corner(u)
@@ -100,7 +105,6 @@ def is_point(self):
 					return False
 
 			return True
-
 		except EmptyDomainException:
 			return False
 		except UnboundedDomainException:
@@ -284,16 +288,16 @@ def _corner(self, p, **kwargs):
 
 		problem.solve(**kwargs)
 		if problem.status in ['infeasible']:
-			self._empty = True
-			self._unbounded = False
-			self._point = False
+			#self._empty = True
+			#self._unbounded = False
+			#self._point = False
 			raise EmptyDomainException	
 		# For some reason, if no constraints are provided CVXPY doesn't note
 		# the domain is unbounded
 		elif problem.status in ['unbounded'] or len(constraints) == 0:
-			self._unbounded = True
-			self._empty = False
-			self._point = False
+			#self._unbounded = True
+			#self._empty = False
+			#self._point = False
 			raise UnboundedDomainException
 		elif problem.status not in ['optimal', 'optimal_inaccurate']:
 			raise SolverError("CVXPY exited with status '%s'" % problem.status)
@@ -675,7 +679,7 @@ def random_direction(self, x):
 			# Orthogonalize against equality constarints constraints
 			p = p - Qeq @ (Qeq.T @ p)
 			# check that the direction isn't hitting a constraint
-			if self.extent(x, p) > 0:
+			if x is None or self.extent(x, p) > 0:
 				break
 		return p	
 
@@ -825,7 +829,10 @@ def _corner_center(self):
 		# at the center of the domain.
 
 		# Generate random orthogonal directions to sample
-		U = ortho_group.rvs(len(self))
+		if len(self) > 1:
+			U = ortho_group.rvs(len(self))
+		else:
+			U = np.ones((1,1))
 		X = []
 		for i in range(len(self)):
 			X += [self.corner(U[:,i])]
@@ -856,7 +863,7 @@ def _hit_and_run(self, _recurse = 2):
 			# Get the current location of where the hit and run sampler is
 			x0 = self._hit_and_run_state
 			if x0 is None: raise AttributeError
-		except (AttributeError,):
+		except AttributeError:
 
 			try:
 				# The simpliest and inexpensive approach is to start hit and run
@@ -866,7 +873,7 @@ def _hit_and_run(self, _recurse = 2):
 
 				# TODO: chebyshev_center breaks when running test_lipschitz_sample yielding a SolverError
 				# It really shouldn't error
-			except (NotImplementedError):
+			except (AttributeError, NotImplementedError):
 				x0 = self._corner_center()
 			self._hit_and_run_state = x0
 

tests/domains/test_convexhull.py

@@ -4,6 +4,26 @@
 from scipy.linalg import orth
 import pytest
 
+from psdr.domains.convexhull import _hull_to_linineq
+
+
+@pytest.mark.parametrize("M", [1,2,3,4,5,6,7,8,9,10])
+@pytest.mark.parametrize("m", [1,2,3,4])
+def test_hull_to_ineq(M, m):
+	np.random.seed(0)
+	X = np.random.randn(M,m)
+	dom_hull = psdr.ConvexHullDomain(X)
+	dom_linineq = dom_hull.to_linineq()
+
+	for it in range(10):
+		if M > 1:
+			# Hit and run will error out with a point domain
+			x = dom_hull._hit_and_run()
+		else:
+			x = dom_hull.sample()
+		assert dom_linineq.isinside(x)
+
+
 def test_sample():
 	np.random.seed(0)
 	m = 4
@@ -107,4 +127,5 @@ def test_A_eq_deficient(m, nullspace_dim):
 if __name__ == '__main__':
 	#test_sample()
 	#test_A_eq(4, 2)
-	test_A_eq_deficient(4, 0)
+	#test_A_eq_deficient(4, 0)
+	test_hull_to_ineq(1, 4)

tests/domains/test_euclidean.py

@@ -28,67 +28,58 @@ def test_empty(m = 3):
 	# Unbounded
 	dom = psdr.LinQuadDomain(lb = -np.inf*np.ones(m), ub = np.inf*np.ones(m))
 	assert dom.is_empty == False
-	assert dom._empty == False
-	assert dom._point == False
-	assert dom._unbounded == True
+	assert dom.is_point == False
+	assert dom.is_unbounded == True
 
 	# Emtpy
 	dom = psdr.LinQuadDomain(lb = 1*np.ones(m), ub = -1*np.ones(m))
 	assert dom.is_empty == True
-	assert dom._empty == True
-	assert dom._point == False
-	assert dom._unbounded == False
+	assert dom.is_point == False
+	assert dom.is_unbounded == False
 
 def test_point(m = 3):
 	# Unbounded
 	dom = psdr.LinQuadDomain(lb = -np.inf*np.ones(m), ub = np.inf*np.ones(m))
 	assert dom.is_point == False
-	assert dom._empty == False
-	assert dom._point == False
-	assert dom._unbounded == True
+	assert dom.is_empty == False
+	assert dom.is_unbounded == True
 
 	# Emtpy
 	dom = psdr.LinQuadDomain(lb = 1*np.ones(m), ub = -1*np.ones(m))
 	assert dom.is_point == False
-	assert dom._empty == True
-	assert dom._point == False
-	assert dom._unbounded == False
+	assert dom.is_empty == True
+	assert dom.is_unbounded == False
 
 	# Point
 	dom = psdr.LinQuadDomain(lb = 1*np.ones(m), ub = 1*np.ones(m))
 	assert dom.is_point == True
-	assert dom._empty == False
-	assert dom._point == True
-	assert dom._unbounded == False
+	assert dom.is_empty == False
+	assert dom.is_unbounded == False
 
 def test_unbounded(m = 3):
 	# Unbounded
 	dom = psdr.LinQuadDomain(lb = -np.inf*np.ones(m), ub = np.inf*np.ones(m))
 	assert dom.is_unbounded == True
-	assert dom._empty == False
-	assert dom._point == False
-	assert dom._unbounded == True
+	assert dom.is_empty == False
+	assert dom.is_point == False
 
 	# Empty	
 	dom = psdr.LinQuadDomain(lb = 1*np.ones(m), ub = -1*np.ones(m))
 	assert dom.is_unbounded == False
-	assert dom._empty == True
-	assert dom._point == False
-	assert dom._unbounded == False
+	assert dom.is_empty == True
+	assert dom.is_point == False
 
 	# Point
 	dom = psdr.LinQuadDomain(lb = 1*np.ones(m), ub = 1*np.ones(m))
 	assert dom.is_unbounded == False
-	assert dom._empty == False
-	assert dom._point == True
-	assert dom._unbounded == False
+	assert dom.is_empty == False
+	assert dom.is_point == True
 
 	# Box
 	dom = psdr.LinQuadDomain(lb = -1*np.ones(m), ub = 1*np.ones(m))
 	assert dom.is_unbounded == False
-	assert dom._empty == False
-	assert dom._point == False
-	assert dom._unbounded == False
+	assert dom.is_empty == False
+	assert dom.is_point == False
 
 
 def test_sweep(m = 5):

tests/sample/test_stretch.py

@@ -82,6 +82,6 @@ def test_stretch_sample_uniform(m,n):
 			assert pvalue > 1e-3
 
 if __name__ == '__main__':
-	test_stretch_sample(5, [3,4, 5])
+	test_stretch_sample(5, [3, 3])
 	#test_stretch_sample_uniform(5, 1)
 	#test_stretch_sample_uniform(5, 2)