From 5321bfc8c9631cb8593d16da7c617435003d82aa Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 3 May 2013 14:22:18 +0100 Subject: [PATCH 1/3] some minor example modifications and cgd adjustments --- GPy/examples/dimensionality_reduction.py | 14 +++++------ GPy/inference/conjugate_gradient_descent.py | 28 +++++++++++---------- GPy/models/Bayesian_GPLVM.py | 16 ++++++------ GPy/testing/cgd_tests.py | 6 ++--- 4 files changed, 33 insertions(+), 31 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 03b041f1..e4fcc234 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -176,13 +176,12 @@ def bgplvm_simulation_matlab_compare(): Y = sim_data['Y'] S = sim_data['S'] mu = sim_data['mu'] - M, [_, Q] = 30, mu.shape - Q = 2 + M, [_, Q] = 20, mu.shape from GPy.models import mrd from GPy import kern reload(mrd); reload(kern) - #k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) + # k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, # X=mu, @@ -191,14 +190,15 @@ def bgplvm_simulation_matlab_compare(): m.ensure_default_constraints() m.auto_scale_factor = True m['noise'] = Y.var() / 100. + m['linear_variance'] = .01 - lscstr = '{}'.format(k.parts[0].name) +# lscstr = '{}'.format(k.parts[0].name) # m[lscstr] = .01 - m.unconstrain(lscstr); m.constrain_fixed(lscstr, 10) +# m.unconstrain(lscstr); m.constrain_fixed(lscstr, 10) - lscstr = 'X_variance' +# lscstr = 'X_variance' # m[lscstr] = .01 - m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1) +# m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1) # cstr = 'white' # m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.) diff --git a/GPy/inference/conjugate_gradient_descent.py b/GPy/inference/conjugate_gradient_descent.py index 93dac6df..c88249c3 100644 --- a/GPy/inference/conjugate_gradient_descent.py +++ b/GPy/inference/conjugate_gradient_descent.py @@ -166,25 +166,26 @@ class Async_Optimize(object): except Empty: pass - def fmin_async(self, f, df, x0, callback, update_rule=FletcherReeves, + def opt_async(self, f, df, x0, callback, update_rule=FletcherReeves, messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6, report_every=10, *args, **kwargs): self.runsignal.set() outqueue = Queue() + c = None if callback: self.callback = callback - c = Thread(target=self.async_callback_collect, args=(outqueue,)) - c.start() + c = Thread(target=self.async_callback_collect, args=(outqueue,)) + c.start() p = _CGDAsync(f, df, x0, update_rule, self.runsignal, self.SENTINEL, report_every=report_every, messages=messages, maxiter=maxiter, max_f_eval=max_f_eval, gtol=gtol, outqueue=outqueue, *args, **kwargs) p.run() return p, c - def fmin(self, f, df, x0, callback=None, update_rule=FletcherReeves, + def opt(self, f, df, x0, callback=None, update_rule=FletcherReeves, messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6, report_every=10, *args, **kwargs): - p, c = self.fmin_async(f, df, x0, callback, update_rule, messages, + p, c = self.opt_async(f, df, x0, callback, update_rule, messages, maxiter, max_f_eval, gtol, report_every, *args, **kwargs) while self.runsignal.is_set(): @@ -195,7 +196,8 @@ class Async_Optimize(object): # print "^C" self.runsignal.clear() p.join() - if c.is_alive(): + c.join() + if c and c.is_alive(): print "WARNING: callback still running, optimisation done!" return p.result @@ -208,11 +210,11 @@ class CGD(Async_Optimize): if df returns tuple (grad, natgrad) it will optimize according to natural gradient rules ''' - name = "Conjugate Gradient Descent" + opt_name = "Conjugate Gradient Descent" - def fmin_async(self, *a, **kw): + def opt_async(self, *a, **kw): """ - fmin_async(self, f, df, x0, callback, update_rule=FletcherReeves, + opt_async(self, f, df, x0, callback, update_rule=FletcherReeves, messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6, report_every=10, *args, **kwargs) @@ -240,11 +242,11 @@ class CGD(Async_Optimize): at end of optimization! """ - return super(CGD, self).fmin_async(*a, **kw) + return super(CGD, self).opt_async(*a, **kw) - def fmin(self, *a, **kw): + def opt(self, *a, **kw): """ - fmin(self, f, df, x0, callback=None, update_rule=FletcherReeves, + opt(self, f, df, x0, callback=None, update_rule=FletcherReeves, messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6, report_every=10, *args, **kwargs) @@ -267,5 +269,5 @@ class CGD(Async_Optimize): at end of optimization """ - return super(CGD, self).fmin(*a, **kw) + return super(CGD, self).opt(*a, **kw) diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index b824bdfe..099a17ea 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -259,28 +259,28 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): ax2.text(.5, .5, r"$\mathbf{X}$", alpha=.5, transform=ax2.transAxes, ha='center', va='center') figs[-1].canvas.draw() - figs[-1].tight_layout(rect=(0, 0, 1, .9)) + figs[-1].tight_layout(rect=(0, 0, 1, .86)) # ax3 = pylab.subplot2grid(splotshape, (3, 0), 2, 4, sharex=ax2) figs.append(pylab.figure("BGPLVM DEBUG S", figsize=(12, 4))) ax3 = self._debug_get_axis(figs) ax3.text(.5, .5, r"$\mathbf{S}$", alpha=.5, transform=ax3.transAxes, ha='center', va='center') figs[-1].canvas.draw() - figs[-1].tight_layout(rect=(0, 0, 1, .9)) + figs[-1].tight_layout(rect=(0, 0, 1, .86)) # ax4 = pylab.subplot2grid(splotshape, (5, 0), 2, 2) figs.append(pylab.figure("BGPLVM DEBUG Z", figsize=(6, 4))) ax4 = self._debug_get_axis(figs) ax4.text(.5, .5, r"$\mathbf{Z}$", alpha=.5, transform=ax4.transAxes, ha='center', va='center') figs[-1].canvas.draw() - figs[-1].tight_layout(rect=(0, 0, 1, .9)) + figs[-1].tight_layout(rect=(0, 0, 1, .86)) # ax5 = pylab.subplot2grid(splotshape, (5, 2), 2, 2) figs.append(pylab.figure("BGPLVM DEBUG theta", figsize=(6, 4))) ax5 = self._debug_get_axis(figs) ax5.text(.5, .5, r"${\theta}$", alpha=.5, transform=ax5.transAxes, ha='center', va='center') figs[-1].canvas.draw() - figs[-1].tight_layout(rect=(0, 0, 1, .9)) + figs[-1].tight_layout(rect=(.15, 0, 1, .86)) figs.append(pylab.figure("BGPLVM DEBUG Kmm", figsize=(12, 6))) fig = figs[-1] ax6 = fig.add_subplot(121) @@ -345,16 +345,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): # loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.15, 1, 1.15), # borderaxespad=0, mode="expand") ax2.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], - loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), + loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1), borderaxespad=0, mode="expand") ax3.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], - loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), + loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1), borderaxespad=0, mode="expand") ax4.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], - loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), + loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1), borderaxespad=0, mode="expand") ax5.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], - loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), + loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1), borderaxespad=0, mode="expand") Lleg = ax1.legend() Lleg.draggable() diff --git a/GPy/testing/cgd_tests.py b/GPy/testing/cgd_tests.py index ecd6f829..57a08511 100644 --- a/GPy/testing/cgd_tests.py +++ b/GPy/testing/cgd_tests.py @@ -26,7 +26,7 @@ class Test(unittest.TestCase): for _ in range(restarts): try: x0 = numpy.random.randn(N) * .5 - res = opt.fmin(f, df, x0, messages=0, + res = opt.opt(f, df, x0, messages=0, maxiter=1000, gtol=1e-10) assert numpy.allclose(res[0], 0, atol=1e-3) break @@ -48,7 +48,7 @@ class Test(unittest.TestCase): for _ in range(restarts): try: x0 = numpy.random.randn(N) * .5 - res = opt.fmin(f, df, x0, messages=0, + res = opt.opt(f, df, x0, messages=0, maxiter=1000, gtol=1e-2) assert numpy.allclose(res[0], 1, atol=.01) break @@ -103,7 +103,7 @@ if __name__ == "__main__": if r[-1] != RUNNING: res[0] = r - p, c = opt.fmin_async(f, df, x0.copy(), callback, messages=True, maxiter=1000, + p, c = opt.opt_async(f, df, x0.copy(), callback, messages=True, maxiter=1000, report_every=20, gtol=1e-12) From f4b997beb8e994a40d60ae0a3675801d845e8b00 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 3 May 2013 14:38:42 +0100 Subject: [PATCH 2/3] last opt updates and tests --- GPy/inference/conjugate_gradient_descent.py | 10 ++++++---- GPy/testing/cgd_tests.py | 11 +++++------ 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/GPy/inference/conjugate_gradient_descent.py b/GPy/inference/conjugate_gradient_descent.py index c88249c3..c9981fe8 100644 --- a/GPy/inference/conjugate_gradient_descent.py +++ b/GPy/inference/conjugate_gradient_descent.py @@ -63,14 +63,15 @@ class _Async_Optimization(Thread): return f_w def callback(self, *a): - self.outq.put(a) + if self.outq is not None: + self.outq.put(a) # self.parent and self.parent.callback(*a, **kw) pass # print "callback done" def callback_return(self, *a): self.callback(*a) - self.outq.put(self.SENTINEL) + self.callback(self.SENTINEL) self.runsignal.clear() def run(self, *args, **kwargs): @@ -170,16 +171,17 @@ class Async_Optimize(object): messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6, report_every=10, *args, **kwargs): self.runsignal.set() - outqueue = Queue() c = None + outqueue = None if callback: + outqueue = Queue() self.callback = callback c = Thread(target=self.async_callback_collect, args=(outqueue,)) c.start() p = _CGDAsync(f, df, x0, update_rule, self.runsignal, self.SENTINEL, report_every=report_every, messages=messages, maxiter=maxiter, max_f_eval=max_f_eval, gtol=gtol, outqueue=outqueue, *args, **kwargs) - p.run() + p.start() return p, c def opt(self, f, df, x0, callback=None, update_rule=FletcherReeves, diff --git a/GPy/testing/cgd_tests.py b/GPy/testing/cgd_tests.py index 57a08511..79e5e08b 100644 --- a/GPy/testing/cgd_tests.py +++ b/GPy/testing/cgd_tests.py @@ -14,7 +14,7 @@ from scipy.optimize.optimize import rosen, rosen_der class Test(unittest.TestCase): def testMinimizeSquare(self): - N = 2 + N = 100 A = numpy.random.rand(N) * numpy.eye(N) b = numpy.random.rand(N) * 0 f = lambda x: numpy.dot(x.T.dot(A), x) - numpy.dot(x.T, b) @@ -25,7 +25,7 @@ class Test(unittest.TestCase): restarts = 10 for _ in range(restarts): try: - x0 = numpy.random.randn(N) * .5 + x0 = numpy.random.randn(N) * 300 res = opt.opt(f, df, x0, messages=0, maxiter=1000, gtol=1e-10) assert numpy.allclose(res[0], 0, atol=1e-3) @@ -37,10 +37,9 @@ class Test(unittest.TestCase): raise AssertionError("Test failed for {} restarts".format(restarts)) def testRosen(self): - N = 2 + N = 20 f = rosen df = rosen_der - x0 = numpy.random.randn(N) * .5 opt = CGD() @@ -49,8 +48,8 @@ class Test(unittest.TestCase): try: x0 = numpy.random.randn(N) * .5 res = opt.opt(f, df, x0, messages=0, - maxiter=1000, gtol=1e-2) - assert numpy.allclose(res[0], 1, atol=.01) + maxiter=5e2, gtol=1e-2) + assert numpy.allclose(res[0], 1, atol=.1) break except: # RESTART From 8fb9ab5610ffd1a90132f6843c4dc6df859be5c2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 3 May 2013 14:57:03 +0100 Subject: [PATCH 3/3] BGPLVM example MATLAB compare --- GPy/examples/dimensionality_reduction.py | 4 ---- GPy/models/Bayesian_GPLVM.py | 6 +++--- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index e4fcc234..97601b0f 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -192,10 +192,6 @@ def bgplvm_simulation_matlab_compare(): m['noise'] = Y.var() / 100. m['linear_variance'] = .01 -# lscstr = '{}'.format(k.parts[0].name) -# m[lscstr] = .01 -# m.unconstrain(lscstr); m.constrain_fixed(lscstr, 10) - # lscstr = 'X_variance' # m[lscstr] = .01 # m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1) diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index 099a17ea..fc7e4ba9 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -96,9 +96,9 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): print "\rWARNING: Caught LinAlgError, continueing without setting " if self._debug: self._savederrors.append(self.f_call) -# if save_count > 10: -# raise -# self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1) + if save_count > 10: + raise + self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1) def dKL_dmuS(self): dKL_dS = (1. - (1. / (self.X_variance))) * 0.5