Single-Task Test#
#install the right fvgp version
#!pip install fvgp~=4.8.6
Setup#
import numpy as np
import matplotlib.pyplot as plt
from fvgp import GP
import time
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#
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(200).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 0x7f892047d4d0>
Customizing a Gaussian Process#
from fvgp.kernels import *
from scipy import sparse
def my_noise(x,hps):
#This is a simple noise function, but can be arbitrarily complex using many hyperparameters.
#The noise can be a vector, a matrix, or a sparse matrix in case gp2Scale is used.
return np.zeros(len(x)) + hps[2]
#stationary
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 ios a simple mean function but it can be arbitrarily complex using many hyperparameters.
return 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 0x7f892049d9d0>]
Initialization and different training options#
st = time.time()
from loguru import logger
logger.disable("fvgp")
my_gp1 = GP(x_data,y_data,
init_hyperparameters = np.ones((2))/10., # we need enough of those for kernel, noise, and prior mean functions
noise_variances=np.ones(y_data.shape) * 0.01, # providing noise variances and a noise function will raise a warning
compute_device='cpu',
kernel_function=skernel,
#kernel_function_grad=None,
#prior_mean_function=meanf,
#prior_mean_function_grad=None,
#noise_function=my_noise,
gp2Scale = False,
#linalg_mode='Inv',
ram_economy=True,
)
print("Initial likelihood: ", my_gp1.log_likelihood())
hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
[0.01,10.], #length scale for the kernel
#[0.001,0.1], #noise
#[0.01,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_gp1.update_gp_data(x_update,
y_update,
noise_variances_new=np.ones(y_update.shape) * 0.05,
append=True, rank_n_update=True)
print("Standard Training (MCMC)")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, info = False)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
print("Bayesian optimization Training")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, method='bo', max_iter=20)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
print("ADAM")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100, method="adam")
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
print("Global Training")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, method='global', max_iter = 20)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
print("Local Training")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, method='local')
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
print("HGDL Training")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, method='hgdl', max_iter=2, dask_client=client)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("ML: ",my_gp1.log_likelihood())
print("")
Initial likelihood: -213.1740936362793
Standard Training (MCMC)
Result= [0.76424892 0.0527507 ] after 20.727948904037476 seconds
ML: 5.06092122839587
Bayesian optimization Training
Result= [0.89956693 0.0509623 ] after 22.176537036895752 seconds
ML: 5.91214654554571
ADAM
fvGP adam iteration 1 out of 100: f(x)= -5.91214654554571, |grad|= 26.74986798447284
fvGP adam iteration 10 out of 100: f(x)= -5.359152300352463, |grad|= 242.87238573555769
fvGP adam iteration 20 out of 100: f(x)= -5.921681444400349, |grad|= 5.716186345237661
fvGP adam iteration 30 out of 100: f(x)= -5.863242734971095, |grad|= 73.61599765159376
fvGP adam iteration 40 out of 100: f(x)= -5.8862360252392705, |grad|= 57.49494172920567
fvGP adam iteration 50 out of 100: f(x)= -5.910705237080009, |grad|= 32.36911565098958
fvGP adam iteration 60 out of 100: f(x)= -5.9181978701654145, |grad|= 18.804750403416932
fvGP adam iteration 70 out of 100: f(x)= -5.920405214054995, |grad|= 12.291675326119792
fvGP adam iteration 80 out of 100: f(x)= -5.9214052817576714, |grad|= 7.70633963085418
fvGP adam iteration 90 out of 100: f(x)= -5.921932539080814, |grad|= 3.2280779898193313
fvGP adam iteration 100 out of 100: f(x)= -5.922042901817321, |grad|= 0.3727405575221578
Result= [0.92111895 0.05079595] after 35.50975775718689 seconds
ML: 5.922001260056021
Global Training
Result= [0.92111895 0.05079595] after 52.05067777633667 seconds
ML: 5.922001260056675
Local Training
Result= [0.92111807 0.0508386 ] after 52.362624168395996 seconds
ML: 5.922044369128827
HGDL Training
Result= [0.92111802 0.05083854] after 54.367964029312134 seconds
ML: 5.922044369232509
#You can always test your gradient like this before running local optimizers
my_gp1.test_log_likelihood_gradient(np.array([1.,1.]), epsilon=1e-6)
(array([ 92.49749701, -262.43370701]), array([ 92.49749498, -262.43357188]))
More advanced: 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_gp1.stop_training(opt_obj)’
HGDL#
my_gp1.set_hyperparameters(np.array([1.,1.]))
print(my_gp1.hyperparameters)
opt_obj = my_gp1.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method='hgdl')
# The result won't change much (or at all) since this is such a simple optimization
for i in range(20):
my_gp1.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gp1.hyperparameters)
time.sleep(0.1)
my_gp1.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1.]
iteration 0 : [1. 1.]
iteration 1 : [1. 1.]
/home/marcus/Coding/fvGP/fvgp/gp.py:1276: UserWarning: Hyperparameter update not successful len(optima list) = 0
hps = self.trainer.update_hyperparameters(opt_obj)
iteration 2 : [1. 1.]
iteration 3 : [1. 1.]
iteration 4 : [1. 1.]
iteration 5 : [1. 1.]
iteration 6 : [0.92097494 0.05083568]
iteration 7 : [0.92097494 0.05083568]
iteration 8 : [0.92097494 0.05083568]
iteration 9 : [0.92097494 0.05083568]
iteration 10 : [0.92097494 0.05083568]
iteration 11 : [0.92097494 0.05083568]
iteration 12 : [0.92097494 0.05083568]
iteration 13 : [0.92097494 0.05083568]
iteration 14 : [0.92097494 0.05083568]
iteration 15 : [0.92097494 0.05083568]
iteration 16 : [0.92097494 0.05083568]
iteration 17 : [0.92097494 0.05083568]
iteration 18 : [0.92097494 0.05083568]
iteration 19 : [0.92097494 0.05083568]
ADAM#
my_gp1.set_hyperparameters(np.array([1.,1.,]))
print(my_gp1.hyperparameters)
opt_obj = my_gp1.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_gp1.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gp1.hyperparameters)
time.sleep(0.1)
my_gp1.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1.]
iteration 0 : [1. 1.]
iteration 1 : [1.02995373 0.97002198]
iteration 2 : [1.07918433 0.92039027]
iteration 3 : [1.1267212 0.87158198]
iteration 4 : [1.18043943 0.8146615 ]
iteration 5 : [1.22183328 0.76892938]
iteration 6 : [1.25991524 0.72488246]
iteration 7 : [1.29470529 0.68248928]
iteration 8 : [1.32641284 0.64156345]
iteration 9 : [1.35538533 0.60176299]
iteration 10 : [1.38206675 0.56257972]
iteration 11 : [1.4069748 0.52330586]
iteration 12 : [1.43070364 0.4829574 ]
iteration 13 : [1.45396385 0.44011846]
iteration 14 : [1.48259568 0.38232212]
iteration 15 : [1.50876588 0.32466938]
iteration 16 : [1.53968842 0.25298503]
iteration 17 : [1.57872597 0.16461964]
iteration 18 : [1.62171292 0.06716527]
iteration 19 : [1.64825418 0.02796812]
MCMC#
my_gp1.set_hyperparameters(np.array([1.,1.]))
print(my_gp1.hyperparameters)
opt_obj = my_gp1.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_gp1.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gp1.hyperparameters)
time.sleep(0.1)
my_gp1.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1.]
iteration 0 : [1. 1.]
iteration 1 : [0.48343583 0.03623791]
iteration 2 : [0.48343583 0.03623791]
iteration 3 : [0.78621225 0.04968284]
iteration 4 : [0.84217543 0.05161471]
iteration 5 : [0.92171792 0.0456737 ]
iteration 6 : [0.97339672 0.05351485]
iteration 7 : [0.94858155 0.0434971 ]
iteration 8 : [0.88565774 0.04844877]
iteration 9 : [0.87348956 0.05260795]
iteration 10 : [0.92558521 0.04948731]
iteration 11 : [0.92558521 0.04948731]
iteration 12 : [0.95630491 0.05351459]
iteration 13 : [0.94708566 0.0487519 ]
iteration 14 : [0.99992396 0.04681474]
iteration 15 : [1.02223552 0.05179775]
iteration 16 : [1.03888557 0.04876278]
iteration 17 : [1.00950004 0.05099945]
iteration 18 : [0.88520292 0.05119419]
iteration 19 : [0.88520292 0.05119419]
BO#
my_gp1.set_hyperparameters(np.array([1.,1.]))
print(my_gp1.hyperparameters)
opt_obj = my_gp1.train(hyperparameter_bounds=hps_bounds, dask_client=client, asynchronous=True, method='bo')
# The result won't change much (or at all) since this is such a simple optimization
for i in range(20):
my_gp1.update_hyperparameters(opt_obj)
print("iteration ", i, " : ",my_gp1.hyperparameters)
time.sleep(0.1)
my_gp1.stop_training(opt_obj) ##this leaves the dask client alive, kill_client() will shut it down.
[1. 1.]
iteration 0 : [1. 1.]
iteration 1 : [0.33328893 0.03835522]
iteration 2 : [0.48182468 0.0370138 ]
iteration 3 : [0.48182468 0.0370138 ]
iteration 4 : [0.48182468 0.0370138 ]
iteration 5 : [0.48182468 0.0370138 ]
iteration 6 : [0.77216794 0.0436141 ]
iteration 7 : [1.18671927 0.05501523]
iteration 8 : [0.98966314 0.04929904]
iteration 9 : [0.98966314 0.04929904]
iteration 10 : [0.86001917 0.04967969]
iteration 11 : [0.86001917 0.04967969]
iteration 12 : [0.86001917 0.04967969]
iteration 13 : [0.92885257 0.050797 ]
iteration 14 : [0.92885257 0.050797 ]
iteration 15 : [0.92885257 0.050797 ]
iteration 16 : [0.92885257 0.050797 ]
iteration 17 : [0.92885257 0.050797 ]
iteration 18 : [0.91440511 0.05056518]
iteration 19 : [0.91440511 0.05056518]
The Result#
#let's make a prediction
x_pred = np.linspace(0,1,1000)
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, info = False)
# different ways to call
var1 = my_gp1.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=False)["v(x)"]
var1 = my_gp1.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=True)["v(x)"]
mean1 = my_gp1.posterior_mean(x_pred.reshape(-1,1))["m(x)"]
var1 = my_gp1.posterior_covariance(x_pred.reshape(-1,1), variance_only=False, add_noise=True)["v(x)"]
mean_grad = my_gp1.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_gp1.x_data,my_gp1.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_gp1.rmse(x_pred1D,f1(x_pred1D).flatten()))
print("NRMSE: ",my_gp1.nrmse(x_pred1D,f1(x_pred1D).flatten()))
print("CRPS (mean, std): ",my_gp1.crps(x_pred1D,f1(x_pred1D).flatten()))
print("R2: ",my_gp1.r2(x_pred1D,f1(x_pred1D).flatten()))
print("NLPD: ",my_gp1.nlpd(x_pred1D,f1(x_pred1D).flatten()))
print("MSLL: ",my_gp1.msll(x_pred1D,f1(x_pred1D).flatten()))
print("MAPE: ",my_gp1.mape(x_pred1D,f1(x_pred1D).flatten()))
print("INTERVAL SCORE: ",my_gp1.interval_score(x_pred1D,f1(x_pred1D).flatten()))
print("MPIW: ",my_gp1.mpiw(x_pred1D))
print("PICP: ",my_gp1.picp(x_pred1D,f1(x_pred1D).flatten()))
print("Coverage Curve:")
cov_curve = my_gp1.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_gp1.plot_observed_vs_predicted(x_pred1D,f1(x_pred1D).flatten())
Posterior Mean and Uncertainty
Posterior Mean Gradient
RMSE: 0.09260647039789782
NRMSE: 0.023480912301903974
CRPS (mean, std): (np.float64(0.053575710460306916), np.float64(0.037424610587455144))
R2: 0.9919448686246046
NLPD: -0.9159794406460198
MSLL: -2.380091675715631
MAPE: 0.858495859845033
INTERVAL SCORE: 0.5106077573570535
MPIW: 0.5106077573570535
PICP: 1.0
Coverage Curve:
predicted vs. observed
Predicted Information Gain#
relative_entropy = my_gp1.gp_relative_information_entropy_set(x_pred.reshape(-1,1))["RIE"]
plt.figure(figsize = (16,10))
plt.plot(x_pred,relative_entropy, label = "relative_entropy", linewidth = 4)
plt.scatter(x_data,y_data, color = 'black')
plt.legend()
<matplotlib.legend.Legend at 0x7f89000d5890>
#We can ask mutual information and total correlation there is given some test data
x_test = np.array([[0.45],[0.45]])
print("MI: ",my_gp1.gp_mutual_information(x_test))
print("TC: ",my_gp1.gp_total_correlation(x_test))
my_gp1.gp_entropy(x_test)
my_gp1.gp_entropy_grad(x_test, 0)
my_gp1.gp_kl_div(x_test, np.ones((len(x_test))), np.identity((len(x_test))))
my_gp1.gp_relative_information_entropy(x_test)
my_gp1.gp_relative_information_entropy_set(x_test)
my_gp1.posterior_covariance(x_test)
my_gp1.posterior_covariance_grad(x_test)
my_gp1.posterior_mean(x_test)
my_gp1.posterior_mean_grad(x_test)
my_gp1.posterior_probability(x_test, np.ones((len(x_test))), np.identity((len(x_test))))
MI: {'x': array([[0.45],
[0.45]]), 'mutual information': np.float64(4.228743563495186)}
TC: {'x': array([[0.45],
[0.45]]), 'total correlation': np.float64(14.374505869618503)}
{'mu': array([0.70488351, 0.70488351]),
'covariance': array([[0.01503793, 0.00455779],
[0.00455779, 0.01503793]]),
'probability': np.float64(0.14343483888312153)}
Running many GPs at once in parallel#
#duplicate data: in practice, this would be different data in every column
y_data = np.broadcast_to(y_data[:, None], (y_data.size, 10))
my_gp1 = GP(x_data,y_data,
init_hyperparameters = np.ones((2))/10., # we need enough of those for kernel, noise, and prior mean functions
noise_variances=np.ones(y_data.shape[0]) * 0.1, # providing noise variances and a noise function will raise a warning
compute_device='cpu',
)
hps_bounds = np.array([[0.01,10.], #signal variance for the kernel
[0.01,10.], #length scale for the kernel
])
print("Standard Training (MCMC)")
hps = my_gp1.train(hyperparameter_bounds=hps_bounds, info = True, max_iter = 100)
print("Result=", hps, "after ", time.time() - st, " seconds")
print("")
Standard Training (MCMC)
Starting likelihood. f(x)= -59.847042329257675
Finished 10 out of 100 iterations. f(x)= -59.847042329257675
Finished 20 out of 100 iterations. f(x)= -35.28901444862643
Finished 30 out of 100 iterations. f(x)= -35.06830309494018
Finished 40 out of 100 iterations. f(x)= -33.108443868916595
Finished 50 out of 100 iterations. f(x)= -32.68669805235254
Finished 60 out of 100 iterations. f(x)= -35.15067991134816
Finished 70 out of 100 iterations. f(x)= -32.96405020660424
Finished 80 out of 100 iterations. f(x)= -33.164071729452104
Finished 90 out of 100 iterations. f(x)= -32.83160745938736
Result= [1.81417019 0.29047259] after 102.68088722229004 seconds
x_pred = np.linspace(0,1,1000).reshape(1000,1)
mean = my_gp1.posterior_mean(x_pred)["m(x)"]
sd = np.sqrt(my_gp1.posterior_covariance(x_pred)["v(x)"])
print("Posterior Means")
plt.figure(figsize = (16,10))
for i in range(10):
plt.plot(x_pred.flatten(),mean[:,i], label = "posterior mean", linewidth = 4)
plt.scatter(my_gp1.x_data,my_gp1.y_data[:,0], color = 'black')
plt.plot(x_pred1D,f1(x_pred1D), label = "latent function", linewidth = 4)
plt.show()
Posterior Means
