Minimum working examples

basic_gp_example.py

@@ -1,5 +1,6 @@
 """
 This code will run an offline example as-is without any input arguments. Simply run the file in python.
+Use the OptOnFly=True flag to optimise on the fly
 
 In order to use it on a specific machine:
 1. Change the *importlib.import_module('machine_interfaces.machine_interface_example')* to your machine interface name.
@@ -9,6 +10,8 @@
 """
 
 import numpy as np
+import os
+import datetime
 import importlib
 import time
 import matplotlib.pyplot as plt
@@ -19,8 +22,8 @@
 from modules.bayes_optimization import BayesOpt, negUCB, negExpImprove
 from modules.OnlineGP import OGP
 
-saveResultsQ = False
-
+saveResultsQ = True
+OptOnFly = False  # use this to switch between using the hyperparameters from the params file, and optimising these on the fly
 mi_module = importlib.import_module('machine_interfaces.machine_interface_example')
 
 
@@ -36,6 +39,7 @@
 start_point = 0*scan_params['start_point'] #if start_point is set to None, the optimizer will start from the current device settings.
 mi = mi_module.machine_interface(dev_ids = dev_ids, start_point = start_point)  #an isotropic n-dimensional gaussian with amplitude=1, centered at the origin, plus gaussian background noise with std dev = 0.1
 
+
 #create the gp
 ndim = len(dev_ids)
 # GP parameters
@@ -45,44 +49,48 @@
 hyperparams = {'precisionMatrix': gp_precisionmat, 'amplitude_covar': gp_amp, 'noise_variance': gp_noise_variance} 
 gp = OGP(ndim, hyperparams)
 
-#create the bayesian optimizer that will use the gp as the model to optimize the machine 
-opt = BayesOpt(gp, mi, acq_func="UCB", start_dev_vals = mi.x, dev_ids = dev_ids)
-opt.ucb_params = scan_params['ucb_params'] #set the acquisition function parameters
-print('ucb_params',opt.ucb_params)
+#create the bayesian optimizer that will use the gp as the model to optimize the machine
+
+if OptOnFly:
+    # optimize_kernel_on_the_fly is the iteration number to start optimize the kernel's hyperparmaters. If None, no optimization of the hypers during BO.
+    optimize_kernel_on_the_fly = 2
+    opt = BayesOpt(gp, mi, acq_func="UCB", start_dev_vals=mi.x, dev_ids=dev_ids,
+                   optimize_kernel_on_the_fly=optimize_kernel_on_the_fly)
+    # opt.ucb_params = scan_params['ucb_params']  # set the acquisition function parameters
+else:
+    optimize_kernel_on_the_fly = None
+    opt = BayesOpt(gp, mi, acq_func="UCB", start_dev_vals = mi.x, dev_ids = dev_ids)
+    opt.ucb_params = scan_params['ucb_params'] #set the acquisition function parameters
 
 #run the gp search for some number of steps
 Obj_state_s=[]
 
-optimize_kernel_on_the_fly = None #optimize_kernel_on_the_fly is the iteration number to start optimize the kernel's hyperparmaters. If None, no optimization of the hypers during BO. 
-
-Niter = 10
+Niter = 100
 for i in range(Niter):
     clear_output(wait=True) 
     print ('iteration =', i)
     print ('current position:', mi.x, 'current objective value:', mi.getState()[1])
 
     Obj_state_s.append(mi.getState()[1][0])
-
-    f = plt.figure(figsize=(20,3))
-    ax = f.add_subplot(121)
-    ax2 = f.add_subplot(122)
-    ax.set_ylabel('Input controls',fontsize=12)
-    ax.set_xlabel('Iteration',fontsize=12)
-    for x, label in zip(opt.X_obs.T, opt.dev_ids):
-        ax.plot(x,'.-',label = label)
-    ax.legend()
-    ax2.set_ylabel('Objective',fontsize=12)
-    ax2.set_xlabel('Iteration',fontsize=12)
-    ax2.plot(Obj_state_s,'.-')
-    plt.show(); 
-
     if optimize_kernel_on_the_fly is not None:
         if i > optimize_kernel_on_the_fly:
-            opt.optimize_kernel_hyperparameters()    
-
+            opt.optimize_kernel_hyperparameters()
     opt.OptIter()
-    time.sleep(acquisition_delay)   
-
+    time.sleep(acquisition_delay)
+
+f = plt.figure(figsize=(20, 3))
+ax = f.add_subplot(121)
+ax2 = f.add_subplot(122)
+ax.set_ylabel('Input controls', fontsize=12)
+ax.set_xlabel('Iteration', fontsize=12)
+for x, label in zip(opt.X_obs.T, opt.dev_ids):
+    ax.plot(x, '.-', label=label)
+ax.legend()
+ax2.set_ylabel('Objective', fontsize=12)
+ax2.set_xlabel('Iteration', fontsize=12)
+ax2.plot(Obj_state_s, '.-')
+plt.show()
+
 #save results if desired
 if saveResultsQ == True:
     timestr = datetime.now().strftime('%Y-%m-%d-%H%M%S')

modules/bayes_optimization.py

@@ -583,6 +583,6 @@ def negUCB(x_new, model, ndim, nsteps, nu = 1., delta = 1.):
         tau = 2.*np.log(nsteps**(0.5*ndim+2.)*(np.pi**2.)/3./delta)
         GPUCB = y_mean + np.sqrt(nu * tau * y_var)
 
-    return -GPUCB
+    return -GPUCB[0]
 
 

modules/parallelstuff.py

@@ -7,6 +7,7 @@
 import numpy as np
 import multiprocessing as mp
 import copy
+import time
 
 # handle 'IOError: [Errno 4] Interrupted system call' errors from multiprocessing.Queue.get
 #https://stackoverflow.com/questions/14136195/what-is-the-proper-way-to-handle-in-python-ioerror-errno-4-interrupted-syst
@@ -34,7 +35,7 @@ def mworker(f,x0,fargs,margs,out_q):
         #print 'worker: margs = ',margs
         res = minimize(f, x0, args = fargs, **margs)
         #return [res.x, res.fun]
-        out_q.put([[res.x, res.fun[0][0]]])
+        out_q.put([[res.x, res.fun[0]]])
 
     # parallelize minimizations using different starting positions using multiprocessing, scipy.optimize.minimize
     def parallelminimize(f,x0s,fargs,margs,v0best=None,relative_bounds=None):
@@ -337,7 +338,7 @@ def parallelmap2(f,fargslist,hostlist=None):
     def eworker(f,x,fargs,out_q):
         # worker invoked in a process puts the results in the output queue out_q
         res = f(x, *fargs)
-        out_q.put(np.hstack((x, res[0][0])))
+        out_q.put(np.hstack((x, res[0])))
 
     # eval function over a range of initial points neval and return the nkeep lowest function evals
     def parallelgridsearch(f,x0,lengths,fargs,neval,nkeep):