Multi-Task Test#
At first we have to install the newest version of fvGP
##First, install the newest version of fvgp
#!pip install fvgp~=4.8.6
#!pip install plotly
#pip install --upgrade kaleido
Setup#
import numpy as np
import matplotlib.pyplot as plt
from fvgp import GP
import plotly.graph_objects as go
from itertools import product
import plotly.io as pio
pio.renderers.default = "png"
%load_ext autoreload
%autoreload 2
Simple 1d Example#
Data#
def f1(x): return 0.5 * x
def f2(x): return (-.25 * x) - 1.
x_pred1d = np.linspace(0,1,50)
plt.plot(x_pred1d,f1(x_pred1d))
plt.plot(x_pred1d,f2(x_pred1d))
x_data = np.random.rand(10)
y_data1 = f1(x_data) + np.random.uniform(low = -0.01, high = 0.01, size =len(x_data))
y_data2 = f2(x_data) + np.random.uniform(low = -0.01, high = 0.01, size =len(x_data))
plt.scatter(x_data,y_data1)
plt.scatter(x_data,y_data2)
plt.show()
GP initialization#
from fvgp import fvGP
my_gp2 = fvGP(x_data.reshape(len(x_data),1), np.column_stack([y_data1, y_data2]))
print("Global Training in progress")
my_gp2.train(max_iter = 20, info=True)
Global Training in progress
Starting likelihood. f(x)= 31.507854514037824
Finished 10 out of 20 iterations. f(x)= 43.90229304384373
/home/marcus/Coding/fvGP/fvgp/gp.py:500: UserWarning: Hyperparameters initialized to a vector of ones.
warnings.warn("Hyperparameters initialized to a vector of ones.")
/home/marcus/Coding/fvGP/fvgp/gp.py:536: UserWarning: No noise function or measurement noise provided. Noise variances will be set to (0.01 * mean(|y_data|))^2.
self.likelihood = GPlikelihood(self.data,
/home/marcus/Coding/fvGP/fvgp/gp.py:1021: UserWarning: Default hyperparameter_bounds initialized because none were provided. This will fail for custom kernel, mean, or noise functions
warnings.warn("Default hyperparameter_bounds initialized because none were provided. "
array([0.92849284, 6.60077476, 1.45475153])
print("Data as seen by fvgp:")
print(my_gp2.fvgp_x_data)
print(my_gp2.fvgp_y_data)
print("Transformed data for the GP:")
print(my_gp2.x_data)
print(my_gp2.y_data)
Data as seen by fvgp:
[[0.55496274]
[0.54980083]
[0.73855146]
[0.51203224]
[0.22668151]
[0.09026322]
[0.01267208]
[0.59715033]
[0.14257962]
[0.1684044 ]]
[[ 0.27415257 -1.13294602]
[ 0.28336197 -1.13220687]
[ 0.36511105 -1.17771369]
[ 0.24927794 -1.1323592 ]
[ 0.1051445 -1.05950872]
[ 0.05075115 -1.03151234]
[ 0.00370759 -0.99451728]
[ 0.29506668 -1.15926706]
[ 0.06531947 -1.02623569]
[ 0.08211646 -1.03361486]]
Transformed data for the GP:
[[0.55496274 0. ]
[0.54980083 0. ]
[0.73855146 0. ]
[0.51203224 0. ]
[0.22668151 0. ]
[0.09026322 0. ]
[0.01267208 0. ]
[0.59715033 0. ]
[0.14257962 0. ]
[0.1684044 0. ]
[0.55496274 1. ]
[0.54980083 1. ]
[0.73855146 1. ]
[0.51203224 1. ]
[0.22668151 1. ]
[0.09026322 1. ]
[0.01267208 1. ]
[0.59715033 1. ]
[0.14257962 1. ]
[0.1684044 1. ]]
[[ 0.27415257]
[ 0.28336197]
[ 0.36511105]
[ 0.24927794]
[ 0.1051445 ]
[ 0.05075115]
[ 0.00370759]
[ 0.29506668]
[ 0.06531947]
[ 0.08211646]
[-1.13294602]
[-1.13220687]
[-1.17771369]
[-1.1323592 ]
[-1.05950872]
[-1.03151234]
[-0.99451728]
[-1.15926706]
[-1.02623569]
[-1.03361486]]
Predictions#
y_test = np.hstack([f1(x_pred1d).reshape(len(x_pred1d),1),f2(x_pred1d).reshape(len(x_pred1d),1)])
print("RMSE: ", my_gp2.rmse(x_pred1d.reshape(len(x_pred1d),1),y_test))
print("NLPD: ", my_gp2.nlpd(x_pred1d.reshape(len(x_pred1d),1),y_test))
print("R2: ", my_gp2.r2(x_pred1d.reshape(len(x_pred1d),1),y_test))
print("CRPS (mean, std): ", my_gp2.crps(x_pred1d.reshape(len(x_pred1d),1),y_test))
RMSE: 0.010609817928863252
NLPD: -3.8817623356303104
R2: 0.999768477462133
CRPS (mean, std): (np.float64(0.004232198868484882), np.float64(0.005020593676119593))
#mean and standard deviation
mean = my_gp2.posterior_mean(x_pred=x_pred1d.reshape(50,1))["m(x)"]
std = np.sqrt(my_gp2.posterior_covariance(x_pred=x_pred1d.reshape(50,1), x_out=np.array([0,1]))["v(x)"])
plt.plot(x_pred1d.reshape(50,1),mean[:,0], label = "mean task 1")
plt.plot(x_pred1d.reshape(50,1),mean[:,1], label = "mean task 2")
plt.scatter(x_data,y_data1)
plt.scatter(x_data,y_data2)
plt.plot(x_pred1d,f1(x_pred1d), label = "task 1 ground truth")
plt.plot(x_pred1d,f2(x_pred1d), label = "task 2 ground truth")
plt.fill_between(x_pred1d, mean[:,0] - 3. * std[:,0], mean[:,0] + 3. * std[:,0], alpha = 0.5, color = "grey")
plt.fill_between(x_pred1d, mean[:,1] - 3. * std[:,1], mean[:,1] + 3. * std[:,1], alpha = 0.5, color = "grey")
plt.legend()
plt.show()
#mean gradient and variance gradient
mean_grad = my_gp2.posterior_mean_grad(x_pred=x_pred1d.reshape(50,1), x_out=np.array([0,1]))["dm/dx"]
var_grad = my_gp2.posterior_covariance_grad(x_pred=x_pred1d.reshape(50,1), x_out=np.array([0,1]))["dv/dx"]
plt.plot(x_pred1d.reshape(50,1),mean_grad[:,0,0], label = "mean gradient task 1")
plt.plot(x_pred1d.reshape(50,1),mean_grad[:,0,1], label = "mean gradient task 2")
plt.plot(x_pred1d,np.gradient(f1(x_pred1d), 1./50.), label = "grad task 1 ground truth")
plt.plot(x_pred1d,np.gradient(f2(x_pred1d), 1./50.), label = "grad task 2 ground truth")
plt.plot(x_pred1d.reshape(50,1),var_grad[:,0,0], label = "var gradient task 1")
plt.plot(x_pred1d.reshape(50,1),var_grad[:,0,1], label = "var gradient task 2")
plt.legend()
plt.show()
What if some tasks are missing from the data#
It works just fine, but we have to insert np.nan at positions of missing data, both for y and the variances.#
y_data = np.column_stack([y_data1, y_data2])
noise_variances = np.zeros(y_data.shape) + 0.01
y_data[2,0] = np.nan
noise_variances[2,0] = np.nan
y_data[6,1] = np.nan
noise_variances[6,1] = np.nan
my_gp2 = fvGP(x_data.reshape(len(x_data),1), y_data, noise_variances=noise_variances)
print("Global Training in progress")
my_gp2.train(max_iter = 20)
Global Training in progress
array([0.12779334, 1.76233351, 0.05770436])
3d Example#
Data#
data = np.load("./data/sim_variable_mod.npy")
sparsification = 4
x_data3 = data[:,5:][::sparsification]
y_data3 = data[:,0:2][::sparsification]
#it is good practice to check the format of the data
print(x_data3.shape)
print(y_data3.shape)
(1583, 3)
(1583, 2)
index = np.where(x_data3[:,2] == 1200.)[0]
x_data3=x_data3[index,0:2]
y_data3=y_data3[index]
for i in range(x_data3.shape[1]):
x_data3[:,i] = x_data3[:,i] - np.min(x_data3[:,i])
x_data3[:,i] = x_data3[:,i] / np.max(x_data3[:,i])
x = np.linspace(0,1,100)
y = np.linspace(0,1,100)
x_pred3D = np.asarray(list(product(x, y)))
def scatter(x,y,z,size=3, color = 1):
#if not color: color = z
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=x, y=y, z=z,mode='markers',marker=dict(color=color, size = size)))
fig.update_layout(autosize=False,
width=800, height=800,
font=dict(size=18,),
margin=dict(l=0, r=0, b=0, t=0))
fig.show()
scatter(x_data3[:,0],x_data3[:,1],y_data3[:,0], size = 5, color = y_data3[:,0])
scatter(x_data3[:,0],x_data3[:,1],y_data3[:,1], size = 5, color = y_data3[:,1])
Initialization#
(a) Default behavior — minimal#
from fvgp import fvGP
my_gp2 = fvGP(x_data3,y_data3)
print("Global Training in progress")
my_gp2.train(max_iter = 2)
Global Training in progress
/home/marcus/Coding/fvGP/fvgp/gp.py:500: UserWarning: Hyperparameters initialized to a vector of ones.
warnings.warn("Hyperparameters initialized to a vector of ones.")
/home/marcus/Coding/fvGP/fvgp/gp.py:536: UserWarning: No noise function or measurement noise provided. Noise variances will be set to (0.01 * mean(|y_data|))^2.
self.likelihood = GPlikelihood(self.data,
/home/marcus/Coding/fvGP/fvgp/gp.py:1021: UserWarning: Default hyperparameter_bounds initialized because none were provided. This will fail for custom kernel, mean, or noise functions
warnings.warn("Default hyperparameter_bounds initialized because none were provided. "
array([0.16782057, 5.9472849 , 9.11743636, 4.47873019])
(b) Custom kernel#
It is vital in the multi-task case to think hard about kernel design. The kernel is now a function over X x X x T x T, where X is the input and T is the output space. Print the input of the kernel, it will have the dimensionality of this cartesian product space. The default kernel in fvgp is just a Matern kernel operating in this new space.
#A simple kernel that won't lead to good performance because it's stationary
from fvgp.kernels import *
def mkernel(x1,x2,hps):
d = get_distance_matrix(x1,x2)
return hps[0] * matern_kernel_diff1(d,hps[1])
my_gp2 = fvGP(x_data3,y_data3,
init_hyperparameters=np.ones((2)), kernel_function=mkernel
)
print("MCMC Training in progress")
bounds = np.array([[0.01,1.],[0.01,1.]])
my_gp2.train(hyperparameter_bounds=bounds,max_iter = 20)
MCMC Training in progress
array([0.87816069, 0.58556637])
(c) A custom deep kernel#
from fvgp.deep_kernel_network import *
from fvgp.kernels import *
iset_dim = 3
gp_deep_kernel_layer_width = 5
n = Network(iset_dim, gp_deep_kernel_layer_width)
print(n.number_of_hps)
def deep_multi_task_kernel(x1, x2, hps): # pragma: no cover
signal_var = hps[0]
length_scale = hps[1]
hps_nn = hps[2:]
w1_indices = np.arange(0, gp_deep_kernel_layer_width * iset_dim)
last = gp_deep_kernel_layer_width * iset_dim
w2_indices = np.arange(last, last + gp_deep_kernel_layer_width ** 2)
last = last + gp_deep_kernel_layer_width ** 2
w3_indices = np.arange(last, last + gp_deep_kernel_layer_width * iset_dim)
last = last + gp_deep_kernel_layer_width * iset_dim
b1_indices = np.arange(last, last + gp_deep_kernel_layer_width)
last = last + gp_deep_kernel_layer_width
b2_indices = np.arange(last, last + gp_deep_kernel_layer_width)
last = last + gp_deep_kernel_layer_width
b3_indices = np.arange(last, last + iset_dim)
n.set_weights(hps_nn[w1_indices].reshape(gp_deep_kernel_layer_width, iset_dim),
hps_nn[w2_indices].reshape(gp_deep_kernel_layer_width, gp_deep_kernel_layer_width),
hps_nn[w3_indices].reshape(iset_dim, gp_deep_kernel_layer_width))
n.set_biases(hps_nn[b1_indices].reshape(gp_deep_kernel_layer_width),
hps_nn[b2_indices].reshape(gp_deep_kernel_layer_width),
hps_nn[b3_indices].reshape(iset_dim))
x1_nn = n.forward(x1)
x2_nn = n.forward(x2)
d = get_distance_matrix(x1_nn, x2_nn)
k = signal_var * matern_kernel_diff1(d, length_scale)
return k
my_gp2 = fvGP(x_data3,y_data3,
init_hyperparameters=np.ones((n.number_of_hps+2))*0.1, kernel_function=deep_multi_task_kernel
)
print("MCMC Training in progress")
bounds = np.zeros((n.number_of_hps+2,2))
bounds[0] = np.array([0.01,1.])
bounds[1] = np.array([0.1,1.])
bounds[2:] = np.array([-1,1])
my_gp2.train(hyperparameter_bounds=bounds,max_iter = 300, method = "mcmc", info = True)
68
MCMC Training in progress
Starting likelihood. f(x)= -204924.16723757322
/home/marcus/Coding/fvGP/fvgp/gp.py:536: UserWarning: No noise function or measurement noise provided. Noise variances will be set to (0.01 * mean(|y_data|))^2.
self.likelihood = GPlikelihood(self.data,
Finished 10 out of 300 iterations. f(x)= -2014.3176487834087
Finished 20 out of 300 iterations. f(x)= -738.365549087808
Finished 30 out of 300 iterations. f(x)= -738.365549087808
Finished 40 out of 300 iterations. f(x)= -738.365549087808
Finished 50 out of 300 iterations. f(x)= -350.78515724388916
Finished 60 out of 300 iterations. f(x)= 62.80042041768576
Finished 70 out of 300 iterations. f(x)= 297.6517366030599
Finished 80 out of 300 iterations. f(x)= 351.78608632540823
Finished 90 out of 300 iterations. f(x)= 671.1634231780963
Finished 100 out of 300 iterations. f(x)= 856.3724746635505
Finished 110 out of 300 iterations. f(x)= 856.3724746635505
Finished 120 out of 300 iterations. f(x)= 900.0185671570468
Finished 130 out of 300 iterations. f(x)= 900.0185671570468
Finished 140 out of 300 iterations. f(x)= 902.1095666763365
Finished 150 out of 300 iterations. f(x)= 908.4551669613095
Finished 160 out of 300 iterations. f(x)= 908.4551669613095
Finished 170 out of 300 iterations. f(x)= 914.5271251575653
Finished 180 out of 300 iterations. f(x)= 927.6184482883712
Finished 190 out of 300 iterations. f(x)= 927.6184482883712
Finished 200 out of 300 iterations. f(x)= 932.9295086565519
Finished 210 out of 300 iterations. f(x)= 937.7466958752291
Finished 220 out of 300 iterations. f(x)= 937.7466958752291
Finished 230 out of 300 iterations. f(x)= 937.7466958752291
Finished 240 out of 300 iterations. f(x)= 937.7466958752291
Finished 250 out of 300 iterations. f(x)= 942.5304816928926
Finished 260 out of 300 iterations. f(x)= 942.5304816928926
Finished 270 out of 300 iterations. f(x)= 942.5304816928926
Finished 280 out of 300 iterations. f(x)= 942.5304816928926
Finished 290 out of 300 iterations. f(x)= 943.4724928195747
array([ 0.25437111, 0.10116515, -0.04178066, -0.6739599 , 0.38963764,
-0.11114264, -0.30660413, -0.27323933, -0.21130059, 0.22761933,
0.06489414, -0.36969247, 0.19846411, 0.61983991, 0.46212535,
0.21324786, -0.04935924, -0.1820413 , -0.78757663, 0.14636997,
0.02496616, 0.13003899, 0.07516544, 0.93365117, -0.0682337 ,
-0.09693454, 0.30227051, -0.06500306, 0.32356108, 0.10489028,
0.30111028, -0.51734482, 0.25470617, 0.36258342, 0.25175426,
0.42536017, 0.21826459, 0.18263785, 0.65670512, -0.27412586,
0.69986775, 0.95709718, 0.48286612, 0.10761125, 0.69169167,
0.23008218, -0.21289665, 0.02511654, 0.47261045, -0.28088074,
-0.0625802 , 0.40291708, -0.67251631, -0.1793064 , -0.33292806,
0.80448875, 0.54553336, -0.44180672, 0.11166209, 0.25168685,
0.1645877 , 0.43237351, 0.28575278, -0.28132167, 0.96462488,
0.52930974, 0.15588001, 0.80152634, -0.66054458, 0.60861895])
Prediction#
mean = my_gp2.posterior_mean(x_pred3D)["m(x)"]
var = my_gp2.posterior_covariance(x_pred3D)["v(x)"]
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=x_pred3D[:,0],y=x_pred3D[:,1], z=mean[:,0],
mode='markers',marker=dict(color=mean[:,0], size = 5)))
fig.add_trace(go.Scatter3d(x=x_data3[:,0], y=x_data3[:,1] , z=y_data3[:,0],
mode='markers',marker=dict(color=y_data3[:,0], size = 5)))
fig.update_layout(autosize=False,
width=800, height=800,
font=dict(size=18,),
margin=dict(l=0, r=0, b=0, t=0))
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=x_pred3D[:,0],y=x_pred3D[:,1], z=mean[:,1],
mode='markers',marker=dict(color=mean[:,0], size = 5)))
fig.add_trace(go.Scatter3d(x=x_data3[:,0], y=x_data3[:,1] , z=y_data3[:,1],
mode='markers',marker=dict(color=y_data3[:,1], size = 5)))
fig.update_layout(autosize=False,
width=800, height=800,
font=dict(size=18,),
margin=dict(l=0, r=0, b=0, t=0))
fig.show()
