jeffrey-hokanson / PSDR / 349

pending completion

Build # 349

Build Type

push

travis-ci-com

Committed by

jeffrey-hokanson

Commit Message

Added covering based approach

Run Details

1155 of 1635 branches covered (70.64%)

Branch coverage included in aggregate %.

2 of 2 new or added lines in 1 file covered. (100.0%)

4 existing lines in 1 file now uncovered.

4325 of 4999 relevant lines covered (86.52%)

2.36 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

91.23

/psdr/sample/minimax.py

from __future__ import print_function, division

import numpy as np
import cvxpy as cp
from scipy.spatial.distance import cdist, pdist, squareform

from ..domains.domain import DEFAULT_CVXPY_KWARGS

from .poisson import poisson_disk_sample


def _cq_center_cvxpy(Y, L, q = 10):
        xhat = cp.Variable(L.shape[1])
        obj = cp.sum([cp.norm(L*xhat - y)**q for y in Y])
        #obj = cp.sum([cp.norm(L*xhat - y) for y in Y])
        prob = cp.Problem(cp.Minimize(obj))
        #prob.solve(**DEFAULT_CVXPY_KWARGS)
        prob.solve()
        return np.array(xhat.value)


def _cq_center_agd(X, L, q = 10, maxiter = 1000, verbose = False, xtol = 1e-7, xhat0 = None):
        r""" This computes Cq center using accelerated gradient descent
        See Algorithm 2 in [MJ18]_.

        """
        assert q>= 4, "Require q >= 4"

        # Initialization of the two variables
        if xhat0 is None:
                xhat0 = np.mean(X, axis = 0)
        z = np.copy(xhat0)
        u = np.copy(xhat0)

        lam = 1

        Y = L.dot(X.T).T

        # Compute Lipschitz-smooth constant
        D = squareform(pdist(Y))
        beta = (q-1)*(q-2)*np.max(np.sum(D**(q-2), axis = 1))/(len(X)*(q-2))
        Lnorm = np.linalg.norm(L)         # TODO: Cache this?
        beta *= Lnorm**2                        # This is added b/c beta is proportional to derivative, and LT*L appears in the derivative

        for it in range(maxiter):
                lam_new = (1 + np.sqrt(1 + 4*lam**2))/2.
                gam = (1 - lam)/lam_new

                dZ = np.tile(z, (len(X),1)) - X
                LLdZ = L.T.dot(L.dot(dZ.T)).T

                d = cdist(L.dot(z).reshape(1,-1), Y).T

                grad = np.sum(LLdZ*d**(q-2), axis = 0)/len(X)
                u_new = z - 1./beta*grad
                z_new = (1 - gam)*u_new + gam * u
                dx = np.linalg.norm(z - z_new)

                # update all the variables                
                z = z_new
                lam = lam_new
                u = u_new
                
                if verbose:
                        d = cdist((L.dot(z)).reshape(1,-1), Y)
                        obj = 1/(len(X)*q)*np.sum(d**q)
                        print("%4d | %7.3e | %7.3e" % (it, obj, dx))
                if dx < xtol and it > 10:
                        break
        
        return z

def minimax_cluster(domain, N, L = None, maxiter = 30, N0 = None, xtol = 1e-5, verbose = True, q = 10):
        r"""Identifies an approximate minimax design using a clustering technique due to Mak and Joseph

        This function implements a clustering based approach for minimax sampling following [MJ18]_.


        References
        ----------
        .. [MJ18] Simon Mak and V. Roshan Joseph.
                Minimax and Minimax Projection Designs Using Clustering.
                 Journal of Computational and Graphical Statistics. 2018, vol 27:1 pp 166-178
                DOI:10.1080/10618600.2017.1302881
        """

        if N0 is None:
                N0 = min(N*100, int(1e4))
        N0 = int(N0)

        if L is None:
                L = np.eye(len(domain))

        # Samples from the domain to cluster
        # TODO: Should these be distributed with respect to the L norm?
        X = domain.sample(N0)
        Y = L.dot(X.T).T

        # Initial cluster centers
        # We use the first N random points so that we don't end up with 
        # empty regions 
        Xhat = X[0:N]
        #Xhat = domain.sample(N)

        dx = 0
        I_old = np.zeros(len(X))
        for it in range(maxiter):
                Yhat = L.dot(Xhat.T).T
                # Assign each point to its nearest neighbor in L-norm
                D = cdist(Yhat, Y)
                I = np.argmin(D, axis = 0)
                if verbose:
                        print("%4d obj: %7.3e; dx %7.3e" % (it, np.max(np.min(D, axis = 0)), dx)) 
                if np.all(I_old == I):
                        if verbose: print("point sets unchanged")
                #        break
                
                dx = 0
                for i in range(N):
                        #xhat = _cq_center_cvxpy(Y[I == i], L, q = 10)
                        #Ii = (I == i)
                        #if np.all(I_old[Ii] == I[Ii]):
                        #        # If we haven't changed the elements of this set between iterations
                        #        xhat = _cq_center_agd(X[Ii], L, q = q, xhat0 = Xhat[i])
                        #else:
                        #        xhat = _cq_center_agd(X[Ii], L, q = q, xhat0 = Xhat[i])
                        xhat = _cq_center_agd(X[I == i], L, q = q, xhat0 = Xhat[i])
                        dx = max(dx, np.linalg.norm(xhat - Xhat[i]))
                        Xhat[i] = xhat

                if dx < xtol:
                        print('stopped due to small movement')
                        break

                I_old = I

        return Xhat

def minimax_covering(domain, r, L = None, **kwargs):
        r"""Approximate a minimax design using a discrete approximation of the domain.

        This function is a convience wrapper around :meth:`psdr.minimax_covering_discrete`:
        it constructs a discrete approximation of the domain using Poisson disk sampling
        and then selects from those samples points such that the maximum distance between
        any point in the domain and the samples is at most :math:`r`.

        *Note* This method will solve an expensive 0-1 linear program and will not scale well
        if more than a few hundred Poisson disk samples are taken

        Parameters
        ----------
        domain: Domain
                Domain from which to sample
        r: float
                The target maximum distance between any point in the domain and
                the returned sample set.
        L: array-like (?,m)
                Weighting matrix on the two-norm
        **kwargs: dictionary
                Additional parameters for :meth:`psdr.minimax_covering_discrete`
        
        Returns
        -------
        X: np.array
                Samples from the domain approximating a minimax sampling
        """
        # This ensures that no point in the domain is more than r/2 away from one of these samples
        X = poisson_disk_sample(domain, r/4., L = L)
        # This then selects 
        I = minimax_covering_discrete(X, r/2., L = L, **kwargs)
        return X[I]


def minimax_covering_discrete(X, r, L = None, **kwargs):
        r""" Constructs a minimax design on discrete domain by solving a covering problem

        This implements an algorithm due to Tan [Tan13]_ which, given a finite number of points,
        
        
        References
        ----------
        .. [Tan13] Minimax Designs for Finite Design Regions
                Matthias H. Y. Tan. 
                Technometrics 55:3, 346-358
                https://doi.org/10.1080/00401706.2013.804439
        """
        X = np.array(X)
        if L is None:
                L = np.eye(X.shape[1])
        # Index set - 1/true if the node is selected
        I = cp.Variable(len(X), boolean = True) 

        # compute the pairwise distance matrix
        D = squareform(pdist(L.dot(X.T).T))
        # we require that every 
        constraints = [ I.__rmatmul__(D < r) >= np.ones(len(X)) ]

        prob = cp.Problem(cp.Minimize(cp.sum(I)), constraints)
        prob.solve(**kwargs)

        # Convert to a boolean array
        I = np.array(I.value > 0.5, dtype = np.bool)
        return I

1	from __future__ import print_function, division	5✔
2
3	import numpy as np	5✔
4	import cvxpy as cp	5✔
5	from scipy.spatial.distance import cdist, pdist, squareform	5✔
6
7	from ..domains.domain import DEFAULT_CVXPY_KWARGS	5✔
8
9	from .poisson import poisson_disk_sample	5✔
10
11
12	def _cq_center_cvxpy(Y, L, q = 10):	5✔
13	xhat = cp.Variable(L.shape[1])	5✔
14	obj = cp.sum([cp.norm(Lxhat - y)*q for y in Y])	5✔
15	#obj = cp.sum([cp.norm(L*xhat - y) for y in Y])
16	prob = cp.Problem(cp.Minimize(obj))	5✔
17	#prob.solve(**DEFAULT_CVXPY_KWARGS)
18	prob.solve()	5✔
19	return np.array(xhat.value)	5✔
20
21
22	def _cq_center_agd(X, L, q = 10, maxiter = 1000, verbose = False, xtol = 1e-7, xhat0 = None):	5✔
23	r""" This computes Cq center using accelerated gradient descent
24	See Algorithm 2 in [MJ18]_.
25
26	"""
27	assert q>= 4, "Require q >= 4"	5✔
28
29	# Initialization of the two variables
30	if xhat0 is None:	5✔
31	xhat0 = np.mean(X, axis = 0)	5✔
32	z = np.copy(xhat0)	5✔
33	u = np.copy(xhat0)	5✔
34
35	lam = 1	5✔
36
37	Y = L.dot(X.T).T	5✔
38
39	# Compute Lipschitz-smooth constant
40	D = squareform(pdist(Y))	5✔
41	beta = (q-1)(q-2)np.max(np.sum(D*(q-2), axis = 1))/(len(X)(q-2))	5✔
42	Lnorm = np.linalg.norm(L) # TODO: Cache this?	5✔
43	beta = Lnorm2 # This is added b/c beta is proportional to derivative, and LTL appears in the derivative	5✔
44
45	for it in range(maxiter):	5✔
46	lam_new = (1 + np.sqrt(1 + 4lam*2))/2.	5✔
47	gam = (1 - lam)/lam_new	5✔
48
49	dZ = np.tile(z, (len(X),1)) - X	5✔
50	LLdZ = L.T.dot(L.dot(dZ.T)).T	5✔
51
52	d = cdist(L.dot(z).reshape(1,-1), Y).T	5✔
53
54	grad = np.sum(LLdZd*(q-2), axis = 0)/len(X)	5✔
55	u_new = z - 1./beta*grad	5✔
56	z_new = (1 - gam)u_new + gam u	5✔
57	dx = np.linalg.norm(z - z_new)	5✔
58
59	# update all the variables
60	z = z_new	5✔
61	lam = lam_new	5✔
62	u = u_new	5✔
63
64	if verbose:	5!
UNCOV 65	d = cdist((L.dot(z)).reshape(1,-1), Y)	×
UNCOV 66	obj = 1/(len(X)q)np.sum(d**q)	×
UNCOV 67	print("%4d \| %7.3e \| %7.3e" % (it, obj, dx))	×
68	if dx < xtol and it > 10:	5✔
69	break	5✔
70
71	return z	5✔
72
73	def minimax_cluster(domain, N, L = None, maxiter = 30, N0 = None, xtol = 1e-5, verbose = True, q = 10):	5✔
74	r"""Identifies an approximate minimax design using a clustering technique due to Mak and Joseph
75
76	This function implements a clustering based approach for minimax sampling following [MJ18]_.
77
78
79	References
80	----------
81	.. [MJ18] Simon Mak and V. Roshan Joseph.
82	Minimax and Minimax Projection Designs Using Clustering.
83	Journal of Computational and Graphical Statistics. 2018, vol 27:1 pp 166-178
84	DOI:10.1080/10618600.2017.1302881
85	"""
86
87	if N0 is None:	5!
UNCOV 88	N0 = min(N*100, int(1e4))	×
89	N0 = int(N0)	5✔
90
91	if L is None:	5!
92	L = np.eye(len(domain))	5✔
93
94	# Samples from the domain to cluster
95	# TODO: Should these be distributed with respect to the L norm?
96	X = domain.sample(N0)	5✔
97	Y = L.dot(X.T).T	5✔
98
99	# Initial cluster centers
100	# We use the first N random points so that we don't end up with
101	# empty regions
102	Xhat = X[0:N]	5✔
103	#Xhat = domain.sample(N)
104
105	dx = 0	5✔
106	I_old = np.zeros(len(X))	5✔
107	for it in range(maxiter):	5!
108	Yhat = L.dot(Xhat.T).T	5✔
109	# Assign each point to its nearest neighbor in L-norm
110	D = cdist(Yhat, Y)	5✔
111	I = np.argmin(D, axis = 0)	5✔
112	if verbose:	5!
113	print("%4d obj: %7.3e; dx %7.3e" % (it, np.max(np.min(D, axis = 0)), dx))	5✔
114	if np.all(I_old == I):	5✔
115	if verbose: print("point sets unchanged")	5✔
116	# break
117
118	dx = 0	5✔
119	for i in range(N):	5✔
120	#xhat = _cq_center_cvxpy(Y[I == i], L, q = 10)
121	#Ii = (I == i)
122	#if np.all(I_old[Ii] == I[Ii]):
123	# # If we haven't changed the elements of this set between iterations
124	# xhat = _cq_center_agd(X[Ii], L, q = q, xhat0 = Xhat[i])
125	#else:
126	# xhat = _cq_center_agd(X[Ii], L, q = q, xhat0 = Xhat[i])
127	xhat = _cq_center_agd(X[I == i], L, q = q, xhat0 = Xhat[i])	5✔
128	dx = max(dx, np.linalg.norm(xhat - Xhat[i]))	5✔
129	Xhat[i] = xhat	5✔
130
131	if dx < xtol:	5✔
132	print('stopped due to small movement')	5✔
133	break	5✔
134
135	I_old = I	5✔
136
137	return Xhat	5✔
138
139	def minimax_covering(domain, r, L = None, **kwargs):	5✔
140	r"""Approximate a minimax design using a discrete approximation of the domain.
141
142	This function is a convience wrapper around :meth:`psdr.minimax_covering_discrete`:
143	it constructs a discrete approximation of the domain using Poisson disk sampling
144	and then selects from those samples points such that the maximum distance between
145	any point in the domain and the samples is at most :math:`r`.
146
147	Note This method will solve an expensive 0-1 linear program and will not scale well
148	if more than a few hundred Poisson disk samples are taken
149
150	Parameters
151	----------
152	domain: Domain
153	Domain from which to sample
154	r: float
155	The target maximum distance between any point in the domain and
156	the returned sample set.
157	L: array-like (?,m)
158	Weighting matrix on the two-norm
159	**kwargs: dictionary
160	Additional parameters for :meth:`psdr.minimax_covering_discrete`
161
162	Returns
163	-------
164	X: np.array
165	Samples from the domain approximating a minimax sampling
166	"""
167	# This ensures that no point in the domain is more than r/2 away from one of these samples
168	X = poisson_disk_sample(domain, r/4., L = L)	5✔
169	# This then selects
170	I = minimax_covering_discrete(X, r/2., L = L, **kwargs)	5✔
171	return X[I]	5✔
172
173
174	def minimax_covering_discrete(X, r, L = None, **kwargs):	5✔
175	r""" Constructs a minimax design on discrete domain by solving a covering problem
176
177	This implements an algorithm due to Tan [Tan13]_ which, given a finite number of points,
178
179
180	References
181	----------
182	.. [Tan13] Minimax Designs for Finite Design Regions
183	Matthias H. Y. Tan.
184	Technometrics 55:3, 346-358
185	https://doi.org/10.1080/00401706.2013.804439
186	"""
187	X = np.array(X)	5✔
188	if L is None:	5!
189	L = np.eye(X.shape[1])	5✔
190	# Index set - 1/true if the node is selected
191	I = cp.Variable(len(X), boolean = True)	5✔
192
193	# compute the pairwise distance matrix
194	D = squareform(pdist(L.dot(X.T).T))	5✔
195	# we require that every
196	constraints = [ I.__rmatmul__(D < r) >= np.ones(len(X)) ]	5✔
197
198	prob = cp.Problem(cp.Minimize(cp.sum(I)), constraints)	5✔
199	prob.solve(**kwargs)	5✔
200
201	# Convert to a boolean array
202	I = np.array(I.value > 0.5, dtype = np.bool)	5✔
203	return I	5✔

jeffrey-hokanson / PSDR / 349

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous