gpOptimizer: Single-Task Acquisition Functions#
#!pip install gpcam==8.4.0
#!pip install matplotlib
Setup#
import numpy as np
import matplotlib.pyplot as plt
from gpcam import GPOptimizer
import time
from loguru import logger
from distributed import Client
client = Client()
%load_ext autoreload
%autoreload 2
from itertools import product
x_pred1D = np.linspace(0,1,1000).reshape(-1,1)
Data Preparation#
x = np.linspace(0,600,1000)
def f1(x):
return np.sin(5. * x) + np.cos(10. * x) + (2.* (x-0.4)**2) * np.cos(100. * x)
x_data = np.random.rand(50).reshape(-1,1)
y_data = f1(x_data[:,0]) + (np.random.rand(len(x_data))-0.5) * 0.5
plt.figure(figsize = (15,5))
plt.xticks([0.,0.5,1.0])
plt.yticks([-2,-1,0.,1])
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.plot(x_pred1D,f1(x_pred1D), color = 'orange', linewidth = 4)
plt.scatter(x_data[:,0],y_data, color = 'black')
<matplotlib.collections.PathCollection at 0x7f3914597290>
Customizing the Gaussian Process#
def my_noise(x,hps):
#This is a simple noise function but can be made arbitrarily complex using many hyperparameters.
#The noise function can return a matrix or a vector
return np.zeros((len(x))) + hps[2]
#stationary
from gpcam.kernels import *
def skernel(x1,x2,hps):
#The kernel follows the mathematical definition of a kernel. This
#means there is no limit to the variety of kernels you can define.
d = get_distance_matrix(x1,x2)
return hps[0] * matern_kernel_diff1(d,hps[1])
def meanf(x, hps):
#This is a simple mean function but it can be arbitrarily complex using many hyperparameters.
return 1.-np.sin(hps[3] * x[:,0])
#it is a good idea to plot the prior mean function to make sure we did not mess up
plt.figure(figsize = (15,5))
plt.plot(x_pred1D,meanf(x_pred1D, np.array([1.,1.,5.0,2.])), color = 'orange', label = 'task1')
[<matplotlib.lines.Line2D at 0x7f39146267d0>]
Initialization and Different Training Options#
my_gpo = GPOptimizer(x_data,y_data,
init_hyperparameters = np.ones((4))/10., # We need enough of those for kernel, noise, and prior mean functions
compute_device='cpu',
kernel_function=skernel,
kernel_function_grad=None,
prior_mean_function=meanf,
prior_mean_function_grad=None,
noise_function=my_noise,
#noise_variances=np.zeros(y_data.shape) + 0.1,
gp2Scale = False,
ram_economy=False,
args={'a': 1.5, 'b':2.},
)
hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
[0.01,10.], #length scale for the kernel
[0.00001,0.1], #noise
[0.00001,1.] #mean
])
#the following is not needed, this is just to show how data is replced or appended
x_update = np.array([0.1,0.2,0.5]).reshape(3,1)
y_update = f1(x_update[:,0]) + (np.random.rand(len(x_update))-0.5) * 0.5
my_gpo.tell(x_update, y_update, append=True) ##append to the data
my_gpo.tell(x_data, y_data, append=False) ## back to normal overwriting the updated data
st = time.time()
print("Standard Training (MCMC)")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
print("ADAM")
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100, method="adam")
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
print("Global Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='global', max_iter = 20)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
print("Local Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='local')
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
print("HGDL Training")
my_gpo.train(hyperparameter_bounds=hps_bounds, method='hgdl', max_iter=2, dask_client=client)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
Standard Training (MCMC)
Starting likelihood. f(x)= -56.32785658372028
Finished 10 out of 100 iterations. f(x)= -36.08840067752926
Finished 20 out of 100 iterations. f(x)= -36.08840067752926
Finished 30 out of 100 iterations. f(x)= -21.2115579270182
Finished 40 out of 100 iterations. f(x)= -13.985568700479519
Finished 50 out of 100 iterations. f(x)= -13.559537891389219
Finished 60 out of 100 iterations. f(x)= -12.86877237412611
Finished 70 out of 100 iterations. f(x)= -14.569584563969997
Finished 80 out of 100 iterations. f(x)= -14.475330777730239
Finished 90 out of 100 iterations. f(x)= -13.737521046849132
Result= [1.61575711 0.26539302 0.08954816 0.05094794] after 0.03347516059875488 seconds
ADAM
Result= [2.29705964 0.46311201 0.0271336 0.87073776] after 0.12983918190002441 seconds
Global Training
Result= [2.29705964 0.46311201 0.0271336 0.87073776] after 0.4582846164703369 seconds
Local Training
Result= [2.29705964 0.46311201 0.0271336 0.87073776] after 0.4632253646850586 seconds
HGDL Training
Result= [2.29705964 0.46311201 0.0271336 0.87073776] after 1.3452765941619873 seconds
Asynchronous Training#
Train asynchronously – via Adam, HGDL, or MCMC – on a remote server or locally. You can also start a bunch of different training runs on different computers. This training will continue without any signs of life until you query the solution via ‘update_hyperparameters(object)’ or call ‘my_gpo.stop_training(opt_obj)’
HGDL#
my_gpo.set_hyperparameters(np.ones((4))/10.)
opt_obj = my_gpo.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method="hgdl")
print(my_gpo.hyperparameters)
for i in range(20):
my_gpo.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gpo.hyperparameters)
time.sleep(0.1)
my_gpo.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[0.1 0.1 0.1 0.1]
iteration 0 : [0.1 0.1 0.1 0.1]
iteration 1 : [4.83429379 0.62844593 0.02707247 1. ]
/home/marcus/Coding/fvGP/fvgp/gp.py:881: UserWarning: Hyperparameter update not successful len(optima list) = 0
hps = self.trainer.update_hyperparameters(opt_obj)
iteration 2 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 3 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 4 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 5 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 6 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 7 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 8 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 9 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 10 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 11 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 12 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 13 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 14 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 15 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 16 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 17 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 18 : [4.83429379 0.62844593 0.02707247 1. ]
iteration 19 : [4.83429379 0.62844593 0.02707247 1. ]
ADAM#
my_gpo.set_hyperparameters(np.array([1.,1.,1.,1.]))
print(my_gpo.hyperparameters)
opt_obj = my_gpo.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method='adam')
# The result won't change much (or at all) since this is such a simple optimization
for i in range(20):
my_gpo.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gpo.hyperparameters)
time.sleep(0.1)
my_gpo.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1. 1. 1.]
iteration 0 : [0.48672518 5.01203075 0.01398589 0.9175439 ]
iteration 1 : [0.80432051 4.65849955 0.18054652 0.97946032]
iteration 2 : [0.99636743 4.33849772 0.23052178 0.96805725]
iteration 3 : [1.16995333 3.95739798 0.26861871 0.94899815]
iteration 4 : [1.34229071 3.46733371 0.29925316 0.92163628]
iteration 5 : [1.53143749 2.75221342 0.32286231 0.88736197]
iteration 6 : [1.80191368 1.23268828 0.33396493 0.9128212 ]
iteration 7 : [1.91748234 0.35630937 0.30976529 1.20047573]
iteration 8 : [1.92584739 0.36431463 0.2772591 1.32277396]
iteration 9 : [1.93649147 0.37028064 0.23624153 1.36300871]
iteration 10 : [1.95080556 0.37985907 0.18372561 1.36945446]
iteration 11 : [1.97223955 0.39965644 0.10660943 1.36058614]
iteration 12 : [2.00800891 0.43718714 0.02824416 1.32807693]
iteration 13 : [2.05539256 0.4423543 0.02712094 1.31179315]
iteration 14 : [2.1026721 0.44665996 0.02711375 1.30805829]
iteration 15 : [2.15057573 0.45094608 0.02711201 1.30540593]
iteration 16 : [2.20190304 0.45546744 0.02711022 1.30270665]
iteration 17 : [2.25275236 0.45987711 0.02710855 1.30009947]
iteration 18 : [2.30115361 0.46401275 0.02710704 1.29767356]
iteration 19 : [2.35233505 0.46832258 0.02710553 1.29516458]
MCMC#
my_gpo.set_hyperparameters(np.array([1.,1.,1.,1.]))
print(my_gpo.hyperparameters)
opt_obj = my_gpo.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method='mcmc')
# The result won't change much (or at all) since this is such a simple optimization
for i in range(20):
my_gpo.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gpo.hyperparameters)
time.sleep(0.1)
my_gpo.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1. 1. 1.]
iteration 0 : [1. 1. 1. 1.]
iteration 1 : [5.01858325 0.71143118 0.04832434 0.7940703 ]
iteration 2 : [8.57830742 1.21681742 0.03310934 0.31085877]
iteration 3 : [1.62706714 0.50502652 0.02752291 0.40039741]
iteration 4 : [5.30963736 0.7304087 0.02918963 0.45076057]
iteration 5 : [7.1664687 0.50012083 0.02342949 0.23741583]
iteration 6 : [7.16628915 0.74647664 0.02704181 0.56496979]
iteration 7 : [9.11034635 0.57559477 0.02041948 0.72403078]
iteration 8 : [8.18844805 0.88808423 0.02141499 0.86949354]
iteration 9 : [7.47179674 0.92895819 0.03001137 0.43229148]
iteration 10 : [9.03150409 0.60726481 0.03688093 0.63894609]
iteration 11 : [8.40712502 0.75367922 0.02758517 0.93494694]
iteration 12 : [8.96426077 0.70782725 0.02599546 0.86242299]
iteration 13 : [8.17511137 0.95154084 0.02622233 0.9785727 ]
iteration 14 : [9.70345414 0.94918753 0.02793636 0.41505202]
iteration 15 : [9.29576049 0.99797824 0.03201542 0.48639497]
iteration 16 : [9.35690593 0.89997536 0.03366795 0.44038127]
iteration 17 : [9.35690593 0.89997536 0.03366795 0.44038127]
iteration 18 : [9.35690593 0.89997536 0.03366795 0.44038127]
iteration 19 : [9.35690593 0.89997536 0.03366795 0.44038127]
Vizualizing the Results#
#let's make a prediction
x_pred = np.linspace(0,1,1000)
hps = my_gpo.train(hyperparameter_bounds=hps_bounds, info = False)
# different ways to call
var1 = my_gpo.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=False)["v(x)"]
var1 = my_gpo.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=True)["v(x)"]
mean1 = my_gpo.posterior_mean(x_pred.reshape(-1,1))["m(x)"]
var1 = my_gpo.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=True)["v(x)"]
mean_grad = my_gpo.posterior_mean_grad(x_pred.reshape(-1,1), direction=0)["dm/dx"]
print("Posterior Mean and Uncertainty")
plt.figure(figsize = (16,10))
plt.plot(x_pred,mean1, label = "posterior mean", linewidth = 4)
plt.plot(x_pred1D,f1(x_pred1D), label = "latent function", linewidth = 4)
plt.fill_between(x_pred, mean1 - 3. * np.sqrt(var1), mean1 + 3. * np.sqrt(var1), alpha = 0.5, color = "grey", label = "var")
plt.scatter(my_gpo.x_data,my_gpo.y_data, color = 'black')
plt.show()
print("Posterior Mean Gradient")
plt.figure(figsize = (16,10))
dx = 1./len(x_pred)
plt.plot(x_pred1D,np.gradient(f1(x_pred1D).flatten(), dx), label = "ground truth gradient", linewidth = 4)
plt.plot(x_pred1D,mean_grad, label = "posterior mean grad", linewidth = 4)
plt.show()
##looking at some validation metrics
print("RMSE: ",my_gpo.rmse(x_pred1D,f1(x_pred1D).flatten()))
print("NRMSE: ",my_gpo.nrmse(x_pred1D,f1(x_pred1D).flatten()))
print("CRPS (mean, std): ",my_gpo.crps(x_pred1D,f1(x_pred1D).flatten()))
print("R2: ",my_gpo.r2(x_pred1D,f1(x_pred1D).flatten()))
print("NLPD: ",my_gpo.nlpd(x_pred1D,f1(x_pred1D).flatten()))
print("MSLL: ",my_gpo.msll(x_pred1D,f1(x_pred1D).flatten()))
print("MAPE: ",my_gpo.mape(x_pred1D,f1(x_pred1D).flatten()))
print("INTERVAL SCORE: ",my_gpo.interval_score(x_pred1D,f1(x_pred1D).flatten()))
print("MPIW: ",my_gpo.mpiw(x_pred1D))
print("PICP: ",my_gpo.picp(x_pred1D,f1(x_pred1D).flatten()))
print("Coverage Curve:")
cov_curve = my_gpo.coverage_curve(x_pred1D,f1(x_pred1D).flatten())
plt.scatter(cov_curve["target_coverage"], cov_curve["measured_coverage"])
plt.show()
print("predicted vs. observed")
my_gpo.plot_observed_vs_predicted(x_pred1D,f1(x_pred1D).flatten())
Posterior Mean and Uncertainty
Posterior Mean Gradient
RMSE: 0.26551684915712864
NRMSE: 0.0673233503333902
CRPS (mean, std): (np.float64(0.129899673562042), np.float64(0.19291581033212993))
R2: 0.9337822933150128
NLPD: 1.3646130832100094
MSLL: -0.09332464915864147
MAPE: 1.0978656035247554
INTERVAL SCORE: 1.7400799163782434
MPIW: 0.775171272254288
PICP: 0.914
Coverage Curve:
predicted vs. observed
#available acquisition function for the single-task case:
acquisition_functions = ["variance","relative information entropy","relative information entropy set",
"ucb","lcb","maximum","minimum","gradient","expected improvement",
"probability of improvement", "target probability", "total correlation"]
plt.figure(figsize=(16,10))
for acq_func in acquisition_functions:
print("Acquisition function ",acq_func)
res = my_gpo.evaluate_acquisition_function(x_pred, acquisition_function=acq_func)
if len(res)==len(x_pred):
res = res - np.min(res)
res = res/np.max(res)
plt.plot(x_pred,res, label = acq_func, linewidth = 2)
else: print("Some acquisition function return a scalar score for the entirety of points. Here: ", acq_func)
plt.legend()
plt.show()
Acquisition function variance
Acquisition function relative information entropy
Some acquisition function return a scalar score for the entirety of points. Here: relative information entropy
Acquisition function relative information entropy set
Acquisition function ucb
Acquisition function lcb
Acquisition function maximum
Acquisition function minimum
Acquisition function gradient
Acquisition function expected improvement
Acquisition function probability of improvement
Acquisition function target probability
Acquisition function total correlation
Some acquisition function return a scalar score for the entirety of points. Here: total correlation
ask()ing for Optimal Evaluations#
with several optimization methods and acquisition functions
#let's test the asks:
bounds = np.array([[0.0,1.0]])
for acq_func in acquisition_functions:
for method in ["global","local","hgdl"]:
print("Acquisition function ", acq_func," and method ",method)
new_suggestion = my_gpo.ask(bounds, acquisition_function=acq_func,
method=method, max_iter = 2, dask_client=client)
print("led to new suggestion: \n", new_suggestion)
print("")
Acquisition function variance and method global
led to new suggestion:
{'x': array([[0.98235171]]), 'f_a(x)': array([0.18126027]), 'opt_obj': None}
Acquisition function variance and method local
led to new suggestion:
{'x': array([[0.63044797]]), 'f_a(x)': array([0.07434635]), 'opt_obj': None}
Acquisition function variance and method hgdl
[[0.77672333]
[0.35066261]
[0.35065813]] [0 1]
led to new suggestion:
{'x': array([[0.77672333]]), 'f_a(x)': array([0.09543538]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f391464d8d0>}
Acquisition function relative information entropy and method global
led to new suggestion:
{'x': array([[0.00070177]]), 'f_a(x)': array([-87.54477625]), 'opt_obj': None}
Acquisition function relative information entropy and method local
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([-153.04357747]), 'opt_obj': None}
Acquisition function relative information entropy and method hgdl
/home/marcus/Coding/gpCAM/gpcam/gp_optimizer_base.py:430: UserWarning: I set vectorized=False for total corr. or rel. inf. entropy.
warnings.warn("I set vectorized=False for total corr. or rel. inf. entropy.")
[[0.]
[0.]
[1.]] [0 2]
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([-85.5491618]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f3914459790>}
Acquisition function relative information entropy set and method global
led to new suggestion:
{'x': array([[0.99943054]]), 'f_a(x)': array([-155.44095197]), 'opt_obj': None}
Acquisition function relative information entropy set and method local
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([-85.5491618]), 'opt_obj': None}
Acquisition function relative information entropy set and method hgdl
[[0.]
[0.]
[1.]] [0 2]
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([-85.5491618]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38f4fbc8d0>}
Acquisition function ucb and method global
led to new suggestion:
{'x': array([[0.00841685]]), 'f_a(x)': array([1.58684509]), 'opt_obj': None}
Acquisition function ucb and method local
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([1.66523528]), 'opt_obj': None}
Acquisition function ucb and method hgdl
[[0.]
[0.]
[0.]] [0]
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([1.66523528]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f392efd2f90>}
Acquisition function lcb and method global
led to new suggestion:
{'x': array([[0.99620377]]), 'f_a(x)': array([3.47121735]), 'opt_obj': None}
Acquisition function lcb and method local
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([3.52975068]), 'opt_obj': None}
Acquisition function lcb and method hgdl
[[1. ]
[0.32229916]
[0.32229963]] [0 1]
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([3.52975068]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38ecf64d50>}
Acquisition function maximum and method global
led to new suggestion:
{'x': array([[0.58132222]]), 'f_a(x)': array([1.13657902]), 'opt_obj': None}
Acquisition function maximum and method local
led to new suggestion:
{'x': array([[0.34663431]]), 'f_a(x)': array([0.10873133]), 'opt_obj': None}
Acquisition function maximum and method hgdl
[[0.58158734]
[0.5815866 ]
[0. ]
[0. ]] [0 2]
led to new suggestion:
{'x': array([[0.58158734]]), 'f_a(x)': array([1.13658189]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38f40b9390>}
Acquisition function minimum and method global
led to new suggestion:
{'x': array([[0.97244745]]), 'f_a(x)': array([2.65457415]), 'opt_obj': None}
Acquisition function minimum and method local
led to new suggestion:
{'x': array([[0.52972925]]), 'f_a(x)': array([-1.04841067]), 'opt_obj': None}
Acquisition function minimum and method hgdl
[[1. ]
[1. ]
[0.30698471]] [0 2]
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([2.82101469]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f3914563610>}
Acquisition function gradient and method global
led to new suggestion:
{'x': array([[0.9996146]]), 'f_a(x)': array([1.24474969]), 'opt_obj': None}
Acquisition function gradient and method local
led to new suggestion:
{'x': array([[0.85467561]]), 'f_a(x)': array([0.91225811]), 'opt_obj': None}
Acquisition function gradient and method hgdl
[[1. ]
[1. ]
[0.77699718]] [0 2]
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([1.24685241]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38ecd92110>}
Acquisition function expected improvement and method global
led to new suggestion:
{'x': array([[0.9794936]]), 'f_a(x)': array([0.06921434]), 'opt_obj': None}
Acquisition function expected improvement and method local
led to new suggestion:
{'x': array([[0.73320862]]), 'f_a(x)': array([0.03282211]), 'opt_obj': None}
Acquisition function expected improvement and method hgdl
[[0. ]
[0. ]
[0.16479685]] [0 2]
led to new suggestion:
{'x': array([[0.]]), 'f_a(x)': array([0.07134305]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f398851ee10>}
Acquisition function probability of improvement and method global
led to new suggestion:
{'x': array([[0.59528619]]), 'f_a(x)': array([0.22813487]), 'opt_obj': None}
Acquisition function probability of improvement and method local
led to new suggestion:
{'x': array([[0.8751773]]), 'f_a(x)': array([0.]), 'opt_obj': None}
Acquisition function probability of improvement and method hgdl
[[0.68081168]
[0.96320561]
[0.91176481]] [0 1 2]
led to new suggestion:
{'x': array([[0.68081168]]), 'f_a(x)': array([1.19764239e-12]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38ecb56e10>}
Acquisition function target probability and method global
led to new suggestion:
{'x': array([[0.00218701]]), 'f_a(x)': array([-0.46971535]), 'opt_obj': None}
Acquisition function target probability and method local
led to new suggestion:
{'x': array([[0.85715855]]), 'f_a(x)': array([-0.5]), 'opt_obj': None}
Acquisition function target probability and method hgdl
[[0.71750198]
[0.66034488]
[0.93660211]] [0 1 2]
led to new suggestion:
{'x': array([[0.71750198]]), 'f_a(x)': array([-0.5]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f39143c6390>}
Acquisition function total correlation and method global
led to new suggestion:
{'x': array([[0.9463929]]), 'f_a(x)': array([-4.62400106]), 'opt_obj': None}
Acquisition function total correlation and method local
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([-3.23830055]), 'opt_obj': None}
Acquisition function total correlation and method hgdl
[[1.]
[1.]
[0.]] [0 2]
led to new suggestion:
{'x': array([[1.]]), 'f_a(x)': array([-3.23830055]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38f40d3f10>}
#here we can test other options of the ask() command
bounds = np.array([[0.0,1.0]])
new_suggestion = my_gpo.ask(bounds, acquisition_function="total_correlation", method="global",
max_iter=10, n = 5, info = True)
my_gpo.ask(bounds, n = 5, acquisition_function="variance", vectorized=True, method = 'global')
my_gpo.ask(bounds, n = 1, acquisition_function="relative information entropy", vectorized=True, method = 'global')
my_gpo.ask(bounds, n = 2, acquisition_function="expected improvement", vectorized=True, method = 'global')
my_gpo.ask(bounds, n = 1, acquisition_function="variance", vectorized=True, method = 'global')
my_gpo.ask(bounds, n = 3, acquisition_function="variance", vectorized=True, method = 'hgdl', dask_client=client)
print(new_suggestion)
/home/marcus/Coding/gpCAM/gpcam/gp_optimizer_base.py:426: UserWarning: You specified n>1 and method != 'hgdl' in ask(). The acquisition function has therefore been changed to 'total correlation'.
warnings.warn("You specified n>1 and method != 'hgdl' in ask(). The acquisition function "
differential_evolution step 1: f(x)= 23.826817006151664
differential_evolution step 2: f(x)= 23.826817006151664
differential_evolution step 3: f(x)= 23.826817006151664
differential_evolution step 4: f(x)= 23.826817006151664
differential_evolution step 5: f(x)= 23.826817006151664
differential_evolution step 6: f(x)= 23.826817006151664
differential_evolution step 7: f(x)= 23.826817006151664
differential_evolution step 8: f(x)= 23.826817006151664
differential_evolution step 9: f(x)= 23.826817006151664
differential_evolution step 10: f(x)= 23.826817006151664
[[0. ]
[0. ]
[0.56647288]
[0.35067383]] [0 2 3]
{'x': array([[0.16745676],
[0.02410761],
[0.8503915 ],
[0.99319883],
[0.76394347]]), 'f_a(x)': array([-23.82681701]), 'opt_obj': None}
#we can evaluate the acqisiiton function on batches of candidates in parallel:
candidates = np.random.uniform(low = bounds[:,0], high=bounds[:,1], size = (30,1))
candidate_list = [entry for entry in candidates]
#ask sequentially
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=False)["x"][0])
#ask in parallel on DASK workers, but sequentially on each worker:
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=False, batch_size = 10, dask_client=client)["x"][0])
#ask in parallel on DASK workers, and vectorized (if possible) on each worker:
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=True, batch_size = 10, dask_client=client)["x"][0])
#ask vectorized (if possible):
print("suggestions=", my_gpo.ask(candidate_list, n = 30, acquisition_function="variance", vectorized=True)["x"][0])
print("They should be the same!")
suggestions= [0.0001488]
suggestions= [0.0001488]
suggestions= [0.0001488]
suggestions= [0.0001488]
They should be the same!
bounds = np.array([[0.0,1.0]])
#You can even start an ask() search asynchronously and check back later what was found
new_suggestion = my_gpo.ask(bounds, acquisition_function=acquisition_functions[0], method="hgdlAsync", dask_client=client)
time.sleep(10)
print(new_suggestion)
new_suggestion["opt_obj"].kill_client()
{'x': array([[0.]]), 'f_a(x)': array([-0.]), 'opt_obj': <hgdl.hgdl.HGDL object at 0x7f38f4117010>}
[{'x': array([0.16479708]),
'f(x)': np.float64(-0.10366144850444385),
'classifier': 'minimum',
'Hessian eigvals': array([18.12576478]),
'df/dx': array([8.39661674e-07]),
'|df/dx|': np.float64(8.396616735240059e-07),
'radius': np.float64(0.055170085910737195)},
{'x': array([0.56647287]),
'f(x)': np.float64(-0.0848648310545249),
'classifier': 'minimum',
'Hessian eigvals': array([21.50722656]),
'df/dx': array([-6.80289158e-08]),
'|df/dx|': np.float64(6.802891583390647e-08),
'radius': np.float64(0.04649599972507205)},
{'x': array([0.65342379]),
'f(x)': np.float64(-0.0757817303393714),
'classifier': 'minimum',
'Hessian eigvals': array([8.66123839]),
'df/dx': array([-5.50851031e-07]),
'|df/dx|': np.float64(5.508510314555792e-07),
'radius': np.float64(0.11545693059202564)}]