From 99571de90f737441f54dfc300946c90a879674ce Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Fri, 8 Jul 2016 22:58:42 -0700 Subject: [PATCH 1/9] Hyperparamopt package added with tests --- brainiak/hyperparamopt/__init__.py | 4 + brainiak/hyperparamopt/hpo.py | 224 +++++++++++++++++++++++++++++ brainiak/hyperparamopt/mcmc.py | 111 ++++++++++++++ brainiak/hyperparamopt/norm.pyx | 29 ++++ requirements-dev.txt | 3 + setup.cfg | 1 + tests/hyperparamopt/test_hpo.py | 37 +++++ tests/hyperparamopt/test_mcmc.py | 35 +++++ 8 files changed, 444 insertions(+) create mode 100644 brainiak/hyperparamopt/__init__.py create mode 100644 brainiak/hyperparamopt/hpo.py create mode 100644 brainiak/hyperparamopt/mcmc.py create mode 100644 brainiak/hyperparamopt/norm.pyx create mode 100644 tests/hyperparamopt/test_hpo.py create mode 100644 tests/hyperparamopt/test_mcmc.py diff --git a/brainiak/hyperparamopt/__init__.py b/brainiak/hyperparamopt/__init__.py new file mode 100644 index 000000000..7867be74a --- /dev/null +++ b/brainiak/hyperparamopt/__init__.py @@ -0,0 +1,4 @@ +""" Hyper parameter optimization package """ + +import pyximport +pyximport.install() diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py new file mode 100644 index 000000000..ed68a5c9f --- /dev/null +++ b/brainiak/hyperparamopt/hpo.py @@ -0,0 +1,224 @@ +"""Hyper Parameter Optimization (HPO) + +This implementation is based on the work: + +.. [Bergstra2011] "Algorithms for Hyper-Parameter Optimization", + James S. Bergstra and Bardenet, R\'{e}mi and Bengio, Yoshua + and Bal\'{a}zs K\'{e}gl. NIPS 2011 + +.. [Bergstra2013] "Making a Science of Model Search: + Hyperparameter Optimization in Hundreds of Dimensions for + Vision Architectures", James Bergstra, Daniel Yamins, David Cox. + JMLR W&CP 28 (1) : 115–123, 2013 + +""" + +# Authors: Narayanan Sundaram (Intel Labs) + +import logging +from .mcmc import get_multichain_samples +from .norm import getgmmpdf +import numpy as np +from scipy.special import erf +from tqdm import tqdm + + +logger = logging.getLogger(__name__) + + +def getsigma(x, minlimit=-np.inf, maxlimit=np.inf): + z = np.append(x, [minlimit, maxlimit]) + sigma = np.ones(x.shape) + for i in range(x.size): + xleft = z[np.argmin([(x[i] - k) if k < x[i] else np.inf for k in z])] + xright = z[np.argmin([(k - x[i]) if k > x[i] else np.inf for k in z])] + sigma[i] = max(x[i] - xleft, xright - x[i]) + if sigma[i] == np.inf: + sigma[i] = min(x[i] - xleft, xright - x[i]) + if (sigma[i] == -np.inf): + sigma[i] = 1.0 + return sigma + + +class gmm_1d_distribution: + """GMM 1D distribution. + + Given a set of points, we create this object so that we + can calculate likelihoods and generate samples from this + 1D Gaussian mixture model. + + Parameters + ---------- + + x : 1D array + Set of points to create the GMM + + minlimit : double, default : -inf + Minimum limit for the distribution + + maxlimit : double, default : +inf + Maximum limit for the distribution + + weights : double scalar or 1D array with same size as x, default 1.0 + Used to weight the points non-uniformly if required + """ + + def __init__(self, x, minlimit=-np.inf, maxlimit=np.inf, weights=1.0): + self.points = x + self.N = x.size + self.minlimit = minlimit + self.maxlimit = maxlimit + self.sigma = getsigma(x, minlimit=minlimit, maxlimit=maxlimit) + self.weights = 2. / (erf((maxlimit - x) + / (np.sqrt(2.) * self.sigma)) + - erf((minlimit - x) + / (np.sqrt(2.) * self.sigma))) * weights + # return self + + def __call__(self, xt): + if (np.isscalar(xt)): + return getgmmpdf(xt, self.points, self.sigma, self.weights, + self.minlimit, self.maxlimit) + else: + return np.array([getgmmpdf(t, self.points, self.sigma, + self.weights, self.minlimit, + self.maxlimit) for t in xt]) + + def get_samples(self, chains=1, points_per_chain=1): + pts = get_multichain_samples(N=points_per_chain, + p=self, nchains=chains) + return pts + + +def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf, show_plot=False): + z = np.array(list(zip(x, y)), dtype=np.dtype([('x', float), ('y', float)])) + z = np.sort(z, order='y') + n = y.shape[0] + g = int(np.round(np.ceil(0.15 * n))) + ldata = z[0:g] + gdata = z[g:n] + lymin = ldata['y'].min() + lymax = ldata['y'].max() + weights = (lymax - ldata['y']) / (lymax - lymin) + lx = gmm_1d_distribution(ldata['x'], minlimit=minlimit, + maxlimit=maxlimit, weights=weights) + gx = gmm_1d_distribution(gdata['x'], minlimit=minlimit, maxlimit=maxlimit) + + samples = lx.get_samples(chains=10, points_per_chain=100) + ei = lx(samples) / gx(samples) + + if show_plot is True: + import pylab as plt + plt.scatter(samples, lx(samples), color='r') + plt.scatter(samples, gx(samples), color='b') + plt.scatter(samples, ei, color='g') + plt.show() + + h = (x.max() - x.min()) / (10 * x.size) + # assumes prior of x is uniform -- should change for different priors + # d = np.abs(x - samples[ei.argmax()]).min() + # CDF(x+d/2) - CDF(x-d/2) < 1/(10*x.size) then reject else accept + s = 0 + while (np.abs(x - samples[ei.argmax()]).min() < h): + ei[ei.argmax()] = 0 + s = s + 1 + if (s == samples.size): + break + xnext = samples[ei.argmax()] + + return xnext + + +def getSample(x, y, dist, minlimit=-np.inf, maxlimit=np.inf): + if (dist == 'GMM'): + return getNextSample(x, y, minlimit, maxlimit) + if (dist == 'uniform'): + return np.random.random() * (maxlimit - minlimit) + minlimit + if (dist == 'loguniform'): + return np.exp(np.random.random() + * (np.log(maxlimit) - np.log(minlimit)) + + np.log(minlimit)) + else: + logger.error('Unsupported distribution for variable') + + +def fmin(lossfn, + space, + algo, + maxevals, + trials, + init_random_evals=30, + explore_prob=0.2, + verbose=False): + """Find the minimum of function through hyper paramter optimization + + Arguments + --------- + + lossfn : function that takes in a dictionary and returns a real value + Function to be minimized + + space : Dictionary specifying the range and distribution of + the hyperparamters + + algo : Algo to be used (ignored, can use None) + + maxevals : int + Maximum number of evaluations of lossfn allowed + + trials : list + Holds the output of the optimization trials + Need not be empty to begin with, new trials are appended + at the end + + init_random_evals : int, default 30 + Number of random trials to initialize the + optimization + + explore_prob : double in [0, 1], default 0.2 + Controls the exploration-vs-exploitation ratio + Currently 20% of trails are random samples + + verbose : bool, default False + Get information on current point being processed + + Returns + ------- + + best : trial entry (dictionary of hyperparameters) + Best hyperparameter setting found + """ + + if (len(trials) > init_random_evals): + init_random_evals = 0 + + for t in tqdm(range(maxevals)): + sdict = {} + + if (t >= init_random_evals and np.random.random() > explore_prob): + search_algo = 'Exploit' + else: + search_algo = 'Explore' + + yarray = np.array([tr['loss'] for tr in trials]) + for s in space: + sarray = np.array([tr[s] for tr in trials]) + dist = 'GMM' if (search_algo == 'Exploit') else space[s]['dist'] + sdict[s] = getSample(sarray, yarray, dist, + minlimit=space[s]['lo'], + maxlimit=space[s]['hi']) + + if (verbose): + logger.info(search_algo) + logger.info('Next point ', t, ' = ', sdict) + + y = lossfn(sdict) + sdict['loss'] = y + trials.append(sdict) + + yarray = np.array([tr['loss'] for tr in trials]) + yargmin = yarray.argmin() + + if (verbose): + logger.info('Best point so far = ', trials[yargmin]) + return trials[yargmin] diff --git a/brainiak/hyperparamopt/mcmc.py b/brainiak/hyperparamopt/mcmc.py new file mode 100644 index 000000000..0318c4fcd --- /dev/null +++ b/brainiak/hyperparamopt/mcmc.py @@ -0,0 +1,111 @@ +import numpy as np +import scipy.stats as st +import logging + + +def candidate(x): + """Generates candidate around point x + + Returns + ------- + + sample from q(x*|x) - Unit Normal distribution around x + """ + + return np.random.standard_normal() + x + + +def candidate_dist(x, xp): # return value of q(xp| x) + return st.norm.pdf(x - xp) + + +def check_accept(xcurr, xprop, p): + return min(1.0, p(xprop) / p(xcurr)) + + +def get_next(x, p, n): + xnext = np.zeros(n) + xnext[-1] = x + pxcurr = p(x) + for i in range(n): + xp = candidate(xnext[i - 1]) + pxp = p(xp) + if (np.random.random() < pxp / pxcurr): + xnext[i] = xp + pxcurr = pxp + else: + xnext[i] = xnext[i - 1] + return xnext + + +def get_chain(N, p, burn_in=2000): + """Get a sequence of numbers sampled from a single chain + of MCMC (Metropolis-Hastings) sampler + + Arguments + --------- + + N : int + Number of samples required + + p : function that returns a pdf value at any real number + Distribution that needs to be sampled + + burn_in : int, default 2000 + Number of burn-in (discarded) samples + + Returns + ------- + + samples : 1D array, shape [N] + Samples generated from MCMC sampler (should resemble samples from p(x)) + """ + + x = np.ones(N) + x0 = np.random.standard_normal() * 100 + while(p(x0) <= np.finfo(np.double).eps * 10): + x0 = np.random.standard_normal() * 100 + if (p(x0) <= 0): + logging.error('Markov chain failed to initialize properly \ + - Values probably very far from origin') + + # burn in iterations + x0 = get_next(x0, p, burn_in)[-1] + + # actual iterations + x = get_next(x0, p, N) + return x + + +def get_multichain_samples(N, p, nchains=3, burn_in=2000): + """Get a sequence of numbers sampled from multiple chains + of MCMC (Metropolis-Hastings) sampler + + Arguments + --------- + + N : int + Number of samples per chain required + + p : function that returns a pdf value at any real number + Distribution that needs to be sampled + + nchains : int, default 3 + Number of independent MCMC chains to sample + + burn_in : int, default 2000 + Number of burn-in (discarded) samples + + Returns + ------- + + samples : 1D array, shape [nchains*N] + Samples generated from MCMC sampler (should resemble samples from p(x)) + """ + + pts = np.zeros(nchains * N) + for c in range(nchains): + xp = get_chain(N * 3, p, burn_in=burn_in) + pts[c * N: (c + 1) * N] =\ + xp[np.random.choice(N * 3, N)] # pick N points at random + return pts diff --git a/brainiak/hyperparamopt/norm.pyx b/brainiak/hyperparamopt/norm.pyx new file mode 100644 index 000000000..6f9cbe72d --- /dev/null +++ b/brainiak/hyperparamopt/norm.pyx @@ -0,0 +1,29 @@ +#cython: embedsignature=True + +import numpy +cimport numpy as np + +cdef extern from "math.h": + double exp(double x) + double sqrt(double x) + + +pi = numpy.pi + + +cpdef double norm_pdf(double x, double mu, double sigma): + cdef double z + z = (x-mu)/sigma + return exp(-0.5*z*z)/sqrt(2.0*pi)/sigma + + +cpdef double getgmmpdf(double xt,np.ndarray[np.float64_t, ndim=1] x, np.ndarray[np.float64_t, ndim=1] sigma, np.ndarray[np.float64_t, ndim=1] weights, double minlimit, double maxlimit): + cdef double y + y = 0 + if (xt < minlimit): + return 0 + if (xt > maxlimit): + return 0 + for _x in range(x.size): + y += norm_pdf(xt, x[_x], sigma[_x])*weights[_x]/x.size + return y diff --git a/requirements-dev.txt b/requirements-dev.txt index 6c8c1dc27..141c80895 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,8 +1,11 @@ +cython flake8 flake8-print notebook pytest pytest-cov +pytest-cython restructuredtext-lint sphinx sphinx_rtd_theme +tqdm diff --git a/setup.cfg b/setup.cfg index c643accec..d924e617c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -7,6 +7,7 @@ addopts = --cov-report=xml --cov-report=html --cov-report=term + --doctest-cython [coverage:run] branch = True diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py new file mode 100644 index 000000000..a26a85fef --- /dev/null +++ b/tests/hyperparamopt/test_hpo.py @@ -0,0 +1,37 @@ +import pytest + + +def test_simple_gmm(): + from brainiak.hyperparamopt.hpo import gmm_1d_distribution + import numpy as np + + x = np.array([1., 1., 2., 3., 1.]) + d = gmm_1d_distribution(x, minlimit=0., maxlimit=4.) + assert d(1.1) > d(3.5), "GMM distribution not behaving correctly" + assert d(2.0) > d(3.0), "GMM distribution not behaving correctly" + + samples = d.get_samples(chains=3, points_per_chain=10) + np.testing.assert_array_less(samples, 4.) + np.testing.assert_array_less(0., samples) + +def test_simple_hpo(): + from brainiak.hyperparamopt.hpo import fmin + import numpy as np + + def f(args): + x = args['x'] + return x*x + + s = {'x': {'dist': 'uniform', 'lo': -10., 'hi': 10.}} + trials = [] + best = fmin(lossfn=f, space=s, algo=None, maxevals=100, trials=trials) + + yarray = np.array([tr['loss'] for tr in trials]) + np.testing.assert_array_less(yarray, 100.) + + xarray = np.array([tr['x'] for tr in trials]) + np.testing.assert_array_less(np.abs(xarray), 10.) + + assert best['loss'] < 100., "HPO out of range" + assert np.abs(best['x']) < 10., "HPO out of range" + assert np.abs(best['x']) < 1., "HPO not accurate" diff --git a/tests/hyperparamopt/test_mcmc.py b/tests/hyperparamopt/test_mcmc.py new file mode 100644 index 000000000..7615938d5 --- /dev/null +++ b/tests/hyperparamopt/test_mcmc.py @@ -0,0 +1,35 @@ +import pytest + +def test_get_chain(): + from brainiak.hyperparamopt.mcmc import get_multichain_samples + import numpy as np + import numpy.testing as npt + import scipy.stats as st + + def normal(mean, std): + def f(x): + return st.norm.pdf(x, loc=mean, scale=std) + return f + + for mean,std in [(0.,1.), (2.,4.), (-5., 3.)]: + p = normal(mean, std) + samples = get_multichain_samples(1000, p, nchains=5) + assert(np.abs(mean - np.mean(samples)) <= 1.) + assert(np.abs(std - np.std(samples)) <= 1.) + # assert(st.skewtest(samples).pvalue >= 0.05) + # assert(st.kurtosistest(samples).pvalue >= 0.05) + + + """ + #import pylab as plt + #plt.hist(samples, 100) + #plt.show() + # Anderson-Darling test + A2, criticalvalues, significancelevel = st.anderson(samples, 'norm') + print(criticalvalues, significancelevel, A2) + + # critical values at [15, 10, 5, 2.5, 1] + for i in range(len(significancelevel)): + if (significancelevel[i] == 5.): # at 5% significance level + assert(A2 <= criticalvalues[i]) + """ From 5822d39d18d7d1dbbc5fa7ed31d65622c3d35ddf Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Mon, 11 Jul 2016 23:06:31 -0700 Subject: [PATCH 2/9] Added example, better tests, more comments --- brainiak/hyperparamopt/hpo.py | 7 ++-- brainiak/hyperparamopt/mcmc.py | 9 +++++ brainiak/hyperparamopt/norm.pyx | 21 +++++++++++- examples/hpo_example.py | 58 +++++++++++++++++++++++++++++++++ tests/hyperparamopt/test_hpo.py | 27 +++++++++++++-- 5 files changed, 113 insertions(+), 9 deletions(-) create mode 100644 examples/hpo_example.py diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index ed68a5c9f..e6e71dd44 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -84,9 +84,9 @@ def __call__(self, xt): self.weights, self.minlimit, self.maxlimit) for t in xt]) - def get_samples(self, chains=1, points_per_chain=1): + def get_samples(self, chains=1, points_per_chain=1, burn_in=2000): pts = get_multichain_samples(N=points_per_chain, - p=self, nchains=chains) + p=self, nchains=chains, burn_in=burn_in) return pts @@ -144,7 +144,6 @@ def getSample(x, y, dist, minlimit=-np.inf, maxlimit=np.inf): def fmin(lossfn, space, - algo, maxevals, trials, init_random_evals=30, @@ -161,8 +160,6 @@ def fmin(lossfn, space : Dictionary specifying the range and distribution of the hyperparamters - algo : Algo to be used (ignored, can use None) - maxevals : int Maximum number of evaluations of lossfn allowed diff --git a/brainiak/hyperparamopt/mcmc.py b/brainiak/hyperparamopt/mcmc.py index 0318c4fcd..0aeb6f5cb 100644 --- a/brainiak/hyperparamopt/mcmc.py +++ b/brainiak/hyperparamopt/mcmc.py @@ -1,3 +1,12 @@ +"""Metropolis-Hasting Random number generator + +This implementation provides random samples from a user-given +probability density function through the Metropolis-Hasting algorithm. + +""" + +# Authors: Narayanan Sundaram (Intel Labs) + import numpy as np import scipy.stats as st import logging diff --git a/brainiak/hyperparamopt/norm.pyx b/brainiak/hyperparamopt/norm.pyx index 6f9cbe72d..6839b1d86 100644 --- a/brainiak/hyperparamopt/norm.pyx +++ b/brainiak/hyperparamopt/norm.pyx @@ -1,5 +1,11 @@ #cython: embedsignature=True +"""Cython file for optimizing GMM likelihood computation + +""" + +# Authors: Narayanan Sundaram (Intel Labs) + import numpy cimport numpy as np @@ -12,18 +18,31 @@ pi = numpy.pi cpdef double norm_pdf(double x, double mu, double sigma): + """Calculate Gaussian pdf + + Given x, returns exp(-0.5*z*z)/(sigma*sqrt(2.*pi)) where + z = (x-mu)/sigma + """ + cdef double z z = (x-mu)/sigma return exp(-0.5*z*z)/sqrt(2.0*pi)/sigma cpdef double getgmmpdf(double xt,np.ndarray[np.float64_t, ndim=1] x, np.ndarray[np.float64_t, ndim=1] sigma, np.ndarray[np.float64_t, ndim=1] weights, double minlimit, double maxlimit): + """Calculates the 1D GMM likelihood + + y = \sum_{i=1}^{N} norm_pdf(x, x_i, sigma_i)/(\sum weight_i) + """ + cdef double y + cdef double w + w = sum(weights) y = 0 if (xt < minlimit): return 0 if (xt > maxlimit): return 0 for _x in range(x.size): - y += norm_pdf(xt, x[_x], sigma[_x])*weights[_x]/x.size + y += norm_pdf(xt, x[_x], sigma[_x])*weights[_x]/w return y diff --git a/examples/hpo_example.py b/examples/hpo_example.py new file mode 100644 index 000000000..f4a0af759 --- /dev/null +++ b/examples/hpo_example.py @@ -0,0 +1,58 @@ +import numpy as np +import matplotlib.pyplot as plt +import brainiak.hyperparamopt.hpo as hpo + +def branin(x1, x2): + a = 1.0 + b = 5.1/(4*np.pi*np.pi) + c = 5.0/np.pi + r = 6.0 + s = 10.0 + t = 1.0/(8*np.pi) + return a*((x2 - b*x1*x1 + c*x1 - r)**2) + s*(1-t)*np.cos(x1) + s + 5*x1 + +def g(args): + x1 = args['x1'] + x2 = args['x2'] + return branin(x1,x2) + +x1lo = -5 +x1hi = 10 +x2lo = 0 +x2hi = 15 + +space = {'x1':{'dist':'uniform', 'lo':x1lo, 'hi':x1hi}, + 'x2':{'dist':'uniform', 'lo':x2lo, 'hi':x2hi}} +trials = [] +n_hpo_samples = 100 + +best = hpo.fmin(lossfn=g, space=space, maxevals=n_hpo_samples, trials=trials, verbose=False) +print("Best obtained through HPO (", n_hpo_samples, " samples) = ", + best['x1'], best['x2'], "; min value = ", best['loss']) + +nt = 100 +x1t = np.linspace(x1lo, x1hi, nt) +x2t = np.linspace(x2lo, x2hi, nt) +x1m, x2m = np.meshgrid(x1t, x2t) +z = branin(x1m, x2m) +print("Best obtained through grid search (", nt*nt, " samples) = ", + x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], + "; min value = ", z.min()) + +x1 = np.array([tr['x1'] for tr in trials]) +x2 = np.array([tr['x2'] for tr in trials]) +y = np.array([tr['loss'] for tr in trials]) + +h = (z.max()-z.min())/25 +plt.contour(x1m, x2m, z, levels=np.linspace(z.min()-h, z.max(), 26)) +plt.scatter(x1, x2, s=10, color='r', label='HPO Samples') +plt.xlabel('x1') +plt.ylabel('x2') +plt.title('Hyperparamter optimization using HPO') + +plt.scatter(best['x1'], best['x2'], s=30, color='b', label='Best HPO') +plt.scatter(x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], + s=30, color='g', label='Best grid search') +plt.legend() +plt.show() + diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py index a26a85fef..2fc7c204c 100644 --- a/tests/hyperparamopt/test_hpo.py +++ b/tests/hyperparamopt/test_hpo.py @@ -9,11 +9,33 @@ def test_simple_gmm(): d = gmm_1d_distribution(x, minlimit=0., maxlimit=4.) assert d(1.1) > d(3.5), "GMM distribution not behaving correctly" assert d(2.0) > d(3.0), "GMM distribution not behaving correctly" + assert d(-1.0) == 0, "GMM distribution out of bounds error" + assert d(9.0) == 0, "GMM distribution out of bounds error" - samples = d.get_samples(chains=3, points_per_chain=10) + samples = d.get_samples(chains=2, points_per_chain=10, burn_in=50) np.testing.assert_array_less(samples, 4.) np.testing.assert_array_less(0., samples) +def test_simple_gmm_weights(): + from brainiak.hyperparamopt.hpo import gmm_1d_distribution + import numpy as np + + x = np.array([1., 1., 2., 3., 1., 3.]) + d = gmm_1d_distribution(x) + + x2 = np.array([1., 2., 3.]) + w = np.array([3., 1., 2.]) + d2 = gmm_1d_distribution(x2, weights=w) + y2 = d2(np.array([1.1, 2.0])) + + assert(d2(1.1) == y2[0], + "GMM distribution array & scalar results don't match") + assert(np.abs(d(1.1) - d2(1.1)) < 1e-5, + "GMM distribution weights not handled correctly") + assert(np.abs(d(2.0) - d2(2.0)) < 1e-5, + "GMM distribution weights not handled correctly") + + def test_simple_hpo(): from brainiak.hyperparamopt.hpo import fmin import numpy as np @@ -24,7 +46,7 @@ def f(args): s = {'x': {'dist': 'uniform', 'lo': -10., 'hi': 10.}} trials = [] - best = fmin(lossfn=f, space=s, algo=None, maxevals=100, trials=trials) + best = fmin(lossfn=f, space=s, maxevals=50, trials=trials, verbose=True) yarray = np.array([tr['loss'] for tr in trials]) np.testing.assert_array_less(yarray, 100.) @@ -34,4 +56,3 @@ def f(args): assert best['loss'] < 100., "HPO out of range" assert np.abs(best['x']) < 10., "HPO out of range" - assert np.abs(best['x']) < 1., "HPO not accurate" From 0d152610243adf8c074140f913bcf1d7939a351d Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Tue, 12 Jul 2016 00:34:39 -0700 Subject: [PATCH 3/9] Added license paragraph to files, coverage > 90%, more tests --- brainiak/hyperparamopt/hpo.py | 30 +++++++++++++++++++---------- brainiak/hyperparamopt/mcmc.py | 13 +++++++++++++ brainiak/hyperparamopt/norm.pyx | 16 ++++++++++++++-- examples/hpo_example.py | 13 +++++++++++++ tests/hyperparamopt/test_hpo.py | 33 +++++++++++++++++++++++++++++++- tests/hyperparamopt/test_mcmc.py | 13 +++++++++++++ 6 files changed, 105 insertions(+), 13 deletions(-) diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index e6e71dd44..38f86e7b5 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -1,3 +1,16 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. """Hyper Parameter Optimization (HPO) This implementation is based on the work: @@ -90,7 +103,7 @@ def get_samples(self, chains=1, points_per_chain=1, burn_in=2000): return pts -def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf, show_plot=False): +def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf): z = np.array(list(zip(x, y)), dtype=np.dtype([('x', float), ('y', float)])) z = np.sort(z, order='y') n = y.shape[0] @@ -107,13 +120,6 @@ def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf, show_plot=False): samples = lx.get_samples(chains=10, points_per_chain=100) ei = lx(samples) / gx(samples) - if show_plot is True: - import pylab as plt - plt.scatter(samples, lx(samples), color='r') - plt.scatter(samples, gx(samples), color='b') - plt.scatter(samples, ei, color='g') - plt.show() - h = (x.max() - x.min()) / (10 * x.size) # assumes prior of x is uniform -- should change for different priors # d = np.abs(x - samples[ei.argmax()]).min() @@ -138,8 +144,6 @@ def getSample(x, y, dist, minlimit=-np.inf, maxlimit=np.inf): return np.exp(np.random.random() * (np.log(maxlimit) - np.log(minlimit)) + np.log(minlimit)) - else: - logger.error('Unsupported distribution for variable') def fmin(lossfn, @@ -186,6 +190,12 @@ def fmin(lossfn, Best hyperparameter setting found """ + for s in space: + if (space[s]['dist'] is not 'uniform' and + space[s]['dist'] is not 'loguniform'): + logger.error('Unsupported distribution for variable') + raise TypeError('Unknown distribution type for variable') + if (len(trials) > init_random_evals): init_random_evals = 0 diff --git a/brainiak/hyperparamopt/mcmc.py b/brainiak/hyperparamopt/mcmc.py index 0aeb6f5cb..aecee9735 100644 --- a/brainiak/hyperparamopt/mcmc.py +++ b/brainiak/hyperparamopt/mcmc.py @@ -1,3 +1,16 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. """Metropolis-Hasting Random number generator This implementation provides random samples from a user-given diff --git a/brainiak/hyperparamopt/norm.pyx b/brainiak/hyperparamopt/norm.pyx index 6839b1d86..a8a3343ca 100644 --- a/brainiak/hyperparamopt/norm.pyx +++ b/brainiak/hyperparamopt/norm.pyx @@ -1,11 +1,23 @@ -#cython: embedsignature=True - +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. """Cython file for optimizing GMM likelihood computation """ # Authors: Narayanan Sundaram (Intel Labs) +#cython: embedsignature=True import numpy cimport numpy as np diff --git a/examples/hpo_example.py b/examples/hpo_example.py index f4a0af759..6cda63e73 100644 --- a/examples/hpo_example.py +++ b/examples/hpo_example.py @@ -1,3 +1,16 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. import numpy as np import matplotlib.pyplot as plt import brainiak.hyperparamopt.hpo as hpo diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py index 2fc7c204c..8672277e0 100644 --- a/tests/hyperparamopt/test_hpo.py +++ b/tests/hyperparamopt/test_hpo.py @@ -1,3 +1,16 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. import pytest @@ -16,6 +29,7 @@ def test_simple_gmm(): np.testing.assert_array_less(samples, 4.) np.testing.assert_array_less(0., samples) + def test_simple_gmm_weights(): from brainiak.hyperparamopt.hpo import gmm_1d_distribution import numpy as np @@ -46,7 +60,15 @@ def f(args): s = {'x': {'dist': 'uniform', 'lo': -10., 'hi': 10.}} trials = [] - best = fmin(lossfn=f, space=s, maxevals=50, trials=trials, verbose=True) + + #Test fmin and ability to continue adding to trials + best = fmin(lossfn=f, space=s, maxevals=40, trials=trials, verbose=True) + best = fmin(lossfn=f, space=s, maxevals=10, trials=trials, verbose=True) + + assert len(trials) == 50, "HPO continuation trials not working" + + # Test verbose flag + best = fmin(lossfn=f, space=s, maxevals=10, trials=trials, verbose=False) yarray = np.array([tr['loss'] for tr in trials]) np.testing.assert_array_less(yarray, 100.) @@ -56,3 +78,12 @@ def f(args): assert best['loss'] < 100., "HPO out of range" assert np.abs(best['x']) < 10., "HPO out of range" + + #Test unknown distributions + s2 = {'x': {'dist': 'normal', 'mu': 0., 'sigma': 1.}} + trials2 = [] + with pytest.raises(TypeError) as excinfo: + best = fmin(lossfn=f, space=s2, maxevals=40, trials=trials2, verbose=False) + assert "Unknown distribution type for variable" in str(excinfo.value) + + diff --git a/tests/hyperparamopt/test_mcmc.py b/tests/hyperparamopt/test_mcmc.py index 7615938d5..7883b22ec 100644 --- a/tests/hyperparamopt/test_mcmc.py +++ b/tests/hyperparamopt/test_mcmc.py @@ -1,3 +1,16 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. import pytest def test_get_chain(): From 2440bebf2c6aa78a1fb4aada1e5e436baf7ef896 Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Wed, 13 Jul 2016 20:12:31 -0700 Subject: [PATCH 4/9] Removed mcmc, using scipy+numpy samplers for GMM; removed norm.pyx and cython; all scipy.stats continuous distributions supported; more comments added --- brainiak/hyperparamopt/hpo.py | 206 ++++++++++++++++++++++++------- brainiak/hyperparamopt/mcmc.py | 133 -------------------- brainiak/hyperparamopt/norm.pyx | 60 --------- examples/hpo_example.py | 9 +- tests/hyperparamopt/test_hpo.py | 38 +++--- tests/hyperparamopt/test_mcmc.py | 48 ------- 6 files changed, 185 insertions(+), 309 deletions(-) delete mode 100644 brainiak/hyperparamopt/mcmc.py delete mode 100644 brainiak/hyperparamopt/norm.pyx delete mode 100644 tests/hyperparamopt/test_mcmc.py diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index 38f86e7b5..c53f1d010 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -29,17 +29,44 @@ # Authors: Narayanan Sundaram (Intel Labs) import logging -from .mcmc import get_multichain_samples -from .norm import getgmmpdf +import math import numpy as np from scipy.special import erf +import scipy.stats as st from tqdm import tqdm logger = logging.getLogger(__name__) -def getsigma(x, minlimit=-np.inf, maxlimit=np.inf): +def get_sigma(x, minlimit=-np.inf, maxlimit=np.inf): + """Computes the standard deviations around the points for a 1D + Gaussian mixture model computation. + + We take the distance from the nearest left and right neighbors + for each point, then use the max as the estimate of standard + deviation for the gaussian around that point. + + Arguments + --------- + + x : 1D array + Set of points to create the GMM + + minlimit : double, default : -np.inf + Minimum limit for the distribution + + maxlimit : double, default : np.inf + maximum limit for the distribution + + Returns + ------- + + sigma : 1D array + Array of standard deviations + + """ + z = np.append(x, [minlimit, maxlimit]) sigma = np.ones(x.shape) for i in range(x.size): @@ -52,7 +79,6 @@ def getsigma(x, minlimit=-np.inf, maxlimit=np.inf): sigma[i] = 1.0 return sigma - class gmm_1d_distribution: """GMM 1D distribution. @@ -81,29 +107,133 @@ def __init__(self, x, minlimit=-np.inf, maxlimit=np.inf, weights=1.0): self.N = x.size self.minlimit = minlimit self.maxlimit = maxlimit - self.sigma = getsigma(x, minlimit=minlimit, maxlimit=maxlimit) + self.sigma = get_sigma(x, minlimit=minlimit, maxlimit=maxlimit) self.weights = 2. / (erf((maxlimit - x) / (np.sqrt(2.) * self.sigma)) - erf((minlimit - x) / (np.sqrt(2.) * self.sigma))) * weights - # return self + self.W_sum = np.sum(self.weights) + + + def get_gmm_pdf(self, xt): + """Calculates the 1D GMM likelihood for a single point + + y = \sum_{i=1}^{N} norm_pdf(x, x_i, sigma_i)/(\sum weight_i) + """ + + def my_norm_pdf(x, mu, sigma): + z = (x - mu) / sigma + return (math.exp(-0.5 * z * z) + / (math.sqrt(2. * np.pi) * sigma)) + + y = 0 + if (xt < self.minlimit): + return 0 + if (xt > self.maxlimit): + return 0 + for _x in range(self.points.size): + y += (my_norm_pdf(xt, self.points[_x], self.sigma[_x]) + * self.weights[_x]) / self.W_sum + return y def __call__(self, xt): + """Returns the likelihood of point(s) belonging to the GMM + distribution. + + Arguments + --------- + + x : scalar (or) 1D array of reals + Point(s) at which likelihood needs to be computed + + Returns + ------- + + l : scalar (or) 1D array + Likelihood values at the given point(s) + + """ + if (np.isscalar(xt)): - return getgmmpdf(xt, self.points, self.sigma, self.weights, - self.minlimit, self.maxlimit) + return self.get_gmm_pdf(xt) else: - return np.array([getgmmpdf(t, self.points, self.sigma, - self.weights, self.minlimit, - self.maxlimit) for t in xt]) + return np.array([self.get_gmm_pdf(t) for t in xt]) + + def get_samples(self, n): + """Samples the GMM distribution. + + Arguments + --------- + + n : int + Number of samples needed + + Returns + ------- + + samples : 1D array + Samples from the distribution + + """ + + normalized_w = self.weights / np.sum(self.weights) + get_rand_index = st.rv_discrete(values=(range(self.N), + normalized_w)).rvs(size=n) + samples = np.zeros(n) + k = 0 + j = 0 + while (k < n): + i = get_rand_index[j] + j = j + 1 + if (j == n): + get_rand_index = st.rv_discrete(values=(range(self.N), + normalized_w)).rvs(size=n) + j = 0 + v = np.random.normal(loc=self.points[i], scale=self.sigma[i]) + if (v > self.maxlimit or v < self.minlimit): + continue + else: + samples[k] = v + k = k + 1 + if (k == n): + break + return samples + + +def get_next_sample(x, y, minlimit=-np.inf, maxlimit=np.inf): + """Returns the point that gives the largest Expected improvement (EI) in the + optimization function. + + We use [Bergstra2013] to compute this. This model fits 2 different GMMs - + one for points that have loss values in the bottom 15% and another + for the rest. Then we sample from the former distribution and estimate + EI as the ratio of the likelihoods of the 2 distributions. We pick the + point with the best EI among the samples that is also not very close to + a point we have sampled earlier. + + Arguments + --------- + + x : 1D array + Samples generated from the distribution so far + + y : 1D array + Loss values at the corresponding samples + + minlimit : double, default : -inf + Minimum limit for the distribution + + maxlimit : double, default : +inf + Maximum limit for the distribution + + Returns + ------- - def get_samples(self, chains=1, points_per_chain=1, burn_in=2000): - pts = get_multichain_samples(N=points_per_chain, - p=self, nchains=chains, burn_in=burn_in) - return pts + x_next : double + Next value to use for HPO + """ -def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf): z = np.array(list(zip(x, y)), dtype=np.dtype([('x', float), ('y', float)])) z = np.sort(z, order='y') n = y.shape[0] @@ -117,7 +247,7 @@ def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf): maxlimit=maxlimit, weights=weights) gx = gmm_1d_distribution(gdata['x'], minlimit=minlimit, maxlimit=maxlimit) - samples = lx.get_samples(chains=10, points_per_chain=100) + samples = lx.get_samples(n=1000) ei = lx(samples) / gx(samples) h = (x.max() - x.min()) / (10 * x.size) @@ -135,24 +265,12 @@ def getNextSample(x, y, minlimit=-np.inf, maxlimit=np.inf): return xnext -def getSample(x, y, dist, minlimit=-np.inf, maxlimit=np.inf): - if (dist == 'GMM'): - return getNextSample(x, y, minlimit, maxlimit) - if (dist == 'uniform'): - return np.random.random() * (maxlimit - minlimit) + minlimit - if (dist == 'loguniform'): - return np.exp(np.random.random() - * (np.log(maxlimit) - np.log(minlimit)) - + np.log(minlimit)) - - def fmin(lossfn, space, maxevals, trials, init_random_evals=30, - explore_prob=0.2, - verbose=False): + explore_prob=0.2): """Find the minimum of function through hyper paramter optimization Arguments @@ -180,9 +298,6 @@ def fmin(lossfn, Controls the exploration-vs-exploitation ratio Currently 20% of trails are random samples - verbose : bool, default False - Get information on current point being processed - Returns ------- @@ -191,10 +306,13 @@ def fmin(lossfn, """ for s in space: - if (space[s]['dist'] is not 'uniform' and - space[s]['dist'] is not 'loguniform'): + if (hasattr(space[s]['dist'], 'rvs') is False): logger.error('Unsupported distribution for variable') raise TypeError('Unknown distribution type for variable') + if ('lo' not in space[s]): + space[s]['lo'] = -np.inf + if ('hi' not in space[s]): + space[s]['hi'] = np.inf if (len(trials) > init_random_evals): init_random_evals = 0 @@ -210,14 +328,15 @@ def fmin(lossfn, yarray = np.array([tr['loss'] for tr in trials]) for s in space: sarray = np.array([tr[s] for tr in trials]) - dist = 'GMM' if (search_algo == 'Exploit') else space[s]['dist'] - sdict[s] = getSample(sarray, yarray, dist, - minlimit=space[s]['lo'], - maxlimit=space[s]['hi']) + if (search_algo == 'Exploit'): + sdict[s] = get_next_sample(sarray, yarray, + minlimit=space[s]['lo'], + maxlimit=space[s]['hi']) + else: + sdict[s] = space[s]['dist'].rvs() - if (verbose): - logger.info(search_algo) - logger.info('Next point ', t, ' = ', sdict) + logger.debug(search_algo) + logger.info('Next point ', t, ' = ', sdict) y = lossfn(sdict) sdict['loss'] = y @@ -226,6 +345,5 @@ def fmin(lossfn, yarray = np.array([tr['loss'] for tr in trials]) yargmin = yarray.argmin() - if (verbose): - logger.info('Best point so far = ', trials[yargmin]) + logger.info('Best point so far = ', trials[yargmin]) return trials[yargmin] diff --git a/brainiak/hyperparamopt/mcmc.py b/brainiak/hyperparamopt/mcmc.py deleted file mode 100644 index aecee9735..000000000 --- a/brainiak/hyperparamopt/mcmc.py +++ /dev/null @@ -1,133 +0,0 @@ -# Copyright 2016 Intel Corporation -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -"""Metropolis-Hasting Random number generator - -This implementation provides random samples from a user-given -probability density function through the Metropolis-Hasting algorithm. - -""" - -# Authors: Narayanan Sundaram (Intel Labs) - -import numpy as np -import scipy.stats as st -import logging - - -def candidate(x): - """Generates candidate around point x - - Returns - ------- - - sample from q(x*|x) - Unit Normal distribution around x - """ - - return np.random.standard_normal() + x - - -def candidate_dist(x, xp): # return value of q(xp| x) - return st.norm.pdf(x - xp) - - -def check_accept(xcurr, xprop, p): - return min(1.0, p(xprop) / p(xcurr)) - - -def get_next(x, p, n): - xnext = np.zeros(n) - xnext[-1] = x - pxcurr = p(x) - for i in range(n): - xp = candidate(xnext[i - 1]) - pxp = p(xp) - if (np.random.random() < pxp / pxcurr): - xnext[i] = xp - pxcurr = pxp - else: - xnext[i] = xnext[i - 1] - return xnext - - -def get_chain(N, p, burn_in=2000): - """Get a sequence of numbers sampled from a single chain - of MCMC (Metropolis-Hastings) sampler - - Arguments - --------- - - N : int - Number of samples required - - p : function that returns a pdf value at any real number - Distribution that needs to be sampled - - burn_in : int, default 2000 - Number of burn-in (discarded) samples - - Returns - ------- - - samples : 1D array, shape [N] - Samples generated from MCMC sampler (should resemble samples from p(x)) - """ - - x = np.ones(N) - x0 = np.random.standard_normal() * 100 - while(p(x0) <= np.finfo(np.double).eps * 10): - x0 = np.random.standard_normal() * 100 - if (p(x0) <= 0): - logging.error('Markov chain failed to initialize properly \ - - Values probably very far from origin') - - # burn in iterations - x0 = get_next(x0, p, burn_in)[-1] - - # actual iterations - x = get_next(x0, p, N) - return x - - -def get_multichain_samples(N, p, nchains=3, burn_in=2000): - """Get a sequence of numbers sampled from multiple chains - of MCMC (Metropolis-Hastings) sampler - - Arguments - --------- - - N : int - Number of samples per chain required - - p : function that returns a pdf value at any real number - Distribution that needs to be sampled - - nchains : int, default 3 - Number of independent MCMC chains to sample - - burn_in : int, default 2000 - Number of burn-in (discarded) samples - - Returns - ------- - - samples : 1D array, shape [nchains*N] - Samples generated from MCMC sampler (should resemble samples from p(x)) - """ - - pts = np.zeros(nchains * N) - for c in range(nchains): - xp = get_chain(N * 3, p, burn_in=burn_in) - pts[c * N: (c + 1) * N] =\ - xp[np.random.choice(N * 3, N)] # pick N points at random - return pts diff --git a/brainiak/hyperparamopt/norm.pyx b/brainiak/hyperparamopt/norm.pyx deleted file mode 100644 index a8a3343ca..000000000 --- a/brainiak/hyperparamopt/norm.pyx +++ /dev/null @@ -1,60 +0,0 @@ -# Copyright 2016 Intel Corporation -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -"""Cython file for optimizing GMM likelihood computation - -""" - -# Authors: Narayanan Sundaram (Intel Labs) - -#cython: embedsignature=True -import numpy -cimport numpy as np - -cdef extern from "math.h": - double exp(double x) - double sqrt(double x) - - -pi = numpy.pi - - -cpdef double norm_pdf(double x, double mu, double sigma): - """Calculate Gaussian pdf - - Given x, returns exp(-0.5*z*z)/(sigma*sqrt(2.*pi)) where - z = (x-mu)/sigma - """ - - cdef double z - z = (x-mu)/sigma - return exp(-0.5*z*z)/sqrt(2.0*pi)/sigma - - -cpdef double getgmmpdf(double xt,np.ndarray[np.float64_t, ndim=1] x, np.ndarray[np.float64_t, ndim=1] sigma, np.ndarray[np.float64_t, ndim=1] weights, double minlimit, double maxlimit): - """Calculates the 1D GMM likelihood - - y = \sum_{i=1}^{N} norm_pdf(x, x_i, sigma_i)/(\sum weight_i) - """ - - cdef double y - cdef double w - w = sum(weights) - y = 0 - if (xt < minlimit): - return 0 - if (xt > maxlimit): - return 0 - for _x in range(x.size): - y += norm_pdf(xt, x[_x], sigma[_x])*weights[_x]/w - return y diff --git a/examples/hpo_example.py b/examples/hpo_example.py index 6cda63e73..de795383d 100644 --- a/examples/hpo_example.py +++ b/examples/hpo_example.py @@ -11,9 +11,10 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. +import brainiak.hyperparamopt.hpo as hpo import numpy as np +import scipy.stats as st import matplotlib.pyplot as plt -import brainiak.hyperparamopt.hpo as hpo def branin(x1, x2): a = 1.0 @@ -34,12 +35,12 @@ def g(args): x2lo = 0 x2hi = 15 -space = {'x1':{'dist':'uniform', 'lo':x1lo, 'hi':x1hi}, - 'x2':{'dist':'uniform', 'lo':x2lo, 'hi':x2hi}} +space = {'x1':{'dist': st.uniform(x1lo, x1hi-x1lo), 'lo':x1lo, 'hi':x1hi}, + 'x2':{'dist': st.uniform(x2lo, x2hi-x2lo), 'lo':x2lo, 'hi':x2hi}} trials = [] n_hpo_samples = 100 -best = hpo.fmin(lossfn=g, space=space, maxevals=n_hpo_samples, trials=trials, verbose=False) +best = hpo.fmin(lossfn=g, space=space, maxevals=n_hpo_samples, trials=trials) print("Best obtained through HPO (", n_hpo_samples, " samples) = ", best['x1'], best['x2'], "; min value = ", best['loss']) diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py index 8672277e0..da316b030 100644 --- a/tests/hyperparamopt/test_hpo.py +++ b/tests/hyperparamopt/test_hpo.py @@ -12,12 +12,12 @@ # See the License for the specific language governing permissions and # limitations under the License. import pytest +import numpy as np +import scipy.stats as st +from brainiak.hyperparamopt.hpo import gmm_1d_distribution, fmin def test_simple_gmm(): - from brainiak.hyperparamopt.hpo import gmm_1d_distribution - import numpy as np - x = np.array([1., 1., 2., 3., 1.]) d = gmm_1d_distribution(x, minlimit=0., maxlimit=4.) assert d(1.1) > d(3.5), "GMM distribution not behaving correctly" @@ -25,15 +25,12 @@ def test_simple_gmm(): assert d(-1.0) == 0, "GMM distribution out of bounds error" assert d(9.0) == 0, "GMM distribution out of bounds error" - samples = d.get_samples(chains=2, points_per_chain=10, burn_in=50) + samples = d.get_samples(n=25) np.testing.assert_array_less(samples, 4.) np.testing.assert_array_less(0., samples) def test_simple_gmm_weights(): - from brainiak.hyperparamopt.hpo import gmm_1d_distribution - import numpy as np - x = np.array([1., 1., 2., 3., 1., 3.]) d = gmm_1d_distribution(x) @@ -42,33 +39,31 @@ def test_simple_gmm_weights(): d2 = gmm_1d_distribution(x2, weights=w) y2 = d2(np.array([1.1, 2.0])) - assert(d2(1.1) == y2[0], - "GMM distribution array & scalar results don't match") - assert(np.abs(d(1.1) - d2(1.1)) < 1e-5, - "GMM distribution weights not handled correctly") - assert(np.abs(d(2.0) - d2(2.0)) < 1e-5, - "GMM distribution weights not handled correctly") + assert d2(1.1) == y2[0],\ + "GMM distribution array & scalar results don't match" + assert np.abs(d(1.1) - d2(1.1)) < 1e-5,\ + "GMM distribution weights not handled correctly" + assert np.abs(d(2.0) - d2(2.0)) < 1e-5,\ + "GMM distribution weights not handled correctly" def test_simple_hpo(): - from brainiak.hyperparamopt.hpo import fmin - import numpy as np def f(args): x = args['x'] return x*x - s = {'x': {'dist': 'uniform', 'lo': -10., 'hi': 10.}} + s = {'x': {'dist': st.uniform(loc=-10., scale=20), 'lo': -10., 'hi': 10.}} trials = [] #Test fmin and ability to continue adding to trials - best = fmin(lossfn=f, space=s, maxevals=40, trials=trials, verbose=True) - best = fmin(lossfn=f, space=s, maxevals=10, trials=trials, verbose=True) + best = fmin(lossfn=f, space=s, maxevals=40, trials=trials) + best = fmin(lossfn=f, space=s, maxevals=10, trials=trials) assert len(trials) == 50, "HPO continuation trials not working" # Test verbose flag - best = fmin(lossfn=f, space=s, maxevals=10, trials=trials, verbose=False) + best = fmin(lossfn=f, space=s, maxevals=10, trials=trials) yarray = np.array([tr['loss'] for tr in trials]) np.testing.assert_array_less(yarray, 100.) @@ -83,7 +78,10 @@ def f(args): s2 = {'x': {'dist': 'normal', 'mu': 0., 'sigma': 1.}} trials2 = [] with pytest.raises(TypeError) as excinfo: - best = fmin(lossfn=f, space=s2, maxevals=40, trials=trials2, verbose=False) + best2 = fmin(lossfn=f, space=s2, maxevals=40, trials=trials2) assert "Unknown distribution type for variable" in str(excinfo.value) + s3 = {'x': {'dist': st.norm(loc=0., scale=1.)}} + trials3 = [] + best3 = fmin(lossfn=f, space=s3, maxevals=40, trials=trials3) diff --git a/tests/hyperparamopt/test_mcmc.py b/tests/hyperparamopt/test_mcmc.py deleted file mode 100644 index 7883b22ec..000000000 --- a/tests/hyperparamopt/test_mcmc.py +++ /dev/null @@ -1,48 +0,0 @@ -# Copyright 2016 Intel Corporation -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -import pytest - -def test_get_chain(): - from brainiak.hyperparamopt.mcmc import get_multichain_samples - import numpy as np - import numpy.testing as npt - import scipy.stats as st - - def normal(mean, std): - def f(x): - return st.norm.pdf(x, loc=mean, scale=std) - return f - - for mean,std in [(0.,1.), (2.,4.), (-5., 3.)]: - p = normal(mean, std) - samples = get_multichain_samples(1000, p, nchains=5) - assert(np.abs(mean - np.mean(samples)) <= 1.) - assert(np.abs(std - np.std(samples)) <= 1.) - # assert(st.skewtest(samples).pvalue >= 0.05) - # assert(st.kurtosistest(samples).pvalue >= 0.05) - - - """ - #import pylab as plt - #plt.hist(samples, 100) - #plt.show() - # Anderson-Darling test - A2, criticalvalues, significancelevel = st.anderson(samples, 'norm') - print(criticalvalues, significancelevel, A2) - - # critical values at [15, 10, 5, 2.5, 1] - for i in range(len(significancelevel)): - if (significancelevel[i] == 5.): # at 5% significance level - assert(A2 <= criticalvalues[i]) - """ From b93e3a709788f44a08d6f01759f234222e99f128 Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Wed, 13 Jul 2016 20:17:17 -0700 Subject: [PATCH 5/9] formatting fixes --- brainiak/hyperparamopt/hpo.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index c53f1d010..7f1deafae 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -79,6 +79,7 @@ def get_sigma(x, minlimit=-np.inf, maxlimit=np.inf): sigma[i] = 1.0 return sigma + class gmm_1d_distribution: """GMM 1D distribution. @@ -114,7 +115,6 @@ def __init__(self, x, minlimit=-np.inf, maxlimit=np.inf, weights=1.0): / (np.sqrt(2.) * self.sigma))) * weights self.W_sum = np.sum(self.weights) - def get_gmm_pdf(self, xt): """Calculates the 1D GMM likelihood for a single point @@ -122,9 +122,9 @@ def get_gmm_pdf(self, xt): """ def my_norm_pdf(x, mu, sigma): - z = (x - mu) / sigma - return (math.exp(-0.5 * z * z) - / (math.sqrt(2. * np.pi) * sigma)) + z = (x - mu) / sigma + return (math.exp(-0.5 * z * z) + / (math.sqrt(2. * np.pi) * sigma)) y = 0 if (xt < self.minlimit): @@ -133,7 +133,7 @@ def my_norm_pdf(x, mu, sigma): return 0 for _x in range(self.points.size): y += (my_norm_pdf(xt, self.points[_x], self.sigma[_x]) - * self.weights[_x]) / self.W_sum + * self.weights[_x]) / self.W_sum return y def __call__(self, xt): @@ -330,8 +330,8 @@ def fmin(lossfn, sarray = np.array([tr[s] for tr in trials]) if (search_algo == 'Exploit'): sdict[s] = get_next_sample(sarray, yarray, - minlimit=space[s]['lo'], - maxlimit=space[s]['hi']) + minlimit=space[s]['lo'], + maxlimit=space[s]['hi']) else: sdict[s] = space[s]['dist'].rvs() From 6525f36089b9cdaad9767a0cc5b7aa8e8b9bf9bf Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Wed, 13 Jul 2016 20:44:57 -0700 Subject: [PATCH 6/9] removed tqdm; added comments to the example --- brainiak/hyperparamopt/hpo.py | 3 +- examples/{ => hyperparamopt}/hpo_example.py | 70 +++++++++++++++++++-- requirements-dev.txt | 1 - 3 files changed, 65 insertions(+), 9 deletions(-) rename examples/{ => hyperparamopt}/hpo_example.py (51%) diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index 7f1deafae..11b665a76 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -33,7 +33,6 @@ import numpy as np from scipy.special import erf import scipy.stats as st -from tqdm import tqdm logger = logging.getLogger(__name__) @@ -317,7 +316,7 @@ def fmin(lossfn, if (len(trials) > init_random_evals): init_random_evals = 0 - for t in tqdm(range(maxevals)): + for t in range(maxevals): sdict = {} if (t >= init_random_evals and np.random.random() > explore_prob): diff --git a/examples/hpo_example.py b/examples/hyperparamopt/hpo_example.py similarity index 51% rename from examples/hpo_example.py rename to examples/hyperparamopt/hpo_example.py index de795383d..8f508a651 100644 --- a/examples/hpo_example.py +++ b/examples/hyperparamopt/hpo_example.py @@ -11,11 +11,21 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. + +"""Example for using hyperparameter optimization (hpo) package. + +In this example, we will try to optimize a function of +2 variables (branin) using both hpo and grid search. + +""" + import brainiak.hyperparamopt.hpo as hpo import numpy as np import scipy.stats as st import matplotlib.pyplot as plt +# Branin is the function we want to minimize. +# It is a function of 2 variables def branin(x1, x2): a = 1.0 b = 5.1/(4*np.pi*np.pi) @@ -25,48 +35,96 @@ def branin(x1, x2): t = 1.0/(8*np.pi) return a*((x2 - b*x1*x1 + c*x1 - r)**2) + s*(1-t)*np.cos(x1) + s + 5*x1 -def g(args): +# This is a wrapper around branin that takes in a dictionary +def branin_wrapper(args): x1 = args['x1'] x2 = args['x2'] return branin(x1,x2) +# Define ranges for the two variables x1lo = -5 x1hi = 10 x2lo = 0 x2hi = 15 +############################## +# Optimization through hpo +############################## + +# Define a space for hpo to use +# The space needs to define +# 1. Name of the variables +# 2. Default samplers for the variables (use scipy.stats objects) +# 3. lo and hi ranges for the variables (will use -inf, inf if not specified) space = {'x1':{'dist': st.uniform(x1lo, x1hi-x1lo), 'lo':x1lo, 'hi':x1hi}, 'x2':{'dist': st.uniform(x2lo, x2hi-x2lo), 'lo':x2lo, 'hi':x2hi}} + +# The trials object is just a list that stores the samples generated and the +# corresponding function values at those sample points. trials = [] + +# Maximum number of samples that will be generated. +# This is the maximum number of function evaluations that will be performed. n_hpo_samples = 100 -best = hpo.fmin(lossfn=g, space=space, maxevals=n_hpo_samples, trials=trials) +# Call the fmin function that does the optimization. +# The function to be optimized should take in a dictionary. You will probably +# need to wrap your function to do this (see branin() and branin_wrapper()). +# You can pass in a non-empty trials object as well e.g. from a previous +# fmin run. We just append to the trials object and will use existing data +# in our optimization. +print("Starting optimization through hpo") +best = hpo.fmin(lossfn=branin_wrapper, space=space, + maxevals=n_hpo_samples, trials=trials) + +# Print out the best value obtained through HPO print("Best obtained through HPO (", n_hpo_samples, " samples) = ", best['x1'], best['x2'], "; min value = ", best['loss']) -nt = 100 +##################################### +# Optimization through grid search +##################################### + +# Divide the space into a uniform grid (meshgrid) +nt = 200 x1t = np.linspace(x1lo, x1hi, nt) x2t = np.linspace(x2lo, x2hi, nt) x1m, x2m = np.meshgrid(x1t, x2t) + +# Calculate the function values along the grid +print("Starting optimization through grid search") z = branin(x1m, x2m) + +# Print out the best value obtained through grid search print("Best obtained through grid search (", nt*nt, " samples) = ", x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], "; min value = ", z.min()) +######## +# Plots +######## + +# Convert trials object data into numpy arrays x1 = np.array([tr['x1'] for tr in trials]) x2 = np.array([tr['x2'] for tr in trials]) y = np.array([tr['loss'] for tr in trials]) +# Plot the function contour using the grid search data h = (z.max()-z.min())/25 plt.contour(x1m, x2m, z, levels=np.linspace(z.min()-h, z.max(), 26)) + +# Mark the points that were sampled through HPO plt.scatter(x1, x2, s=10, color='r', label='HPO Samples') -plt.xlabel('x1') -plt.ylabel('x2') -plt.title('Hyperparamter optimization using HPO') +# Mark the best points obtained through both methods plt.scatter(best['x1'], best['x2'], s=30, color='b', label='Best HPO') plt.scatter(x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], s=30, color='g', label='Best grid search') + +# Labels +plt.xlabel('x1') +plt.ylabel('x2') +plt.title('Hyperparamter optimization using HPO') plt.legend() plt.show() diff --git a/requirements-dev.txt b/requirements-dev.txt index af02c3f87..bd99b5fb3 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -7,4 +7,3 @@ pytest-cython restructuredtext-lint sphinx sphinx_rtd_theme -tqdm From 1511e31589b14d38d96240e08b0c028266391513 Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Thu, 14 Jul 2016 12:38:50 -0700 Subject: [PATCH 7/9] Updated according to comments --- brainiak/hyperparamopt/hpo.py | 87 +++++++++++++------------ examples/hyperparamopt/hpo_example.py | 24 +++---- examples/hyperparamopt/requirements.txt | 1 + tests/hyperparamopt/test_hpo.py | 12 ++-- 4 files changed, 63 insertions(+), 61 deletions(-) create mode 100644 examples/hyperparamopt/requirements.txt diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index 11b665a76..6f401fc7f 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -38,7 +38,7 @@ logger = logging.getLogger(__name__) -def get_sigma(x, minlimit=-np.inf, maxlimit=np.inf): +def get_sigma(x, min_limit=-np.inf, max_limit=np.inf): """Computes the standard deviations around the points for a 1D Gaussian mixture model computation. @@ -52,10 +52,10 @@ def get_sigma(x, minlimit=-np.inf, maxlimit=np.inf): x : 1D array Set of points to create the GMM - minlimit : double, default : -np.inf + min_limit : double, default : -np.inf Minimum limit for the distribution - maxlimit : double, default : np.inf + max_limit : double, default : np.inf maximum limit for the distribution Returns @@ -66,7 +66,7 @@ def get_sigma(x, minlimit=-np.inf, maxlimit=np.inf): """ - z = np.append(x, [minlimit, maxlimit]) + z = np.append(x, [min_limit, max_limit]) sigma = np.ones(x.shape) for i in range(x.size): xleft = z[np.argmin([(x[i] - k) if k < x[i] else np.inf for k in z])] @@ -92,50 +92,50 @@ class gmm_1d_distribution: x : 1D array Set of points to create the GMM - minlimit : double, default : -inf + min_limit : double, default : -inf Minimum limit for the distribution - maxlimit : double, default : +inf + max_limit : double, default : +inf Maximum limit for the distribution weights : double scalar or 1D array with same size as x, default 1.0 Used to weight the points non-uniformly if required """ - def __init__(self, x, minlimit=-np.inf, maxlimit=np.inf, weights=1.0): + def __init__(self, x, min_limit=-np.inf, max_limit=np.inf, weights=1.0): self.points = x self.N = x.size - self.minlimit = minlimit - self.maxlimit = maxlimit - self.sigma = get_sigma(x, minlimit=minlimit, maxlimit=maxlimit) - self.weights = 2. / (erf((maxlimit - x) + self.min_limit = min_limit + self.max_limit = max_limit + self.sigma = get_sigma(x, min_limit=min_limit, max_limit=max_limit) + self.weights = 2. / (erf((max_limit - x) / (np.sqrt(2.) * self.sigma)) - - erf((minlimit - x) + - erf((min_limit - x) / (np.sqrt(2.) * self.sigma))) * weights self.W_sum = np.sum(self.weights) - def get_gmm_pdf(self, xt): + def get_gmm_pdf(self, x): """Calculates the 1D GMM likelihood for a single point y = \sum_{i=1}^{N} norm_pdf(x, x_i, sigma_i)/(\sum weight_i) """ - def my_norm_pdf(x, mu, sigma): - z = (x - mu) / sigma + def my_norm_pdf(xt, mu, sigma): + z = (xt - mu) / sigma return (math.exp(-0.5 * z * z) / (math.sqrt(2. * np.pi) * sigma)) y = 0 - if (xt < self.minlimit): + if (x < self.min_limit): return 0 - if (xt > self.maxlimit): + if (x > self.max_limit): return 0 for _x in range(self.points.size): - y += (my_norm_pdf(xt, self.points[_x], self.sigma[_x]) + y += (my_norm_pdf(x, self.points[_x], self.sigma[_x]) * self.weights[_x]) / self.W_sum return y - def __call__(self, xt): + def __call__(self, x): """Returns the likelihood of point(s) belonging to the GMM distribution. @@ -153,10 +153,10 @@ def __call__(self, xt): """ - if (np.isscalar(xt)): - return self.get_gmm_pdf(xt) + if np.isscalar(x): + return self.get_gmm_pdf(x) else: - return np.array([self.get_gmm_pdf(t) for t in xt]) + return np.array([self.get_gmm_pdf(t) for t in x]) def get_samples(self, n): """Samples the GMM distribution. @@ -189,7 +189,7 @@ def get_samples(self, n): normalized_w)).rvs(size=n) j = 0 v = np.random.normal(loc=self.points[i], scale=self.sigma[i]) - if (v > self.maxlimit or v < self.minlimit): + if (v > self.max_limit or v < self.min_limit): continue else: samples[k] = v @@ -199,7 +199,7 @@ def get_samples(self, n): return samples -def get_next_sample(x, y, minlimit=-np.inf, maxlimit=np.inf): +def get_next_sample(x, y, min_limit=-np.inf, max_limit=np.inf): """Returns the point that gives the largest Expected improvement (EI) in the optimization function. @@ -219,10 +219,10 @@ def get_next_sample(x, y, minlimit=-np.inf, maxlimit=np.inf): y : 1D array Loss values at the corresponding samples - minlimit : double, default : -inf + min_limit : double, default : -inf Minimum limit for the distribution - maxlimit : double, default : +inf + max_limit : double, default : +inf Maximum limit for the distribution Returns @@ -242,9 +242,10 @@ def get_next_sample(x, y, minlimit=-np.inf, maxlimit=np.inf): lymin = ldata['y'].min() lymax = ldata['y'].max() weights = (lymax - ldata['y']) / (lymax - lymin) - lx = gmm_1d_distribution(ldata['x'], minlimit=minlimit, - maxlimit=maxlimit, weights=weights) - gx = gmm_1d_distribution(gdata['x'], minlimit=minlimit, maxlimit=maxlimit) + lx = gmm_1d_distribution(ldata['x'], min_limit=min_limit, + max_limit=max_limit, weights=weights) + gx = gmm_1d_distribution(gdata['x'], min_limit=min_limit, + max_limit=max_limit) samples = lx.get_samples(n=1000) ei = lx(samples) / gx(samples) @@ -264,9 +265,9 @@ def get_next_sample(x, y, minlimit=-np.inf, maxlimit=np.inf): return xnext -def fmin(lossfn, +def fmin(loss_fn, space, - maxevals, + max_evals, trials, init_random_evals=30, explore_prob=0.2): @@ -275,14 +276,14 @@ def fmin(lossfn, Arguments --------- - lossfn : function that takes in a dictionary and returns a real value + loss_fn : function that takes in a dictionary and returns a real value Function to be minimized space : Dictionary specifying the range and distribution of the hyperparamters - maxevals : int - Maximum number of evaluations of lossfn allowed + max_evals : int + Maximum number of evaluations of loss_fn allowed trials : list Holds the output of the optimization trials @@ -305,21 +306,21 @@ def fmin(lossfn, """ for s in space: - if (hasattr(space[s]['dist'], 'rvs') is False): + if not hasattr(space[s]['dist'], 'rvs'): logger.error('Unsupported distribution for variable') raise TypeError('Unknown distribution type for variable') - if ('lo' not in space[s]): + if 'lo' not in space[s]: space[s]['lo'] = -np.inf - if ('hi' not in space[s]): + if 'hi' not in space[s]: space[s]['hi'] = np.inf - if (len(trials) > init_random_evals): + if len(trials) > init_random_evals: init_random_evals = 0 - for t in range(maxevals): + for t in range(max_evals): sdict = {} - if (t >= init_random_evals and np.random.random() > explore_prob): + if t >= init_random_evals and np.random.random() > explore_prob: search_algo = 'Exploit' else: search_algo = 'Explore' @@ -329,15 +330,15 @@ def fmin(lossfn, sarray = np.array([tr[s] for tr in trials]) if (search_algo == 'Exploit'): sdict[s] = get_next_sample(sarray, yarray, - minlimit=space[s]['lo'], - maxlimit=space[s]['hi']) + min_limit=space[s]['lo'], + max_limit=space[s]['hi']) else: sdict[s] = space[s]['dist'].rvs() logger.debug(search_algo) logger.info('Next point ', t, ' = ', sdict) - y = lossfn(sdict) + y = loss_fn(sdict) sdict['loss'] = y trials.append(sdict) diff --git a/examples/hyperparamopt/hpo_example.py b/examples/hyperparamopt/hpo_example.py index 8f508a651..ef13c19e3 100644 --- a/examples/hyperparamopt/hpo_example.py +++ b/examples/hyperparamopt/hpo_example.py @@ -74,8 +74,8 @@ def branin_wrapper(args): # fmin run. We just append to the trials object and will use existing data # in our optimization. print("Starting optimization through hpo") -best = hpo.fmin(lossfn=branin_wrapper, space=space, - maxevals=n_hpo_samples, trials=trials) +best = hpo.fmin(loss_fn=branin_wrapper, space=space, + max_evals=n_hpo_samples, trials=trials) # Print out the best value obtained through HPO print("Best obtained through HPO (", n_hpo_samples, " samples) = ", @@ -86,18 +86,18 @@ def branin_wrapper(args): ##################################### # Divide the space into a uniform grid (meshgrid) -nt = 200 -x1t = np.linspace(x1lo, x1hi, nt) -x2t = np.linspace(x2lo, x2hi, nt) -x1m, x2m = np.meshgrid(x1t, x2t) +n = 200 +x1 = np.linspace(x1lo, x1hi, n) +x2 = np.linspace(x2lo, x2hi, n) +x1_grid, x2_grid = np.meshgrid(x1, x2) # Calculate the function values along the grid print("Starting optimization through grid search") -z = branin(x1m, x2m) +z = branin(x1_grid, x2_grid) # Print out the best value obtained through grid search -print("Best obtained through grid search (", nt*nt, " samples) = ", - x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], +print("Best obtained through grid search (", n*n, " samples) = ", + x1_grid.flatten()[z.argmin()], x2_grid.flatten()[z.argmin()], "; min value = ", z.min()) ######## @@ -111,20 +111,20 @@ def branin_wrapper(args): # Plot the function contour using the grid search data h = (z.max()-z.min())/25 -plt.contour(x1m, x2m, z, levels=np.linspace(z.min()-h, z.max(), 26)) +plt.contour(x1_grid, x2_grid, z, levels=np.linspace(z.min()-h, z.max(), 26)) # Mark the points that were sampled through HPO plt.scatter(x1, x2, s=10, color='r', label='HPO Samples') # Mark the best points obtained through both methods plt.scatter(best['x1'], best['x2'], s=30, color='b', label='Best HPO') -plt.scatter(x1m.flatten()[z.argmin()], x2m.flatten()[z.argmin()], +plt.scatter(x1_grid.flatten()[z.argmin()], x2_grid.flatten()[z.argmin()], s=30, color='g', label='Best grid search') # Labels plt.xlabel('x1') plt.ylabel('x2') -plt.title('Hyperparamter optimization using HPO') +plt.title('Hyperparamter optimization using HPO (Branin function)') plt.legend() plt.show() diff --git a/examples/hyperparamopt/requirements.txt b/examples/hyperparamopt/requirements.txt new file mode 100644 index 000000000..6ccafc3f9 --- /dev/null +++ b/examples/hyperparamopt/requirements.txt @@ -0,0 +1 @@ +matplotlib diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py index da316b030..ea4ac5aae 100644 --- a/tests/hyperparamopt/test_hpo.py +++ b/tests/hyperparamopt/test_hpo.py @@ -19,7 +19,7 @@ def test_simple_gmm(): x = np.array([1., 1., 2., 3., 1.]) - d = gmm_1d_distribution(x, minlimit=0., maxlimit=4.) + d = gmm_1d_distribution(x, min_limit=0., max_limit=4.) assert d(1.1) > d(3.5), "GMM distribution not behaving correctly" assert d(2.0) > d(3.0), "GMM distribution not behaving correctly" assert d(-1.0) == 0, "GMM distribution out of bounds error" @@ -57,13 +57,13 @@ def f(args): trials = [] #Test fmin and ability to continue adding to trials - best = fmin(lossfn=f, space=s, maxevals=40, trials=trials) - best = fmin(lossfn=f, space=s, maxevals=10, trials=trials) + best = fmin(loss_fn=f, space=s, max_evals=40, trials=trials) + best = fmin(loss_fn=f, space=s, max_evals=10, trials=trials) assert len(trials) == 50, "HPO continuation trials not working" # Test verbose flag - best = fmin(lossfn=f, space=s, maxevals=10, trials=trials) + best = fmin(loss_fn=f, space=s, max_evals=10, trials=trials) yarray = np.array([tr['loss'] for tr in trials]) np.testing.assert_array_less(yarray, 100.) @@ -78,10 +78,10 @@ def f(args): s2 = {'x': {'dist': 'normal', 'mu': 0., 'sigma': 1.}} trials2 = [] with pytest.raises(TypeError) as excinfo: - best2 = fmin(lossfn=f, space=s2, maxevals=40, trials=trials2) + best2 = fmin(loss_fn=f, space=s2, max_evals=40, trials=trials2) assert "Unknown distribution type for variable" in str(excinfo.value) s3 = {'x': {'dist': st.norm(loc=0., scale=1.)}} trials3 = [] - best3 = fmin(lossfn=f, space=s3, maxevals=40, trials=trials3) + best3 = fmin(loss_fn=f, space=s3, max_evals=40, trials=trials3) From 54f179edafb5020bc635eaf884e5121021a01c8f Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Thu, 14 Jul 2016 14:08:39 -0700 Subject: [PATCH 8/9] Added more notes on the branin function --- examples/hyperparamopt/hpo_example.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/examples/hyperparamopt/hpo_example.py b/examples/hyperparamopt/hpo_example.py index ef13c19e3..ba3635c2e 100644 --- a/examples/hyperparamopt/hpo_example.py +++ b/examples/hyperparamopt/hpo_example.py @@ -25,7 +25,14 @@ import matplotlib.pyplot as plt # Branin is the function we want to minimize. -# It is a function of 2 variables +# It is a function of 2 variables. +# In the range x1 in [-5, 10] and x2 in [0, 15], +# this function has 2 local minima and 1 global minima. +# Global minima of -16.6 at (-3.7, 13.7). +# This is the modified version (Branin-Hoo) of the standard branin function. +# If you want the standard version (which has 3 global minima), +# you can omit the "+5*x1" term at the end +# For more details, see http://www.sfu.ca/~ssurjano/branin.html def branin(x1, x2): a = 1.0 b = 5.1/(4*np.pi*np.pi) From 78d1b0aac3569d8aa30afd31b735293b09459bca Mon Sep 17 00:00:00 2001 From: Narayanan Sundaram Date: Fri, 15 Jul 2016 15:44:13 -0700 Subject: [PATCH 9/9] updated according to comments --- brainiak/hyperparamopt/hpo.py | 202 ++++++++++++++------------ examples/hyperparamopt/hpo_example.py | 2 +- tests/hyperparamopt/test_hpo.py | 2 +- 3 files changed, 113 insertions(+), 93 deletions(-) diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py index 6f401fc7f..19d13130b 100644 --- a/brainiak/hyperparamopt/hpo.py +++ b/brainiak/hyperparamopt/hpo.py @@ -39,42 +39,44 @@ def get_sigma(x, min_limit=-np.inf, max_limit=np.inf): - """Computes the standard deviations around the points for a 1D - Gaussian mixture model computation. + """Compute the standard deviations around the points for a 1D GMM. We take the distance from the nearest left and right neighbors for each point, then use the max as the estimate of standard - deviation for the gaussian around that point. + deviation for the gaussian mixture around that point. Arguments --------- - x : 1D array - Set of points to create the GMM + Set of points to create the GMM - min_limit : double, default : -np.inf - Minimum limit for the distribution + min_limit : Optional[float], default : -inf + Minimum limit for the distribution - max_limit : double, default : np.inf - maximum limit for the distribution + max_limit : Optional[float], default : inf + maximum limit for the distribution Returns ------- - - sigma : 1D array - Array of standard deviations - + 1D array + Array of standard deviations """ z = np.append(x, [min_limit, max_limit]) sigma = np.ones(x.shape) for i in range(x.size): + # Calculate the nearest left neighbor of x[i] + # Find the minimum of (x[i] - k) for k < x[i] xleft = z[np.argmin([(x[i] - k) if k < x[i] else np.inf for k in z])] + + # Calculate the nearest right neighbor of x[i] + # Find the minimum of (k - x[i]) for k > x[i] xright = z[np.argmin([(k - x[i]) if k > x[i] else np.inf for k in z])] + sigma[i] = max(x[i] - xleft, xright - x[i]) if sigma[i] == np.inf: sigma[i] = min(x[i] - xleft, xright - x[i]) - if (sigma[i] == -np.inf): + if (sigma[i] == -np.inf): # should never happen sigma[i] = 1.0 return sigma @@ -86,20 +88,22 @@ class gmm_1d_distribution: can calculate likelihoods and generate samples from this 1D Gaussian mixture model. - Parameters + Attributes ---------- + points : 1D array + Set of points to create the GMM - x : 1D array - Set of points to create the GMM + N : int + Number of points to create the GMM - min_limit : double, default : -inf - Minimum limit for the distribution + min_limit : Optional[float], default : -inf + Minimum limit for the distribution - max_limit : double, default : +inf - Maximum limit for the distribution + max_limit : Optional[float], default : inf + Maximum limit for the distribution - weights : double scalar or 1D array with same size as x, default 1.0 - Used to weight the points non-uniformly if required + weights : Optional[1D array], default : array of ones + Used to weight the points non-uniformly if required """ def __init__(self, x, min_limit=-np.inf, max_limit=np.inf, weights=1.0): @@ -108,16 +112,27 @@ def __init__(self, x, min_limit=-np.inf, max_limit=np.inf, weights=1.0): self.min_limit = min_limit self.max_limit = max_limit self.sigma = get_sigma(x, min_limit=min_limit, max_limit=max_limit) - self.weights = 2. / (erf((max_limit - x) - / (np.sqrt(2.) * self.sigma)) - - erf((min_limit - x) - / (np.sqrt(2.) * self.sigma))) * weights + self.weights = (2 + / (erf((max_limit - x) / (np.sqrt(2.) * self.sigma)) + - erf((min_limit - x) / (np.sqrt(2.) * self.sigma))) + * weights) self.W_sum = np.sum(self.weights) def get_gmm_pdf(self, x): - """Calculates the 1D GMM likelihood for a single point + """Calculate the GMM likelihood for a single point. + + .. math:: + y = \sum_{i=1}^{N} w_i*normpdf(x, x_i, \sigma_i)/\sum_{i=1}^{N} w_i - y = \sum_{i=1}^{N} norm_pdf(x, x_i, sigma_i)/(\sum weight_i) + Arguments + --------- + x : float + Point at which likelihood needs to be computed + + Returns + ------- + float + Likelihood value at x """ def my_norm_pdf(xt, mu, sigma): @@ -136,21 +151,20 @@ def my_norm_pdf(xt, mu, sigma): return y def __call__(self, x): - """Returns the likelihood of point(s) belonging to the GMM - distribution. + """Return the GMM likelihood for given point(s). + + .. math:: + y = \sum_{i=1}^{N} w_i*normpdf(x, x_i, \sigma_i)/\sum_{i=1}^{N} w_i Arguments --------- - x : scalar (or) 1D array of reals - Point(s) at which likelihood needs to be computed + Point(s) at which likelihood needs to be computed Returns ------- - - l : scalar (or) 1D array - Likelihood values at the given point(s) - + scalar (or) 1D array + Likelihood values at the given point(s) """ if np.isscalar(x): @@ -159,20 +173,17 @@ def __call__(self, x): return np.array([self.get_gmm_pdf(t) for t in x]) def get_samples(self, n): - """Samples the GMM distribution. + """Sample the GMM distribution. Arguments --------- - n : int - Number of samples needed + Number of samples needed Returns ------- - - samples : 1D array - Samples from the distribution - + 1D array + Samples from the distribution """ normalized_w = self.weights / np.sum(self.weights) @@ -200,37 +211,34 @@ def get_samples(self, n): def get_next_sample(x, y, min_limit=-np.inf, max_limit=np.inf): - """Returns the point that gives the largest Expected improvement (EI) in the - optimization function. + """Get the next point to try, given the previous samples. - We use [Bergstra2013] to compute this. This model fits 2 different GMMs - - one for points that have loss values in the bottom 15% and another - for the rest. Then we sample from the former distribution and estimate - EI as the ratio of the likelihoods of the 2 distributions. We pick the - point with the best EI among the samples that is also not very close to - a point we have sampled earlier. + We use [Bergstra2013]_ to compute the point that gives the largest + Expected improvement (EI) in the optimization function. This model fits 2 + different GMMs - one for points that have loss values in the bottom 15% + and another for the rest. Then we sample from the former distribution + and estimate EI as the ratio of the likelihoods of the 2 distributions. + We pick the point with the best EI among the samples that is also not + very close to a point we have sampled earlier. Arguments --------- - x : 1D array - Samples generated from the distribution so far + Samples generated from the distribution so far y : 1D array - Loss values at the corresponding samples + Loss values at the corresponding samples - min_limit : double, default : -inf - Minimum limit for the distribution + min_limit : float, default : -inf + Minimum limit for the distribution - max_limit : double, default : +inf - Maximum limit for the distribution + max_limit : float, default : +inf + Maximum limit for the distribution Returns ------- - - x_next : double - Next value to use for HPO - + float + Next value to use for HPO """ z = np.array(list(zip(x, y)), dtype=np.dtype([('x', float), ('y', float)])) @@ -251,7 +259,8 @@ def get_next_sample(x, y, min_limit=-np.inf, max_limit=np.inf): ei = lx(samples) / gx(samples) h = (x.max() - x.min()) / (10 * x.size) - # assumes prior of x is uniform -- should change for different priors + # TODO + # assumes prior of x is uniform; should ideally change for other priors # d = np.abs(x - samples[ei.argmax()]).min() # CDF(x+d/2) - CDF(x-d/2) < 1/(10*x.size) then reject else accept s = 0 @@ -271,44 +280,55 @@ def fmin(loss_fn, trials, init_random_evals=30, explore_prob=0.2): - """Find the minimum of function through hyper paramter optimization + """Find the minimum of function through hyper parameter optimization. Arguments --------- + loss_fn : ``function(*args) -> float`` + Function that takes in a dictionary and returns a real value. + This is the function to be minimized. - loss_fn : function that takes in a dictionary and returns a real value - Function to be minimized - - space : Dictionary specifying the range and distribution of - the hyperparamters + space : dictionary + Custom dictionary specifying the range and distribution of + the hyperparamters. + E.g. ``space = {'x': {'dist':scipy.stats.uniform(0,1), + 'lo':0, 'hi':1}}`` + for a 1-dimensional space with variable x in range [0,1] max_evals : int - Maximum number of evaluations of loss_fn allowed + Maximum number of evaluations of loss_fn allowed trials : list - Holds the output of the optimization trials - Need not be empty to begin with, new trials are appended - at the end + Holds the output of the optimization trials. + Need not be empty to begin with, new trials are appended + at the end. - init_random_evals : int, default 30 - Number of random trials to initialize the - optimization + init_random_evals : Optional[int], default 30 + Number of random trials to initialize the + optimization. - explore_prob : double in [0, 1], default 0.2 - Controls the exploration-vs-exploitation ratio - Currently 20% of trails are random samples + explore_prob : Optional[float], default 0.2 + Controls the exploration-vs-exploitation ratio. Value should + be in [0,1]. By default, 20% of trails are random samples. Returns ------- - - best : trial entry (dictionary of hyperparameters) - Best hyperparameter setting found + trial entry (dictionary of hyperparameters) + Best hyperparameter setting found. + E.g. {'x': 5.6, 'loss' : 0.5} where x is the best hyparameter + value found and loss is the value of the function for the + best hyperparameter value(s). + + Raises + ------ + ValueError + If the distribution specified in space does not support a ``rvs()`` + method to generate random numbers, a ValueError is raised. """ for s in space: if not hasattr(space[s]['dist'], 'rvs'): - logger.error('Unsupported distribution for variable') - raise TypeError('Unknown distribution type for variable') + raise ValueError('Unknown distribution type for variable') if 'lo' not in space[s]: space[s]['lo'] = -np.inf if 'hi' not in space[s]: @@ -321,21 +341,21 @@ def fmin(loss_fn, sdict = {} if t >= init_random_evals and np.random.random() > explore_prob: - search_algo = 'Exploit' + use_random_sampling = False else: - search_algo = 'Explore' + use_random_sampling = True yarray = np.array([tr['loss'] for tr in trials]) for s in space: sarray = np.array([tr[s] for tr in trials]) - if (search_algo == 'Exploit'): + if use_random_sampling: + sdict[s] = space[s]['dist'].rvs() + else: sdict[s] = get_next_sample(sarray, yarray, min_limit=space[s]['lo'], max_limit=space[s]['hi']) - else: - sdict[s] = space[s]['dist'].rvs() - logger.debug(search_algo) + logger.debug('Explore' if use_random_sampling else 'Exploit') logger.info('Next point ', t, ' = ', sdict) y = loss_fn(sdict) diff --git a/examples/hyperparamopt/hpo_example.py b/examples/hyperparamopt/hpo_example.py index ba3635c2e..e5fc030f8 100644 --- a/examples/hyperparamopt/hpo_example.py +++ b/examples/hyperparamopt/hpo_example.py @@ -131,7 +131,7 @@ def branin_wrapper(args): # Labels plt.xlabel('x1') plt.ylabel('x2') -plt.title('Hyperparamter optimization using HPO (Branin function)') +plt.title('Hyperparameter optimization using HPO (Branin function)') plt.legend() plt.show() diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py index ea4ac5aae..64b9b7eab 100644 --- a/tests/hyperparamopt/test_hpo.py +++ b/tests/hyperparamopt/test_hpo.py @@ -77,7 +77,7 @@ def f(args): #Test unknown distributions s2 = {'x': {'dist': 'normal', 'mu': 0., 'sigma': 1.}} trials2 = [] - with pytest.raises(TypeError) as excinfo: + with pytest.raises(ValueError) as excinfo: best2 = fmin(loss_fn=f, space=s2, max_evals=40, trials=trials2) assert "Unknown distribution type for variable" in str(excinfo.value)