In [ ]:

    

import copy
import pickle
import numpy as np
import sklearn.gaussian_process as gp
import matplotlib.pyplot as plt


    

In [ ]:

    

import sys
sys.path.append('..')
from optimisation.bayesian_utils import store_GP, restore_GP


    

In [ ]:

    

np.random.seed(100)


    

In [ ]:

    

f = lambda x: np.exp(-(x - 2)**2) + np.exp(-(x - 6)**2/10) + 1/ (x**2 + 1)
x_min = -2
x_max = 10

xs = np.linspace(x_min, x_max, 100)
ys = f(xs)

rxs = np.random.uniform(x_min, x_max, size=(6,1))
rys = f(rxs)


    

In [ ]:

    

plt.figure(figsize=(16,6))
plt.plot(xs, ys, 'g-')
plt.plot(rxs, rys, 'bo')
plt.show()


    

In [ ]:

    

gp_params = dict(
    alpha = 1e-10, # larger => more noise. Default = 1e-10
    # nu=1.5 assumes the target function is once-differentiable
    kernel = 1.0 * gp.kernels.Matern(nu=1.5),# + gp.kernels.WhiteKernel(),
    #kernel = 1.0 * gp.kernels.RBF(),
    n_restarts_optimizer = 10,
    # make the mean 0 (theoretically a bad thing, see docs, but can help)
    # with the constant offset in the kernel this shouldn't be required
    #normalize_y = True,
    copy_X_train = True # whether to make a copy of the training data (in-case it is modified)
)

gp_model = gp.GaussianProcessRegressor(**gp_params)


    

In [ ]:

    

gp_model.fit(rxs, rys)


    

In [ ]:

    

gp_mu, gp_sigma = gp_model.predict(xs.reshape(-1,1), return_std=True)
gp_sigma = gp_sigma.reshape(-1,1)


    

In [ ]:

    

plt.figure(figsize=(16,6))
plt.plot(xs, ys, 'g-')
plt.plot(rxs, rys, 'bo')
plt.plot(xs, gp_mu, 'r-')
plt.fill_between(xs, gp_mu.flatten()-2*gp_sigma.flatten(), gp_mu.flatten()+2*gp_sigma.flatten(), alpha=0.3)
#plt.plot(xs, sample.flatten(), 'bo', color='black')
#plt.plot(xs2, sample2.flatten(), '-', color='black')
plt.show()


    

In [ ]:

    

def GPs_equal(a, b):
    ad, bd = a.__dict__, b.__dict__
    if ad.keys() != bd.keys():
        print('dicts different')
        return False
    for k in ad.keys():
        if k in {'rng'}: # exclude these items
            continue
        x, y = ad[k], bd[k]
        if isinstance(x, np.ndarray):
            if not np.all(x == y):
                print('element {} is different:\n{}\n{}'.format(k, x, y))
                return False
        else:
            if x != y:
                print('element {} is different:\n{}\n{}'.format(k, x, y))
                return False
    return True


    

In [ ]:

    

gp_before = copy.deepcopy(gp_model)
gp_pickled = pickle.loads(pickle.dumps(gp_model))


    

In [ ]:

    

gp_stored = store_GP(gp_model)
gp_stored


    

In [ ]:

    

restored = restore_GP(gp_stored, gp_params, rxs, rys)


    

In [ ]:

    

if GPs_equal(gp_pickled, gp_model):
    print('GP pickled == GP after')
else:
    print('NOT EQUAL')
    
if GPs_equal(gp_before, gp_model):
    print('GP before == GP after')
else:
    print('NOT EQUAL')
    
if GPs_equal(gp_model, restored):
    print('original GP == restored GP')
else:
    print('NOT EQUAL')


    

In [ ]:

    

gp2_mu, gp2_sigma = restored.predict(xs.reshape(-1,1), return_std=True)
gp2_sigma = gp2_sigma.reshape(-1,1)


    

In [ ]:

    

np.all(gp2_mu == gp_mu)


    

In [ ]:

    

np.all(gp2_sigma == gp_sigma)


    

In [ ]:

    

plt.figure(figsize=(16,6))
plt.plot(xs, ys, 'g-')
plt.plot(rxs, rys, 'bo')
plt.plot(xs, gp_mu, 'r-')
plt.fill_between(xs, gp_mu.flatten()-2*gp_sigma.flatten(), gp_mu.flatten()+2*gp_sigma.flatten(), hatch='/', alpha=0.1)
plt.fill_between(xs, gp2_mu.flatten()-2*gp2_sigma.flatten(), gp2_mu.flatten()+2*gp2_sigma.flatten(), hatch='\\', alpha=0.2)
#plt.plot(xs, sample.flatten(), 'bo', color='black')
#plt.plot(xs2, sample2.flatten(), '-', color='black')
plt.show()


    

In [ ]: