Preliminary stretched sampling code

psdr/sample/stretch.py

@@ -2,12 +2,16 @@
 """
 import numpy as np
 from .initial import initial_sample
+from .maximin_coffeehouse import maximin_coffeehouse 
 from ..geometry import voronoi_vertex_sample, cdist
+from cvxpy.error import SolverError
 
 def _maximin_block_candidate(domain, Xf, L):
 	r"""
 	"""
 	X0 = initial_sample(domain, L, Nsamp = 2*len(Xf))
+	#X0 = maximin_coffeehouse(domain, 2*len(Xf), L = L, N0 = 1, full = False, verbose = False)
+	#X0 = domain.sample(2*len(Xf))
 	X = voronoi_vertex_sample(domain, Xf, X0, L = L, randomize = False) 
 	# find the point furthest away 
 	D = cdist(Xf, X, L)
@@ -17,7 +21,7 @@ def _maximin_block_candidate(domain, Xf, L):
 
 
 
-def stretch_sample_domain(domain, X, Ls):
+def stretch_sample_domain(domain, X, Ls, verbose = False):
 	r"""
 
 
@@ -34,7 +38,8 @@ def stretch_sample_domain(domain, X, Ls):
 	Ls = Ls.copy()
 	# A copy of the domain we'll keep adding constraints to	
 	domain = domain.copy()
-
+
+	it = 0	
 	while len(Ls) > 0:
 		# Compute candidate points and distances
 		xd = [list(_maximin_block_candidate(domain, X, L)) + [L, k] for k, L in enumerate(Ls)]
@@ -44,19 +49,32 @@ def stretch_sample_domain(domain, X, Ls):
 		# Remove this L from the list
 		Ls.pop(k)
 
+
 		# Try to add the new constraint
 		domain_new = domain.add_constraints(A_eq = L, b_eq = L @ x)
+
 
 		# If we have an empty domain, stop and return the existing domain
-		# otherwise continue	
-		if domain_new.is_empty: 	break
-		else: domain_new = domain
+		# otherwise continue
+		try:
+			if domain_new.is_empty: 	break
+			else: domain = domain_new
+		except SolverError:
+			break
 
+		it += 1
+		if verbose:
+			print(f'iter {it:2d}, domain size {len(domain)}, eq constraints {len(domain.A_eq)}')
+
+		# If we have more constraints than the dimension of the space, stop
+		if len(domain.A_eq) > len(domain): break
 
 	return domain
 
-def stretch_sample(domain, Xf, Ls):
-	domain_stretch = stretch_sample_domain(domain, Xf, Ls)
+def stretch_sample(domain, Xf, Ls, verbose = False):
+	domain_stretch = stretch_sample_domain(domain, Xf, Ls, verbose = verbose)
+	#print(domain_stretch.A_eq)
+	#print(domain_stretch.b_eq)
 	X0 = domain_stretch.sample(100)
 	X = voronoi_vertex_sample(domain_stretch, Xf, X0)
 	D = cdist(Xf, X)
@@ -65,4 +83,8 @@ def stretch_sample(domain, Xf, Ls):
 	return X[k]
 
 
-
+
+class StretchSample:
+	def __init__(self, fun, X = None, fX = None):
+		self.fun = fun
+

tests/sample/test_initial.py

@@ -0,0 +1,16 @@
+import numpy as np
+import psdr
+
+
+def test_initial():
+	np.random.seed(0)
+	m = 5
+	L = np.random.randn(1,m)
+	domain = psdr.BoxDomain(-np.ones(m), np.ones(m))
+
+	X0 = psdr.initial_sample(domain, L, 10)
+	print(X0)	
+
+
+if __name__ == '__main__':
+	test_initial()

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()
-	test_stretch_sample_uniform(5, 1)
-	test_stretch_sample_uniform(5, 2)
+	test_stretch_sample(5, [3,4, 5])
+	#test_stretch_sample_uniform(5, 1)
+	#test_stretch_sample_uniform(5, 2)