# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# + [markdown] id="W7w7pNUH6K_9"
# ## Bayesian Optimisation Verification

# + id="c53z-Yqq6LAM" pycharm={"is_executing": true}
# import relevant packages
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernel, Matern
from sklearn.metrics import mean_squared_error
from itertools import product
from gp_para import gp_tuning
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer
from bayesian_opt_min import GPUCB, RGPUCB, EI
import GPy
import GPyOpt

# stop showing warning messages
def warn(*args, **kwargs):
    pass
import warnings
warnings.warn = warn

# + id="Ao0rGjXV6LAY"
### Impot data for training and evaluation
# Seed for reproducibility
rng = np.random.default_rng(12345)
# import data sheet (time:0~5000s) with temp 120, 140, 160, 180
address = 'data/degradation.xlsx'
y_experiment = np.empty(shape=[0,1])
i = 3
temp_list = np.array([120,140,160,180])
for temp in temp_list:
    df = pd.read_excel(address,sheet_name = 'normal data',usecols = [i],names = None,nrows = 5000)
    df = df.values.tolist()
    df = np.array(df)
    y_experiment =  np.append(y_experiment,df)
    i+=1

# + id="kqwnSBWD6LAZ"
### Data setup for training
# NUM_OF_DATA_POINTS must be an even number
NUM_OF_DATA_POINTS = 20
NUM_OF_PLOT_POINTS = 100
### Kernel setting and prediction
LOWEST_NOISE_LEVEL = 1e-16
HIGHEST_NOISE_LEVEL = 1e-12
KERNEL_ALPHA = 1e-15

# choosing data selection mode. See comments for what each mode does.
data_partition_mode = 1

# Mode 1: Data setup with half points evenly spread in 0-999, another half in 1000-4999
if data_partition_mode == 1:
    X_num = np.hstack([ np.linspace(0,999,num=int(3*NUM_OF_DATA_POINTS/4),dtype=int),
                        np.linspace(1000,4999,num=int(NUM_OF_DATA_POINTS/4),dtype=int)])

# Mode 2: Data setup with1/4 in 0-199, 2/4 in 200-999, 1/4 1000-4999
elif data_partition_mode == 2:
    X_num = np.hstack([ np.linspace(0,199,num=int(NUM_OF_DATA_POINTS/4),dtype=int),
                    np.linspace(200,999,num=int(NUM_OF_DATA_POINTS/2),dtype=int),
                    np.linspace(1000,4999,num=int(NUM_OF_DATA_POINTS/4),dtype=int)])

# Mode 3: Data setup with 2/5 in 0-99, 2/5 in 100-999, 1/5 1000-4999
elif data_partition_mode == 3: 
    X_num = np.hstack([ np.linspace(0,99,num=int(NUM_OF_DATA_POINTS*2/5),dtype=int),
                    np.linspace(199,999,num=int(NUM_OF_DATA_POINTS*2/5),dtype=int),
                    np.linspace(1000,4999,num=int(NUM_OF_DATA_POINTS*1/5),dtype=int)])

# Mode 4: Data setup with  points evenly spread 
elif data_partition_mode == 4:
    X_num = np.linspace(0,4999,num=int(NUM_OF_DATA_POINTS),dtype=int)

# Mode 5: Data setup with points spread in log space
else:
    X_num = np.logspace(0,np.log10(4999), num=NUM_OF_DATA_POINTS, base=10.0, dtype=int)
    HIGHEST_NOISE_LEVEL = 1e-15
    KERNEL_ALPHA = 8e-5

# Zipping training data
X_ = np.arange(0,5000).astype(int)
x_normal = np.empty(shape=[0,2])
temp_list = np.array([120,140,160,180])

for temp in temp_list:
    df_temp = np.full((5000,1),temp)
    x_normal_partial = np.hstack([X_.reshape(5000,1),df_temp])
    x_normal =  np.concatenate([x_normal,x_normal_partial])

Y = y_experiment[np.hstack([X_num,X_num+5000,X_num+5000*2,X_num+5000*3])]
X = np.empty(shape=[0,2])
for temp in temp_list:
    df_temp = np.full((int(NUM_OF_DATA_POINTS),1),temp)
    X_partial = np.hstack([X_num.reshape(int(NUM_OF_DATA_POINTS),1),df_temp])
    X = np.concatenate([X, X_partial])


# -

def plot(gp,X,temp_list,NUM_OF_PLOT_POINTS):
    X_5000 = np.arange(0,5000).astype(int)
    x1 = np.linspace(X[:,0].min(), X[:,0].max(),num = NUM_OF_PLOT_POINTS)
    x2 = np.linspace(X[:,1].min(), X[:,1].max(),num = NUM_OF_PLOT_POINTS)
    x1x2 = np.array(list(product(x1, x2)))
    y_pred, sigma = gp.predict(x1x2, return_std=True)
    X0p, X1p = x1x2[:,0].reshape(NUM_OF_PLOT_POINTS,NUM_OF_PLOT_POINTS), x1x2[:,1].reshape(NUM_OF_PLOT_POINTS,NUM_OF_PLOT_POINTS)
    Zp = np.reshape(y_pred,(NUM_OF_PLOT_POINTS,NUM_OF_PLOT_POINTS))
    sigma = np.reshape(sigma,(NUM_OF_PLOT_POINTS,NUM_OF_PLOT_POINTS))

    fig = plt.figure(figsize=(16,24),facecolor='white')
    ax = fig.add_subplot(3,2,1, projection='3d')
    surf = ax.plot_surface(X0p, X1p, Zp, rstride=1, cstride=1, cmap='jet', linewidth=0, antialiased=False)
    ax.scatter(X[:,0],X[:,1],Y)
    fig.colorbar(surf, aspect=5)

    ax = fig.add_subplot(3,2,2)
    pcm = ax.pcolormesh(X0p, X1p, Zp)
    fig.colorbar(pcm, ax=ax)

    i = 1
    j = int(X.size/8)
    k = 0
    plot_index = 3
    for temp in temp_list:
        plt.subplot(3,2,plot_index)
        plt.plot(x1,Zp[:,i],label = str(temp)+' Prediction',c='blue',lw=1,zorder=3)
        plt.scatter (X[j-int(X.size/8):j,0],Y[j-int(X.size/8):j],label = str(temp)+'Trained data',c='blue',zorder=2, s=3)
        plt.fill_between(x1, Zp[:,i] - 1.9600 * sigma[:,i],Zp[:,i] + 1.9600 * sigma[:,i],
                        alpha=0.5,color='grey',zorder=2, label='95% confidence interval')
        plt.scatter(X_5000,y_experiment[k:k+5000],label = str(temp)+' Experiment data',s=1,c='lightblue',zorder=1),
        plt.legend()
        plt.title('Plot for '+str(temp))
        i += int(NUM_OF_PLOT_POINTS/3-1)
        j += int(X.size/8)
        k+=5000
        plot_index+=1
    plt.show()


# Plot function to show confidence bounds 
def plot1D(x_normal,y_normal,X,y_mean, sigma,gp):
    y_mean = y_mean.reshape(y_mean.size)
    plt.figure(facecolor='white')
    plt.plot(x_normal, y_mean, c='blue', lw=1, zorder=4, label='Prediction')
    plt.scatter(X[:, 0], Y, c='r', s=5, zorder=3,label='Trained data')
    plt.fill_between(x_normal, y_mean - 1.9600 * sigma,y_mean + 1.9600 * sigma,
                        alpha=0.5,color='grey',zorder=2, label='95% confidence interval')
    plt.scatter(x_normal,y_normal,c='lightblue',s=1,zorder=1, label='Experiment data')
    plt.tick_params(axis='y')
    plt.tick_params(axis='x')
    plt.ylabel('Lifetime')
    plt.xlabel('Time')
    plt.title("Number of data points: %s\nInitial: %s\nOptimum: %s\nLog marginal likelihood: %s\nMSE: %s" % 
        (X.size,kernel, gp.kernel_,gp.log_marginal_likelihood(gp.kernel_.theta),mean_squared_error(y_normal, y_mean)))
    plt.legend()
    plt.tight_layout()
    plt.show()


# + id="Gxfv1lUA6LAl"
### finding best GP model (already done)
# Scoring metrics
# def weighted_mean_sq_error(estimator, x, y):
#     y_pred = estimator.predict(x)
#     # more weight for points with x<1000
#     weight = np.where(x[:,0]<1000,0.75,0.25)
#     weighted_sq = weight*(np.square(y - y_pred))
#     return -weighted_sq.mean()

# def error_lowest_pt(estimator, x, y):
#     # error for the smallest point
#     y_pred = estimator.predict(x)
#     position = np.argmin(y)
#     return -abs(y[position]-y_pred[position])

# # Defining search space
# alpha_space = np.logspace(-8, -5, 7)
# length_scale_space = np.arange(3,4,0.01)
# length_scale_space = np.append(length_scale_space,[2.35,0.27,4.24,6.17,17,0.32]) # starting from 0.32 num_data = 100
# alpha_space = np.append(alpha_space,3e-7)

# param_space = {'alpha': np.logspace(-8, -5, 7),
#                    'kernel': [1.0 * RBF(length_scale = i,length_scale_bounds=[0,1e5]) 
#                                 + WhiteKernel(noise_level=1e-16,noise_level_bounds=(1e-15,1e-11)) for i in length_scale_space]
#                                 + [1.0 * Matern(length_scale = i,length_scale_bounds=[0,1e5]) 
#                                 + WhiteKernel(noise_level=1e-16,noise_level_bounds=(1e-15,1e-11)) for i in length_scale_space]}

# # Performing grid search
# gp = GaussianProcessRegressor()
# gp_search = GridSearchCV(gp, param_space,
#                         scoring={ "weighted_mean_squared_error": weighted_mean_sq_error,
#                                     "log_marginal_likelihood": lambda estimator,
#                                     x,y: estimator.log_marginal_likelihood_value_,
#                                     "error_lowest_point": error_lowest_pt},            
#                         refit = "log_marginal_likelihood"
#                         )
# gp_search.fit(X,Y)
# print(gp_search)
# # ranking models by three different metrics
# models = pd.DataFrame(gp_search.cv_results_)[['param_kernel',
#                             'param_alpha', 'rank_test_weighted_mean_squared_error',
#                             'rank_test_log_marginal_likelihood',
#                             'rank_test_error_lowest_point']].sort_values(by="rank_test_weighted_mean_squared_error")
# models["rank"] = models["rank_test_weighted_mean_squared_error"] + models["rank_test_log_marginal_likelihood"] + models["rank_test_error_lowest_point"]
# models.sort_values(by="rank", inplace=True)
# print(models.iloc[1,:])
# models

# + colab={"base_uri": "https://localhost:8080/"} id="79FF7cW--trM" outputId="789cfd25-4628-4ec8-9e12-acd4ea3b6828"
# The best model from 
gp = GaussianProcessRegressor(alpha=0.000003 ,kernel = 1**2 * RBF(length_scale=3.42) + WhiteKernel(noise_level=1e-16))
#2.35,0.27,4.24,6.17,17,0.32 (length scale parameters that can be considered)
gp.fit(X,Y)
print(gp)
# -

plot(gp,X,temp_list,NUM_OF_PLOT_POINTS)


# # Using LCB

# +
def lcb(X_exp,Y_exp,bounds):
    opt = GPyOpt.methods.BayesianOptimization(f=None,
                                                                   domain=bounds,
                                                                   acquisition_type="LCB",
                                                                   normalize_Y=True,
                                                                   X=X_exp,
                                                                   Y=Y_exp                                                               )
    return opt.suggest_next_locations(), opt

def denormalize(vec, mu, std):
    return vec*std+mu


# -

x_normal = (x_normal - x_normal.min()) / (x_normal.max() - x_normal.min())


# +
bounds = [ {"name":"time", "type": "discrete", "domain":x_, "dimensionality":1},{"name":"temp","type":"discrete","domain":[
                120,140,160,180],"dimensionality":1}]
rnd = rng.integers(0,x_normal.shape[0])
start = x_normal[rnd]
# X_exp = x_normal[start].reshape(1,-1)
X_exp = start.reshape(1,-1)
Y_exp = y_experiment[rnd].reshape(1,-1)
lookup = {tuple(x_normal[i]):y_experiment[i] for i in range(x_normal.shape[0])}

for _ in range(79):
    nxt_pt, opt = lcb(X_exp, Y_exp, bounds)
    X_exp = np.vstack((X_exp,nxt_pt))
    Y_exp = np.vstack((Y_exp, lookup[tuple(nxt_pt[0])]))

# +
global_mu = np.mean(y_experiment)
global_sigma = np.std(y_experiment)
X_120, Y_120 = X_exp[X_exp[:,1]==120], Y_exp[X_exp[:,1]==120]
X_140, Y_140 = X_exp[X_exp[:,1]==140], Y_exp[X_exp[:,1]==140]
X_160, Y_160 = X_exp[X_exp[:,1]==160], Y_exp[X_exp[:,1]==160]
X_180, Y_180 = X_exp[X_exp[:,1]==180], Y_exp[X_exp[:,1]==180]

y_pred, sigma = opt.model.predict(x_normal[:5000,:])
y_pred = denormalize(y_pred,global_mu,global_sigma)
sigma = sigma*global_sigma
# Plot the function, the prediction and the 95% confidence interval based on
# the MSE
plt.figure(figsize=(20,10))
plt.plot(x, y_experiment[:5000], "r:", label="True data")
plt.plot(X_120[:,0], Y_120, "o", markersize=10, label="Observations")
plt.plot(x, y_pred, "b-", label="Prediction")
plt.fill(
    np.concatenate([x, x[::-1]]),
    np.concatenate([y_pred - 1.9600 * sigma, (y_pred + 1.9600 * sigma)[::-1]]),
    alpha=0.5,
    fc="b",
    ec="None",
    label="95% confidence interval",
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.ylim(0.68,0.80)
plt.legend(loc="upper left")
plt.title("Temp:120");

# +
y_pred, sigma = opt.model.predict(x_normal[5000:10000,:])
y_pred = denormalize(y_pred,global_mu,global_sigma)
sigma = sigma*global_sigma
# Plot the function, the prediction and the 95% confidence interval based on
# the MSE
x = np.arange(5000)
plt.figure(figsize=(20,10))
plt.plot(x, y_experiment[5000:10000], "r:", label="True data")
plt.plot(X_140[:,0], Y_140, "g.", markersize=10, label="Observations")
plt.plot(x, y_pred, "b-", label="Prediction")
plt.fill(
    np.concatenate([x, x[::-1]]),
    np.concatenate([y_pred - 1.9600 * sigma, (y_pred + 1.9600 * sigma)[::-1]]),
    alpha=0.5,
    fc="b",
    ec="None",
    label="95% confidence interval",
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.legend(loc="upper left")
plt.ylim(0.700,0.825)

plt.title("Temp:140");

# +
y_pred, sigma = opt.model.predict(x_normal[10000:15000,:])
y_pred = denormalize(y_pred,global_mu,global_sigma)
sigma = sigma*global_sigma
# Plot the function, the prediction and the 95% confidence interval based on
# the MSE
x = np.arange(5000)
plt.figure(figsize=(20,10))
plt.plot(x, y_experiment[10000:15000], "r:", label="True data")
plt.plot(X_160[:,0], Y_160, "g.", markersize=10, label="Observations")
plt.plot(x, y_pred, "b-", label="Prediction")
plt.fill(
    np.concatenate([x, x[::-1]]),
    np.concatenate([y_pred - 1.9600 * sigma, (y_pred + 1.9600 * sigma)[::-1]]),
    alpha=0.5,
    fc="b",
    ec="None",
    label="95% confidence interval",
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.legend(loc="upper left")
plt.ylim(0.68,0.825)

plt.title("Temp:160");
# -

y_pred, sigma = opt.model.predict(x_normal[15000:,:])
y_pred = denormalize(y_pred,global_mu,global_sigma)
sigma = sigma*global_sigma
# Plot the function, the prediction and the 95% confidence interval based on
# the MSE
x = np.arange(5000)
plt.figure(figsize=(20,10))
plt.plot(x, y_experiment[15000:], "r:", label="True data")
plt.plot(X_180[:,0], Y_180, "g.", markersize=10, label="Observations")
plt.plot(x, y_pred, "b-", label="Prediction")
plt.fill(
    np.concatenate([x, x[::-1]]),
    np.concatenate([y_pred - 1.9600 * sigma, (y_pred + 1.9600 * sigma)[::-1]]),
    alpha=0.5,
    fc="b",
    ec="None",
    label="95% confidence interval",
)
plt.xlabel("$x$")
plt.ylabel("$f(x)$")
plt.legend(loc="upper left")
plt.ylim(0.68,0.825)
plt.title("Temp:180");

y_pred

#

StarcoderdataJupyterScriptsDedupFiltered