<jupyter_start><jupyter_text>Preferential Bayesian optimisation with Skew Gaussian Processes PBO-SkewGP benchmarkThis notebook allows to reproduce the experiments in _"Preferential Bayesian optimisation with SkewGaussian Processes"_ , Section 6.1. **Note:** The GP-Laplace implementation is in Matlab. In order to run it you will need a working Matlab licence and Python needs to be able to call it with the `matlab.engine` API.<jupyter_code># Options
nameBenchmark= 'goldstein' # 'Forrester' 'sixhump' 'levy' 'rosenbrock' 'hartman6'
acquisitionselected=  'Thompson' #'EPI_IGAIN' # 'UCB' #'Thompson'
num_repetitions = 1 # in the paper 20
iterations = 20 # in the paper 100
runGPL = False # NOTE THAT GP-Laplace is run on a MATLAB implementation.
# Required Python packages
%load_ext autoreload
%autoreload 2
from scipy.optimize import minimize
import os.path
import numpy as np
import pandas as pd
import GPy as GPy
from matplotlib import pyplot as plt
from scipy.stats import multivariate_normal, bernoulli
from scipy.spatial.distance import cdist
from sklearn.preprocessing import StandardScaler
from noisyopt import minimizeSPSA
import pymc3 as pm
import pickle
import time
import seaborn as snb



# Local packages
import sys
sys.path.append("..")
import SkewGP as SkewGP
import BO as BO_SGP 
import commoncode as commoncode

%matplotlib inline
from get_oracle import get_oracle
from get_oracle import valid
oracle, bounds, optimum_points = get_oracle(nameBenchmark)
# Functions to save and load results
def save_obj(obj, name):
    with open('obj/'+ name + '.pkl','wb') as file:
        pickle.dump(obj, file)
        
def load_obj(name):
    with open('obj/'+ name + '.pkl','rb') as file:
        return pickle.load(file)
# Create the directory obj
path='./obj'


if not os.path.isdir(path):   
    try:
        os.mkdir(path)
    except OSError:
        print ("Creation of the directory %s failed" % path)
    else:
        print ("Successfully created the directory %s " % path)
else:
    print ("The directory %s already exits." % path)
def queryf(x_new,x_old,i,j,f,valid):
    if valid(x_new,f)==1.0:
        v1=f(x_new)
        v2=f(x_old)
        if v1>=v2:
            return  1.0, [i,j]
        elif v1<v2:
            return  1.0, [j,i]
        else:
            print("error")
    else:
        return  -1.0,[]


class RandomSearch():
    def __init__(self,X,Preference, Class,bounds,oracle,valid):
        self.bounds=bounds
        self.iref = 0
        self.Xref = X[0:1,:]
        self.X=X.copy()
        self.Preference=Preference.copy()
        self.Class=Class.copy()
        self.oracle=oracle
        self.valid=valid  


        
    def find_next(self,oracle):
                
        bb = np.vstack(self.bounds)
        Xnew = (bb[:,0:1]+(bb[:,1:2]-bb[:,0:1])*np.random.rand(self.X.shape[1],1)).T     
        self.X=np.vstack([self.X,Xnew])
        inew = self.X.shape[0]-1
        cl,pref = queryf(Xnew,self.Xref,inew,self.iref,self.oracle,self.valid)
        print("RandomSearch:", Xnew,self.Xref,self.oracle(Xnew),self.oracle(self.Xref),cl,pref)
        print(" ")
        if cl==1.0:
            if pref[0]==inew:
                self.iref = inew #we have found a better point
                self.Xref = Xnew #X[-1:,:]
        self.Class.append(cl)    
        if len(pref)>0:
            self.Preference.append(pref)
        print("RandomSearch Xref=",self.Xref)
        print(" ")<jupyter_output><empty_output><jupyter_text>BO Here the kernel used is RBF defined by the function below<jupyter_code>#define RBF kernel function
def Kernel(X1,X2,params,diag_=False):
    lengthscale=params['lengthscale']['value']
    variance   =params['variance']['value']
    if diag_==False:
        diffs = cdist(np.atleast_2d(X1)/ lengthscale, np.atleast_2d(X2) / lengthscale, metric='sqeuclidean')
    else:
        diffs = np.sum((np.atleast_2d(X1)/ lengthscale-np.atleast_2d(X2)/ lengthscale)*(np.atleast_2d(X1)/ lengthscale-np.atleast_2d(X2)/ lengthscale),axis=1)
    return variance * np.exp(-0.5 * diffs)<jupyter_output><empty_output><jupyter_text>The cell below does the BO experiment<jupyter_code>mBOSGP_all=[]
mBOGPL_all=[]
RS_all=[]
ALL_RT=[]
for trial in range(0,num_repetitions):
    
    # Check if trial already exists 
    if os.path.isfile('obj/'+ nameBenchmark+'_'+acquisitionselected+'_mBOSGP_trial'+str(trial)+'.pkl'):
        mBOSGP_Xref = load_obj(nameBenchmark+'_'+acquisitionselected+'_mBOSGP_trial'+str(trial))
        mBOGPL_Xref = load_obj(nameBenchmark+'_'+acquisitionselected+'_mBOGPL_trial'+str(trial))
        RS_Xref = load_obj(nameBenchmark+'_'+acquisitionselected+'_RS_trial'+str(trial))
        temp_RT = load_obj(nameBenchmark+'_'+acquisitionselected+'_running_time_SGP_GPL_trial'+str(trial))
        
        mBOSGP_all.append(mBOSGP_Xref)
        mBOGPL_all.append(mBOGPL_Xref)
        RS_all.append(RS_Xref)
        ALL_RT.append(temp_RT)
        save_obj(mBOSGP_all,nameBenchmark+'_'+acquisitionselected+'_mBOSGP_all')
        save_obj(mBOGPL_all,nameBenchmark+'_'+acquisitionselected+'_mBOGPL_all')
        save_obj(RS_all,nameBenchmark+'_'+acquisitionselected+'_RS_all')
        save_obj(ALL_RT,nameBenchmark+'_'+acquisitionselected+'_running_time_SGP_GPL')
        
        print('Trial ',trial,' already done, skipping it\n\n\n\n')
        continue
        
    #Initial design
    np.random.seed(trial*101+42)
    bb = np.vstack(bounds)
    X= bb[:,0]+(bb[:,1]-bb[:,0])*np.random.rand(1,len(bounds))
    Class=[1.0]
    Preference= []
    iref = 0
    Xref = X[0:1,:]
    initial_points = 10

    for ii in range(initial_points):    
        Xnew = bb[:,0]+(bb[:,1]-bb[:,0])*np.random.rand(1,len(bounds))
        X=np.vstack([X,Xnew])
        inew = X.shape[0]-1
        cl,pref = BO_SGP.queryf(Xnew,Xref,inew,iref,oracle,valid)
        Class.append(cl)    
        if len(pref)>0:
            Xref=Xnew
            iref=inew
            Preference.append(pref)

    maX=X[np.argmax(oracle(X))]       

    # Initialize kernel
    kernel = Kernel
    
    #Hyperparameters of the kernel
    logexp=commoncode.logexp()
    params={'lengthscale': {'value':np.ones((1,X.shape[1]))*0.3, 
                    'range':np.vstack([[np.exp(-5.0), np.exp(5.0)]]*X.shape[1]),
                    'transform': logexp},
         'variance': {'value':np.array([30.0]), 
                    'range':np.vstack([[np.exp(-5.0), np.exp(4.1)]]),
                    'transform': logexp},
            'noise_variance': {'value':np.array([1.0]),
                               'range':np.vstack([[1.0, 1.0001]]),
                               'transform': logexp}
      }
    
    
    
    Xinit=X.copy()
    
    # Initialize models 
    mBOSGP = BO_SGP.BO(X,Preference, Class, bounds, kernel,params,oracle,valid,maX,
                       alternate_optim=30,nsamples=2000,surrogateM='SGP',acquisition=acquisitionselected)
    mBOGPL = BO_SGP.BO(X,Preference, Class, bounds, kernel,params,oracle,valid,maX,
                       alternate_optim=30,nsamples=2000,surrogateM='GPL',acquisition=acquisitionselected)

    RS = RandomSearch(X,Preference, Class, bounds,oracle,valid)
    
    print("============Initial input points ============")
    print(X)

    ##
    update_model=True
    max_iter=iterations
    ϵx_toll=0.02
    mBOSGP_Xref=[]
    mBOGPL_Xref=[]
    RS_Xref=[]
    SGP_time=0
    GPL_time=0
    for i in range(1,max_iter):
        np.random.seed(i*(trial+1))
        
        print("============Iteration number ",i," ============")
        
        
        print("============1: PBO-SkewGP ============")
        minv=np.inf
        for opt in optimum_points:
            dd=np.max(np.abs(mBOSGP.Xref-np.array(opt)))
            if minv>dd:
                minv=dd
        
        if  dd>ϵx_toll:
            start_time = time.time()
            mBOSGP.find_next(oracle,i,update_model=update_model)    
            mBOSGP_Xref.append(mBOSGP.Xref)
            SGP_time=SGP_time+time.time() - start_time
        else:
            mBOSGP_Xref.append(mBOSGP.Xref)
            
        print("============2: RandSampl ============")
        RS.find_next(oracle)
        RS_Xref.append(RS.Xref)
        
        if runGPL:
            print("============3: PBO-GPL ============")
            minv=np.inf
            for opt in optimum_points:
                dd=np.max(np.abs(mBOGPL.Xref-np.array(opt)))
                if minv>dd:
                    minv=dd
            if  minv>ϵx_toll:
                start_time = time.time()
                mBOGPL.find_next(oracle,i,update_model=update_model)    
                mBOGPL_Xref.append(mBOGPL.Xref)
                GPL_time=GPL_time+time.time() - start_time
            else:
                mBOGPL_Xref.append(mBOGPL.Xref)
        
    print(mBOSGP.Xref)
    mBOSGP_all.append(mBOSGP_Xref)
    RS_all.append(RS_Xref)
    
    if runGPL:
        mBOGPL_all.append(mBOGPL_Xref)
        temp_RT = [SGP_time,GPL_time]
    else:
        temp_RT = [SGP_time]
    
    ALL_RT.append(temp_RT)
    save_obj(mBOSGP_all,nameBenchmark+'_'+acquisitionselected+'_mBOSGP_all')
    save_obj(RS_all,nameBenchmark+'_'+acquisitionselected+'_RS_all')
    save_obj(ALL_RT,nameBenchmark+'_'+acquisitionselected+'_running_time_SGP_GPL')
    
    
    save_obj(mBOSGP_Xref,nameBenchmark+'_'+acquisitionselected+'_mBOSGP_trial'+str(trial))
    save_obj(RS_Xref,nameBenchmark+'_'+acquisitionselected+'_RS_trial'+str(trial))
    save_obj(temp_RT,nameBenchmark+'_'+acquisitionselected+'_running_time_SGP_GPL_trial'+str(trial))
    
    if runGPL:
        save_obj(mBOGPL_Xref,nameBenchmark+'_'+acquisitionselected+'_mBOGPL_trial'+str(trial))
        save_obj(mBOGPL_all,nameBenchmark+'_'+acquisitionselected+'_mBOGPL_all')

    
    
    del mBOSGP
    del mBOGPL<jupyter_output><empty_output><jupyter_text>Re-load the results and plot<jupyter_code>mBOSGP_all=load_obj(nameBenchmark+'_'+acquisitionselected+'_mBOSGP_all')
RS_all=load_obj(nameBenchmark+'_'+acquisitionselected+'_RS_all')

if runGPL:
    mBOGPL_all=load_obj(nameBenchmark+'_'+acquisitionselected+'_mBOGPL_all')




i=0
valmin=1.0 #np.min(np.vstack([-oracle(np.vstack(RS_all[i])),-oracle(np.vstack(mBOSGP_all[i])),-oracle(np.vstack(mBOGPL_all[i]))]))
SGP=[-oracle(np.vstack(mBOSGP_all[i])).T/np.abs(valmin)]
RND=[-oracle(np.vstack(RS_all[i])).T/np.abs(valmin)]

if runGPL:
    GPL=[-oracle(np.vstack(mBOGPL_all[i])).T/np.abs(valmin)]
    

for i in range(0,len(mBOSGP_all)):
    valmin=1.0 #np.min(np.vstack([-oracle(np.vstack(RS_all[i])),-oracle(np.vstack(mBOSGP_all[i])),-oracle(np.vstack(mBOGPL_all[i]))]))
    SGP.append(-oracle(np.vstack(mBOSGP_all[i])).T/np.abs(valmin))
    RND.append(-oracle(np.vstack(RS_all[i])).T/np.abs(valmin))
    if runGPL:
        GPL.append(-oracle(np.vstack(mBOGPL_all[i])).T/np.abs(valmin))

SGP=np.vstack(SGP).T
RND=np.vstack(RND).T

if runGPL:
    GPL=np.vstack(GPL).T


plt.rc('font', family='serif')
plt.rc('xtick', labelsize=10)
plt.rc('ytick', labelsize=10)
plt.figure(figsize=(7,3))
plt.title(nameBenchmark,fontsize=14)
acqu=acquisitionselected

xx = np.arange(0,len(mBOSGP_all[i]))
plt.plot(xx,np.mean(SGP,axis=1),label='SkewGP '+acqu,color='C0')
plt.plot(xx,np.mean(RND,axis=1),label='random', color='C1',linestyle='dashed')
if runGPL:
    plt.plot(xx,np.mean(GPL,axis=1),label='GPL '+acqu,color='C2',linestyle='dashdot')
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14934436328434392, 'fraction_numerical': 0.013552758954501452, 'mean_word_length': 3.168745414526779, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '87289', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4405, 'n_tokens_neox': 4173, 'n_words': 626, 'path': 'notebooks/Benchmarks PBO-SkewGP.ipynb', 'repo_name': 'benavoli/SkewGP', 'size': '11363'}