Merge branch 'master' of github.com:jeffrey-hokanson/PSDR

Reproducibility/Lipschitz/fig_design.py

@@ -2,6 +2,7 @@
 import psdr, psdr.demos
 from psdr.pgf import PGF
 import seaborn as sns
+import matplotlib.pyplot as plt
 
 funs = [psdr.demos.OTLCircuit(), 
 		psdr.demos.Borehole(),
@@ -13,10 +14,17 @@
 		lambda dom, M, L: psdr.minimax_lloyd(dom, M, L = L)	
 		]
 alg_names = ['random', 'LHS', 'minimax']
+
+
 # Number of repetitions
 Ms = [100,100, 10]
 #Ms = [1,1,1]
 
+algs = algs[2:]
+alg_names = alg_names[2:]
+Ms = Ms[2:]
+
+funs = funs[1:2]
 
 Nsamp = 20
 
@@ -34,24 +42,25 @@
 
 	L = lip.L
 
-	# Samples to use when estimating dispersion
-	X0 = psdr.maximin_coffeehouse(fun.domain, 5000, L = L, N0 = 10)
-
+	#X0 = psdr.maximin_coffeehouse(fun.domain, 5000, L = L, N0 = 50)
+	X0 = np.vstack([psdr.random_sample(fun.domain, 5000), fun.domain.sample_grid(2)])
+
+	plt.clf()	
+
 	# Now perform designs
 	for alg, alg_name, M in zip(algs, alg_names, Ms):
 		dispersion = []
 		# We can get a very sloppy fit with more points
 		for i in range(M):
 			np.random.seed(i)
 			X = alg(fun.domain, Nsamp, L)
-			dist = psdr.fill_distance_estimate(fun.domain, X, L = L, X0 = X0) 
+			dist = psdr.fill_distance_estimate(fun.domain, X, L = L, X0 = np.copy(X0)) 
 			dispersion.append(dist)
 			print(f'{alg_name:20s} : {i:4d} dispersion {dist:10.5e}')
 
 		ax = sns.swarmplot(dispersion)
-		x, y = np.array(ax.collections[0].get_offsets()).T
+		x, y = np.array(ax.collections[-1].get_offsets()).T
 		pgf = PGF()
 		pgf.add('x', x)
 		pgf.add('y', y)
 		pgf.write(f'data/fig_design_{name}_{alg_name}.dat')	
-

Reproducibility/Lipschitz/fig_epsilon_rank.py

@@ -16,6 +16,7 @@
 fun = psdr.demos.OTLCircuit()
 
 np.random.seed(0)
+
 @memory.cache
 def generate_X(fun, M):
 	np.random.seed(0)
@@ -28,7 +29,20 @@ def generate_X(fun, M):
 	X = psdr.minimax_lloyd(fun.domain, M, L = L, verbose = True)
 	return X
 
-X = generate_X(fun, 200)
+@memory.cache
+def generate_X_maximin(fun, M):
+	np.random.seed(0)
+	X = fun.domain.sample(1000)
+	grads = fun.grad(X)
+	lip = psdr.LipschitzMatrix(verbose = True, abstol = 1e-7, reltol = 1e-7, feastol = 1e-7)
+	lip.fit(grads = grads)
+	L = lip.L
+
+	X0 = psdr.maximin_coffeehouse(fun.domain, M, L = L, verbose =True)		
+	X = psdr.maximin_block(fun.domain, M, L = L, verbose = True, Xhat = X0)
+	return X
+
+X = generate_X_maximin(fun, 200)
 fX = fun(X)
 
 print(len(fX))
@@ -50,6 +64,8 @@ def lipschitz(X, fX, epsilon):
 lam2 = np.zeros(epsilon.shape)
 lam3 = np.zeros(epsilon.shape)
 lam4 = np.zeros(epsilon.shape)
+lam5 = np.zeros(epsilon.shape)
+lam6 = np.zeros(epsilon.shape)
 
 
 for k, eps in enumerate(epsilon):
@@ -62,6 +78,8 @@ def lipschitz(X, fX, epsilon):
 	lam2[k] = ew[1]
 	lam3[k] = ew[2]
 	lam4[k] = ew[3]
+	lam5[k] = ew[4]
+	lam6[k] = ew[5]
 	print(f"=====> epsilon {eps:8.2e}, rank {rank[k]:2d}, obj {obj[k]:10.5e}") 
 
 pgf = PGF()
@@ -72,4 +90,6 @@ def lipschitz(X, fX, epsilon):
 pgf.add('lam2', lam2)
 pgf.add('lam3', lam3)
 pgf.add('lam4', lam4)
+pgf.add('lam5', lam5)
+pgf.add('lam6', lam6)
 pgf.write('data/fig_epsilon_rank.dat')

psdr/gp.py

@@ -360,6 +360,7 @@ def _fit(self, ell0):
 		self.L = self._make_L(ell)
 		self.alpha, self.beta = self._log_marginal_likelihood(ell, 
 			return_obj = False, return_grad = False, return_alpha_beta = True)
+		self._ell = ell
 
 	def eval(self, Xnew, return_cov = False):
 		Y = np.dot(self.L, self.X.T).T

psdr/sample/minimax_lloyd.py

@@ -28,21 +28,25 @@ def _update_voronoi_sample(domain, Xhat, X0, L, M0):
 	"""
 
 	# Update the estimate Voronoi vertices
-	X0 = voronoi_vertex_sample(domain, Xhat, X0, L = L)
-	I = unique_points(X0)
-	V = np.zeros((M0, len(domain)))
-
-	k = np.sum(I)
-	V[:k] = X0[I]
-
-	# Perform one step of trying to fill missing data
-	V[k:] = domain.sample(M0 - k)
-	V[k:] = voronoi_vertex_sample(domain, Xhat, V[k:], L = L)
+	V = voronoi_vertex_sample(domain, Xhat, X0, L = L)
 
-	I = unique_points(V)
+	# Remove duplicates
+	V = V[unique_points(V)]
+
+	# we need to ensure every node has at least one voronoi vertex associated with it
+	D = cdist(Xhat, V, L =L)
+	d = np.min(D, axis = 0)
+	for k in range(len(Xhat)):
+		I = np.isclose(D[k], d)
+		if np.sum(I) < len(domain):
+			# Generate perturbations of nearby point
+			X0 = np.outer(np.ones(2*len(domain)), Xhat[k])  + 1e-7*np.random.randn(2*len(domain), len(domain))
+			# This pushes these points to nearby bounded Voronoi vertices
+			Vnew = voronoi_vertex_sample(domain, Xhat, X0, L = L)
+			V = np.vstack([V, Vnew])
 
 	# Remove duplicates
-	return V[I]
+	return V[unique_points(V)]
 
 
 def _update_voronoi_full(domain, Xhat, X0, L, M0):