ActiveBOToytask/BayesianOptimization/BOtest_for_Preference.py

import numpy as np
from scipy.optimize import minimize
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from AcquistionFunctions.ExpectedImprovementwithPreference import expected_improvement

from sklearn.exceptions import ConvergenceWarning
import warnings

warnings.filterwarnings("ignore", category=ConvergenceWarning)



# Define the Branin function
def branin(X):
    x, y = X
    result = (y - (5.1 / (4*np.pi**2)) * x**2 +
              5 * x / np.pi - 6)**2 + 10 * (1 - 1 / (8*np.pi)) * np.cos(x) + 10
    return result


# Define the Rosenbrock function
def rosenbrock(X):
    x, y = X
    return (1 - x)**2 + 100 * (y - x**2)**2


# Function to plot the Test function and the samples
def plot_function_and_samples(func, X_samples, Y_samples, iteration):
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    # Create a grid of points
    x = np.linspace(-2, 2, 100)
    y = np.linspace(-1, 3, 100)
    X, Y = np.meshgrid(x, y)
    Z = func([X, Y])

    # Plot the function
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, color='b', alpha=0.2, linewidth=0)

    # Plot the samples
    ax.scatter(X_samples[:, 0], X_samples[:, 1], Y_samples, color='r', s=50)

    ax.set_title(f"Iteration {iteration}")
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('f(X, Y)')
    ax.view_init(elev=30, azim=-30)

    plt.show()

def bayesian_optimization(n_iter, change_to_gaussian_after=10):
    # # The bounds of the function
    # bounds = np.array([[-5.0, 10.0], [0.0, 15.0]]) #brainin
    bounds = np.array([[-2.0, 2.0], [-1.0, 3.0]]) #rosenbrock

    # Initial samples
    X_init = np.random.uniform([b[0] for b in bounds], [b[1] for b in bounds], size=(2, len(bounds)))
    Y_init = np.array([[rosenbrock(x)] for x in X_init])

    # Gaussian process with Matérn kernel as surrogate model
    m52 = Matern(length_scale=[1.0, 1.0], nu=2.5)
    gp = GaussianProcessRegressor(kernel=m52, alpha=1e-5)

    # Boolean array indicating which dimensions should be sampled from a Gaussian distribution
    use_gaussian = np.array([False, False])

    # Means and variances for the Gaussian distributions
    means = np.array([0.0, 0.0])
    variances = np.array([(b[1] - b[0]) * 0.25 for b in bounds])**2

    # Array to store the best result from each iteration
    best_results = np.empty(n_iter)

    # Bayesian optimization loop
    for i_ in range(n_iter):
        # Update Gaussian process with existing samples
        gp.fit(X_init, Y_init)

        # Obtain next sampling point from the acquisition function (expected_improvement)
        n_samples = 1000
        X_next = expected_improvement(X_init, n_samples, gp, use_gaussian, means, variances)

        # Obtain next noisy sample from the objective function
        Y_next = rosenbrock(X_next)

        # Add sample to previous samples
        X_init = np.vstack((X_init, X_next))
        Y_init = np.vstack((Y_init, Y_next))

        # After change_to_gaussian_after iterations, start using Gaussian distribution for sampling
        if i_ >= change_to_gaussian_after:
            use_gaussian = np.array([True, True])
            # Update means with the mean of the best 5 results
            best_5_idx = np.argsort(Y_init, axis=0)[:5].flatten()
            means = np.mean(X_init[best_5_idx], axis=0)
            # Reduce variances by 5 percent
            variances *= 0.95

        # Store the best result from this iteration
        best_results[i_] = np.min(Y_init)

    return best_results


# Run the benchmarking test and store the results
results_with_gaussian = np.empty((50, 20))
results_without_gaussian = np.empty((50, 20))

for i in range(50):
    results_with_gaussian[i, :] = bayesian_optimization(20, 10)
    results_without_gaussian[i, :] = bayesian_optimization(20, 50)

# Compute the mean and standard deviation of the results
mean_with_gaussian = np.mean(results_with_gaussian, axis=0)
std_with_gaussian = np.std(results_with_gaussian, axis=0)
mean_without_gaussian = np.mean(results_without_gaussian, axis=0)
std_without_gaussian = np.std(results_without_gaussian, axis=0)

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(mean_with_gaussian, label='With Gaussian after 10 iterations')
plt.fill_between(range(20), mean_with_gaussian-std_with_gaussian, mean_with_gaussian+std_with_gaussian, alpha=0.2)
plt.plot(mean_without_gaussian, label='Without Gaussian')
plt.fill_between(range(20), mean_without_gaussian-std_without_gaussian, mean_without_gaussian+std_without_gaussian, alpha=0.2)
plt.xlabel('Iteration')
plt.ylabel('Best Result')
plt.legend()
plt.show()


#
# # Initial samples
# X_init = np.random.uniform([b[0] for b in bounds], [b[1] for b in bounds], size=(2, len(bounds)))
# # Y_init = np.array([[branin(x)] for x in X_init])
# Y_init = np.array([[rosenbrock(x)] for x in X_init])
#
# # Gaussian process with Matérn kernel as surrogate model
# m52 = Matern(length_scale=[1.0, 1.0], nu=2.5)
# gp = GaussianProcessRegressor(kernel=m52, alpha=1e-5)
#
# # Number of iterations
# n_iter = 50
#
# # Boolean array indicating which dimensions should be sampled from a Gaussian distribution
# use_gaussian = np.array([False, False])
#
# # Means and variances for the Gaussian distributions
# means = np.array([0.0, 0.0])
# variances = np.array([(b[1] - b[0]) * 0.25 for b in bounds])**2
#
# # Bayesian optimization loop
# for i in range(n_iter):
#     # Update Gaussian process with existing samples
#     gp.fit(X_init, Y_init)
#
#     # Obtain next sampling point from the acquisition function (expected_improvement)
#     n_samples = 1000
#     X_next = expected_improvement(X_init, n_samples, gp, use_gaussian, means, variances)
#
#     # Obtain next noisy sample from the objective function
#     # Y_next = branin(X_next)
#     Y_next = rosenbrock(X_next)
#
#     # Add sample to previous samples
#     X_init = np.vstack((X_init, X_next))
#     Y_init = np.vstack((Y_init, Y_next))
#
#
#     # After 10 iterations, start using Gaussian distribution for sampling
#     if i >= 29:
#         use_gaussian = np.array([True, True])
#         # Update means with the mean of the best 5 results
#         best_5_idx = np.argsort(Y_init, axis=0)[:5].flatten()
#         means = np.mean(X_init[best_5_idx], axis=0)
#         # Reduce variances by 10 percent
#         variances *= 0.9
#
#     # Plot the function and the samples
#     # plot_function_and_samples(rosenbrock, X_init, Y_init, i)
#
# print("Final samples: ", X_init)
# print("Final function values: ", Y_init)
# print(f"Minimum function values: {np.min(Y_init)}, idx: {np.argmin(Y_init)}")