From 15635ad33402f53dfc0efd37dc0d99474990c78e Mon Sep 17 00:00:00 2001
From: "Jeffrey Hokanson @ Thor" <jeffrey@hokanson.us>
Date: Thu, 18 Apr 2019 13:21:38 -0600
Subject: [PATCH] Got multi-objective sampling working

---
 psdr/sample.py                | 25 ++++++++++++++-------
 tests/test_multiobj_sample.py | 41 +++++++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+), 8 deletions(-)
 create mode 100644 tests/test_multiobj_sample.py

diff --git a/psdr/sample.py b/psdr/sample.py
index 6f70651..833b83b 100644
--- a/psdr/sample.py
+++ b/psdr/sample.py
@@ -53,7 +53,6 @@ def initial_sample(domain, L, Nsamp = int(1e4), Nboundary = 50):
 
 	Lrank = U.shape[1]
 
-
 	if Lrank != len(domain):
 		# Attempt to determine the effective dimension
 		if Lrank == 1:
@@ -86,6 +85,7 @@ def initial_sample(domain, L, Nsamp = int(1e4), Nboundary = 50):
 			# in the active direction
 	else:
 		dim = len(domain)
+		Lcs = []
 
 	if dim == 1:
 		assert len(Lcs) == 2
@@ -96,9 +96,11 @@ def initial_sample(domain, L, Nsamp = int(1e4), Nboundary = 50):
 		X0 = np.vstack([alpha*cs[0] + (1-alpha)*cs[1] for alpha in alphas])	
 		return X0
 		
-	elif dim in [2,3]:
+	elif dim in [2,3] and len(Lcs) > dim:
 		# Construct a convex hull of low-dimensional points
-		assert len(Lcs) > dim, "Insufficient points to describe hyperplane in %d dimensions" % dim
+		#if len(Lcs) > dim:
+		#	# Insufficient points to describe hyperplane in m dimensions, instead return corners
+		#	return cs
 		Ldom = ConvexHullDomain(Lcs, solver = 'CVXOPT')
 		# and sample uniformly here
 		LX0 = Ldom.sample(Nsamp)
@@ -168,15 +170,15 @@ def seq_maximin_sample(domain, Xhat, L = None, Nsamp = int(1e3), X0 = None):
 		Sample from inside the domain
 	"""
 	Xhat = np.array(Xhat)
-	Xhat = np.atleast_2d(Xhat)
-
-	if len(Xhat) < 1:
+	if len(Xhat) == 0:
 		# If we don't have any samples, pick one of the corners
 		if L is None:
 			return domain.corner(np.random.randn(len(domain)))
 		else:
 			_, s, VT = scipy.linalg.svd(L)
 			return domain.corner(VT.T[:,0])
+	
+	Xhat = np.atleast_2d(Xhat)
 
 	# Generate candidate points from the Voronoi vertices
 	if X0 is None:
@@ -187,7 +189,6 @@ def seq_maximin_sample(domain, Xhat, L = None, Nsamp = int(1e3), X0 = None):
 
 	Xcan = voronoi_vertex(domain, Xhat, X0, L = L, randomize = True)
 
-
 	# Compute the Euclidean distance between candidates Xcan and current samples Xhat
 	De = cdist(Xcan, Xhat)
 	if L is not None:
@@ -219,11 +220,19 @@ def seq_maximin_sample(domain, Xhat, L = None, Nsamp = int(1e3), X0 = None):
 def multiobj_seq_maximin_sample(domain, Xhat, Ls, Nsamp = int(1e3)):
 	r""" A multi-objective sequential maximin sampling 
 
+	The goal of this algorithm is to return a new sample that maximizes
+	the distance between samples in *several* different metrics.
+
 
 	A typical use case will have Ls that are of size (1,m)
 	
 	"""
 
+	Xhat = np.array(Xhat)
+	if len(Xhat) == 0:
+		Lall = np.vstack(Ls)
+		return seq_maximin_sample(domain, Xhat, L = Lall, Nsamp = Nsamp) 
+
 	vertices = []
 	it = 0
 	queue = []
@@ -278,7 +287,7 @@ def multiobj_seq_maximin_sample(domain, Xhat, Ls, Nsamp = int(1e3)):
 	return seq_maximin_sample(domain_samp, Xhat, L = Lall, Nsamp = Nsamp)
 
 
-def fill_distance_estimate(domain, Xhat, L = None, Nsamp = int(1e4), X0 = None ):
+def fill_distance_estimate(domain, Xhat, L = None, Nsamp = int(1e3), X0 = None ):
 	r""" Estimate the fill distance of the points Xhat in the domain
 
 	The *fill distance* (Def. 1.4 of [Wen04]_) or *dispersion* [LC05]_
diff --git a/tests/test_multiobj_sample.py b/tests/test_multiobj_sample.py
new file mode 100644
index 0000000..1646a37
--- /dev/null
+++ b/tests/test_multiobj_sample.py
@@ -0,0 +1,41 @@
+from __future__ import print_function
+import numpy as np
+import matplotlib.pyplot as plt
+import psdr
+
+
+def test_multiobj_sample():
+	m = 2
+	L1 = np.ones((1,m))
+	L2 = np.ones((1,m))
+	L2[0,1] = 0
+	dom = psdr.BoxDomain(-np.ones(m), np.ones(m))
+	Ls = [L1, L2]
+
+	Xhat = []
+	for i in range(10):
+		x = psdr.multiobj_seq_maximin_sample(dom, Xhat, Ls, Nsamp = 10)
+		Xhat.append(x)
+	
+	Xhat = np.array(Xhat)
+
+	# Now check fill distance
+	for L in Ls:
+		assert L.shape[0] == 1, "Only 1-d tests implemented"
+		c1 =  dom.corner(L.flatten())
+		c2 =  dom.corner(-L.flatten())
+		lb, ub = sorted([L.dot(c1), L.dot(c2)])
+		vol = float(ub - lb)
+		d = psdr.fill_distance_estimate(dom, Xhat, L = L)
+		print("")
+		print("ideal fill distance ", 0.5*vol/(len(Xhat) - 1) ) 
+		print("actual fill distance", d)
+
+		# we add a fudge factor to ensure the suboptimal sampling passes 
+		assert 0.25*d < 0.5*vol/(len(Xhat)-1), "Sampling not efficient enough"
+
+#	print(np.sort(L1.dot(Xhat.T)).T)
+#	print(np.sort(L2.dot(Xhat.T)).T)
+
+if __name__ == '__main__':
+	test_multiobj_sample()