# ---
# 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
# ---

# # Introduction to Gaussian Process Models
# Gaussian process (GP) models serve as approximations of computationally expensive (time-consuming) black-box functions. To reduce the number of times the expensive function must be queried during optimization, the GP is used to guide the sampling decisions in the parameter space and only the "most promising" parameters are selected for evaluation. 

# A GP model treats the function it approximates like the realization of a stochastic process: 
# $m_{GP}(\theta) = \mu + Z(\theta)$,
# where $\mu$ represents the mean of the stochastic process and $Z(\theta) \sim \mathcal{N}(0,\sigma^2)$ is the deviation from the mean. 

# The correlation between two random variables $Z(\theta_k)$ and $Z(\theta_l)$ is defined by a kernel, e.g., the squared exponential (also Radial basis function) kernel: 
# \begin{equation}
# Corr(Z(\theta_k),Z(\theta_l)) = \exp(-\sum_{i=1}^d \gamma_i|\theta_k^{(i)}-\theta_l^{(i)}|^{q_i})
# \end{equation}, 
# with $\gamma_i$ determining how quickly the correlation in dimension $i$ decreases, and $q_i$ refelcts the smoothness of the function in dimension $i$

# Denoting $\mathbf{R}$ as the matrix whose $(k,l)$-th element is given as the correlation above, maximum likelihood estimation is used to determine the GP parameters $\mu$, $\sigma^2$, and $\gamma_i$. Then, at an unsampled point $\theta^{new}$, the GP prediction is \begin{equation}
#     m_{\text{GP}}(\theta^{\text{new}})=\hat{\mu}+\mathbf{r}^T\mathbf{R}^{-1}(\mathbf{f}-\mathbf{1}\hat\mu),
# \end{equation}
# where $\mathbf{1}$ is a vector of ones of appropriate dimension and $\mathbf{f}$ is the vector of function values obtained so far, and
# \begin{equation}
#     \boldsymbol{r}=
#     \begin{bmatrix}
#     Corr\left(Z(\theta^{\text{new}}), Z(\theta_1)\right)\\
#     \vdots\\
#     Corr\left(Z(\theta^{\text{new}}
#     ), Z(\theta_n)\right)
#     \end{bmatrix}.
# \end{equation}
# The  corresponding  mean squared error is 
# \begin{equation}
#     s^2(\theta^{\text{new}})=\hat{\sigma}^2\left(   1-\boldsymbol{r}^T\boldsymbol{R}^{-1}\boldsymbol{r} +\frac{(1-\boldsymbol{1}^T\boldsymbol{R}^{-1}\boldsymbol{r})^2}{\mathbf{1}^T\boldsymbol{R}^{-1}\mathbf{1}}\right)
# \end{equation}
# with 
#  \begin{equation}
#      \hat{\mu} = \frac{\mathbf{1}^T\boldsymbol{R}^{-1}\mathbf{f}}{\mathbf{1}^T\boldsymbol{R}^{-1}\mathbf{1}}
#  \end{equation}
#  and
#  \begin{equation}
#      \hat{\sigma}^2=\frac{(\mathbf{f}-\mathbf{1}\hat{\mu})^T\boldsymbol{R}^{-1}(\mathbf{f}-\mathbf{1}\hat{\mu})}{n}.
#  \end{equation}

# Python has a good implementation of GPs where you can choose different kernels. 

# First, we need (input, output) data pairs. Inputs are parameters where we query the function (for simplicity, the example has an inexpensive function). From the Sckit-Learn website: https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_noisy_targets.html

import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, Matern, RationalQuadratic, ExpSineSquared, WhiteKernel
from scipy.optimize import minimize
from scipy.spatial import distance
import scipy.spatial as scp
from scipy.stats import norm
from pyDOE import * #needed if Latin hypercube design is used
import warnings
warnings.filterwarnings("ignore")


def f(x):
    """The function we want to approximate."""
    return x * np.sin(x)


xlow = 0 #lower bound on x
xup = 10 #upper bound on x
dim = 1 #dimension of the problem
lhs_wanted = False
np.random.seed(420)

if not(lhs_wanted): #when not using space-filling design
    X = np.atleast_2d([1., 3., 7., 8.]).T #select some points where we evaluate the function

    # Function evaluations
    y = f(X).ravel()

# Other options for creating space filling designs is latin hypercube sampling:

if lhs_wanted:
    ninit=6 #6 initial evaluations
    init_design = lhs(dim, samples =ninit, criterion='maximin') #initial design in [0,1]^dim
    X = xlow+(xup-xlow)*init_design #scale to [xlow,xup]
    # Function evaluations
    y = f(X).ravel()

# **Exercise:** run the code with different initial samples, i.e., try lhs_wanted = False and lhs_wanted = True and compare the sampling history

# +
# Select a GP kernel (here RBF or squared exponential)
kernel = RBF()
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True,n_restarts_optimizer=9)

# Fit the GP to the input-output data 
gp.fit(X, y)

# -

# Make some good-looking plots

def plot_the_gp(X, y, gp, xnew):
    #select a bunch of points where we want to make predictions wioth the GP
    x = np.atleast_2d(np.linspace(0, 10, 1000)).T 
    # Make the GP prediction at the points where no evaluations were taken - also return predicted uncertainty
    y_pred, sigma = gp.predict(x, return_std=True)
    plt.figure()
    plt.plot(x, f(x), 'r:', label=r'$f(x) = x\,\sin(x)$')
    plt.plot(X, y, 'r.', markersize=10, label='Observations')
    plt.plot(x, y_pred, 'b-', label='Prediction')
    if len(xnew)>0:
        plt.plot(X[-1], y[-1], 'gs', markersize=10, label='Newest sample')
    plt.fill(np.concatenate([x, x[::-1]]),
         np.concatenate([y_pred - 1.9600 * sigma,
                        (y_pred + 1.9600 * sigma)[::-1]]),
         alpha=.5, fc='b', ec='None', label='95% confidence interval')
    plt.xlabel('$x$')
    plt.ylabel('$f(x)$')
    plt.ylim(-10, 20)
    plt.legend(loc='upper left')
    plt.show()


plot_the_gp(X, y, gp, [])


# **Optional Exercise:** check out the Scikit-Learn website https://scikit-learn.org/stable/modules/gaussian_process.html#kernels-for-gaussian-processes and experiment around with different basic kernels, kernel parameters and kernel combinations, e.g., 
# - does using "kernel = RBF(10, (1e-2, 1e2))" change anything?
# - what happens when you use "kernel = Matern(length_scale=1.0, nu=1.5)"
# - try "kernel = 1.0 * RationalQuadratic(length_scale=1.0, alpha=0.1, alpha_bounds=(1e-5, 1e15))"
# - "kernel = 1.0 * ExpSineSquared(
#     length_scale=1.0,
#     periodicity=3.0,
#     length_scale_bounds=(0.1, 10.0),
#     periodicity_bounds=(1.0, 10.0),)"
# - use a combination of kernels: "kernel = RBF()+WhiteKernel(noise_level=.001)" using different noise_levels
#

# **Exercise:** Change the inputs of the GP (i.e., the training samples) and see how the GP predictions change (use fewer or more points, use different points in [0,10], e.g., "X=np.atleast_2d(np.random.uniform(0,10,5)).T"

# Takeaway: the quality and accuracy of the GP highly depends on the trianing data and the kernel used

# # Adaptive Optimization with the GP

# GP models are often used in optimization algorithms. In each iteration of the optimization, a new sample point is selected by maximizing the expected improvement (EI):
# \begin{equation}
#     \mathbb{E}(I)(\theta) = s(\theta)\left(v\Phi(v)+\phi(v)   \right),
# \end{equation}
# where 
# \begin{equation}
#     v=\frac{f^{\text{best}}-m_{\text{GP}}(\theta)}{s(\theta)}
# \end{equation}
# and $\Phi$ and $\phi$ are the normal cumulative distribution and density functions, respectively, and $s(\theta)=\sqrt{s^2(\theta)}$ is the square root of the mean squared error. 
#

# The function $\mathbb{E}(I)(\theta)$ can be maximized with any python optimization library. The point $\theta^{\text{new}}$ where it reaches its maximum will be the new point where $f$ is evaluated.

# define expected improvement function
def ei(x, gpr_obj, Xsamples, Ysamples): #expected improvement
    dim = len(x)
    x= x.reshape(1, -1)

    min_dist=np.min(scp.distance.cdist(x, Xsamples))
    if min_dist<1e-6: #threshold for when points are so close that they are considered indistinguishable
        expected_improvement=0.0
        return expected_improvement

    mu, sigma = gpr_obj.predict(x.reshape(1, -1), return_std=True)
    mu_sample = gpr_obj.predict(Xsamples)
    mu_sample_opt = np.min(Ysamples)

    # In case sigma equals zero
    with np.errstate(divide='ignore'):
        Z = (mu_sample_opt-mu) / sigma
        expected_improvement = (mu_sample_opt-mu) * norm.cdf(Z) + sigma * norm.pdf(Z)
        expected_improvement[sigma == 0.0] == 0.0
    answer=-1.*expected_improvement #to maximize EI, you minimize the negative of it
    return answer


def plot_the_ei(gpr_obj, X, Y):
    x = np.atleast_2d(np.linspace(0, 10, 1000)).T 
    expimp=np.zeros(1000)
    for ii in range(1000):
        expimp[ii] = -ei(x[ii], gpr_obj, X, Y)
    plt.figure()
    plt.plot(x, expimp, 'k--', label='Expected improvement')
    plt.plot(X, np.zeros(X.shape[0]), 'rx', markersize=10, label='Observation sites')
    #plt.plot(X[-1],0, 'gs', markersize=10, label='Newest sample')
    plt.xlabel('$x$')
    plt.ylabel('$EI(x)$')
    #plt.ylim(-10, 20)
    plt.legend(loc='upper left')
    plt.show()


# do your GP iterations: maximize EI, select new point, evaluate new point, update GP, maximize EI, ....
n_GP_samples = 20 # allow 50 evaluations of f
bound_list = np.array([[xlow, xup]])
xnew=[]
while  X.shape[0]< n_GP_samples: 
    gpr_obj = GaussianProcessRegressor(kernel=kernel, random_state=0,normalize_y=True, n_restarts_optimizer=10).fit(X, y) #create the GP
    plot_the_gp(X, y, gpr_obj, xnew)
    plot_the_ei(gpr_obj, X, y)
    #compute next point by maximizing expected improvement, multi-start optimization
    xnew = []
    fnew =np.inf
    for ii in range(10):
        x0 = xlow + (xup-xlow) * np.random.rand(1,dim) #random starting point for optimizing expected improvement
        res= minimize(ei,np.ravel(x0),method='SLSQP',bounds=bound_list, args=(gpr_obj, X, y))	
        dist = np.min(scp.distance.cdist(np.asmatrix(res.x), X)) #make sure new point is sufficiently far away from already sampled points
        if np.min(dist)>1e-6 and res.success: #1e-3 is tunable
            x_ = np.asmatrix(res.x)
            if res.fun< fnew:
                xnew = x_
                fnew = res.fun
        else: #use random point as new point
            x_ = np.asarray(xlow) + np.asarray(xup-xlow) * np.asarray(np.random.rand(1,dim)) #random starting point
            fv= ei(x_, gpr_obj, X, y)
            if len(xnew)== 0 or fv < fnew:
                xnew = np.asmatrix(x_)
                fnew= fv

                
    fval = f(np.ravel(xnew))

    #update Xsamples and Ysamples arrays
    X=np.concatenate((X, np.asmatrix(xnew)), axis = 0)
    Y_ = np.zeros(len(y)+1)    
    Y_[0:len(y)]= y
    Y_[-1]=fval
    y =Y_
    minID=np.argmin(y) #find index of best point
    print('best point: ', X[minID])
    print('best value: ', y[minID])
    print('Number evaluations: ', X.shape[0])

# From the images of the expected improvement we can see that the peaks are becoming increasingly narrow, to almost looking like jump-discontinuities. This means that for an optimizer that tries to find the maximum of the expected improvement function, it becomes increasingly harder to find the optimum and sampling becomes more "random" in the space because the function is flat and EI values are the same everywhere except at the jumps.

# Takeaways: 
# - GPs can be useful to guide the search during optimization
# - They shine when the number of function evaluations is severely limited
# - The expected improvement function helps to select points that are the "most promising" next evaluations
# - The expected improvement function is multimodal and becomes increasingly harder to optimize

StarcoderdataJupyterScriptsDedupFiltered