<jupyter_start><jupyter_text>Define the function<jupyter_code>oneD_function =lambda x: cos(1000*x-500)+(abs(x*100_000))/(x**4+1000)
plt.plot(np.linspace(-10,10),[oneD_function(_) for _ in np.linspace(-10,10)])
plt.title('1D Demo Function');<jupyter_output><empty_output><jupyter_text>Setup the package BayesianOptimization<jupyter_code>bo = BayesianOptimization(oneD_function,
                          {'x': (0, 50)})
bo.explore({'x': np.linspace(0,50)})<jupyter_output><empty_output><jupyter_text>Setup my implementation<jupyter_code>test_x = np.array(np.linspace(0,50,1000)).reshape(-1,1)
train_x = np.array(np.linspace(0,50,1)).reshape(-1,1)
train_y_numbers = np.array([oneD_function(_) for _ in train_x]).reshape(-1,1)
y_var = 10<jupyter_output><empty_output><jupyter_text>My implementation to run Bayesian Optimization First, sample points and update the posterior<jupyter_code>def test_gaussian_process(test_x, train_x, train_y_numbers, y_var, kernel, return_cov = False,
                          return_sample = False, covariance_noise = 5e-5, n_posteriors = 5):
    """Test one new point in the Gaussian process and update the posterior or an array of points
    Returns the mean, var from normal distribution from the sampled point.
    Return cov = True will return the full covariance matrix.
    
    If return_sample= True 
    returns samples ( a vector) from the informed posterior and the uninformed prior distribution
    
    Covariance diagonal noise is used to help enforce positive definite matrices
    
    N_posteriors indicates the number of posterior functions to create when sampling"""
    
    from numpy import linalg as la
    #https://stackoverflow.com/questions/43238173/python-convert-matrix-to-positive-semi-definite/43244194#43244194
    def nearestPD(A):
        """Find the nearest positive-definite matrix to input

        A Python/Numpy port of <NAME>'s `nearestSPD` MATLAB code [1], which
        credits [2].

        [1] https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd

        [2] <NAME>, "Computing a nearest symmetric positive semidefinite
        matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6
        """

        B = (A + A.T) / 2
        _, s, V = la.svd(B)

        H = np.dot(V.T, np.dot(np.diag(s), V))

        A2 = (B + H) / 2

        A3 = (A2 + A2.T) / 2

        if isPD(A3):
            return A3

        spacing = np.spacing(la.norm(A))
        # The above is different from [1]. It appears that MATLAB's `chol` Cholesky
        # decomposition will accept matrixes with exactly 0-eigenvalue, whereas
        # Numpy's will not. So where [1] uses `eps(mineig)` (where `eps` is Matlab
        # for `np.spacing`), we use the above definition. CAVEAT: our `spacing`
        # will be much larger than [1]'s `eps(mineig)`, since `mineig` is usually on
        # the order of 1e-16, and `eps(1e-16)` is on the order of 1e-34, whereas
        # `spacing` will, for Gaussian random matrixes of small dimension, be on
        # othe order of 1e-16. In practice, both ways converge, as the unit test
        # below suggests.
        I = np.eye(A.shape[0])
        k = 1
        while not isPD(A3):
            mineig = np.min(np.real(la.eigvals(A3)))
            A3 += I * (-mineig * k**2 + spacing)
            k += 1

        return A3

    def isPD(B):
        """Returns true when input is positive-definite, via Cholesky"""
        try:
            _ = la.cholesky(B)
            return True
        except la.LinAlgError:
            return False
    
    
    
    # define the covariance matrices
    covariance_train_train = kernel(train_x,train_x,train=True)
    covariance_test_train  = kernel(test_x,train_x)
    covariance_test_test  = kernel(test_x,test_x)
    

    # Use cholskey decomposition to increase speed for calculating mean
    try :# First try, 
        L_test_test = np.linalg.cholesky(covariance_test_test + covariance_noise*np.eye(len(covariance_test_test)))
        L_train_train = np.linalg.cholesky(covariance_train_train + covariance_noise*np.eye(len(covariance_train_train)))
        Lk = np.linalg.solve(L_train_train, covariance_test_train.T)
        mus = np.dot(Lk.T, np.linalg.solve(L_train_train, train_y_numbers)).reshape((len(test_x),))
        # Compute the standard deviation so we can plot it
        s2 = np.diag(covariance_test_test) - np.sum(Lk**2, axis=0)
        stdv = np.sqrt(abs(s2))
        

    # Full matrix calculation of mean and covariance, much slower than cholesky decomposition
    except Exception as e:
        print(e)#LinAlgError: # In case the covariance matrix is not positive definite 
        # Find the near positive definite matrix to decompose
        decompose_train_train = nearestPD(covariance_train_train + covariance_noise*np.eye(len(train_x)))
        decompose_test_test = nearestPD(covariance_test_test + covariance_noise*np.eye(len(test_x)))
        
        # cholskey decomposition within the try except block
        L_train_train = np.linalg.cholesky(decompose_train_train  )
        L_test_test = np.linalg.cholesky(decompose_test_test  )
        Lk = np.linalg.solve(L_train_train, covariance_test_train.T)
        mus = np.dot(Lk.T, np.linalg.solve(L_train_train, train_y_numbers)).reshape((len(test_x),))
        # Compute the standard deviation so we can plot it
        s2 = np.diag(covariance_test_test) - np.sum(Lk**2, axis=0)
        stdv = np.sqrt(abs(s2))
        
        
#         ##### FULL INVERSION ####
#         mus = covariance_test_train  @ np.linalg.pinv(covariance_train_train) @ train_y_numbers
#         s2 = covariance_test_test - covariance_test_train @ np.linalg.pinv(covariance_train_train ) \
#                      @ covariance_test_train.T
        

    def sample_from_posterior(n_priors=3, n_post=5):
        """Draw samples from the prior distribution of the GP.
        len(test_x) is the number of samplese to draw.
        Resource: http://katbailey.github.io/post/gaussian-processes-for-dummies/.
        
        N-Posteriors / N-Priors tells the number of functions to samples from the dsitribution"""

        # Draw samples from the posterior at our test points.
        covariance_test_test  = kernel(test_x,test_x)
        try: # try inside sample from posterior function
            L = np.linalg.cholesky(covariance_test_test +  covariance_noise*np.eye(len(test_x))- np.dot(Lk.T, Lk))
        except Exception as e:
            # Find the neareset Positive Definite Matrix
            near_decompose = nearestPD(covariance_test_test +  covariance_noise*np.eye(len(test_x))- np.dot(Lk.T, Lk))
            L = np.linalg.cholesky(near_decompose.astype(float) ) # within posterior
        # sample from the posterior
        f_post = mus.reshape(-1,1) + np.dot(L, np.random.normal(size=(len(test_x),n_posteriors)))

        # Sample X sets of standard normals for our test points,
        # multiply them by the square root of the covariance matrix
        f_prior_uninformed = np.dot(L_test_test, np.random.normal(size=(len(test_x),n_priors)))
        # For the posterior, the columns are the vector for that function
        return (f_prior_uninformed, f_post)


    if return_cov == True:
        return y_pred_mean.ravel(), var_y_pred_diag.ravel(), var_y_pred
        
    if return_sample == True:
        f_prior, f_post = sample_from_posterior(n_post = n_posteriors)
        return mus.ravel(), s2.ravel(), f_prior, f_post
    else:
        return mus.ravel(), s2.ravel()<jupyter_output><empty_output><jupyter_text>Second, define the squarred kernel<jupyter_code>def squarred_kernel(a, b, param=2.0, train=False, train_noise = 5e-3, vertical_scale=50):
    """Calculated the squarred exponential kernel. 
    Adds a noise term for the covariance of the training data
    Adjusting the param changes the difference where points will have a positive covariance
    Returns a covaraince Matrix.
    Vertical scale controls the vertical scale of the function"""
    if train == False:
        a = np.array(a)
        b = np.array(b)
        
        sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a, b.T)
        return vertical_scale*np.exp(-.5 * (1/param) * sqdist)

    else:
        noisy_observations = train_noise*np.eye(len(a))
        a = np.array(a)
        b = np.array(b)
        
        sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a, b.T) 
        return vertical_scale*np.exp(-.5 * (1/param) * sqdist) + noisy_observations<jupyter_output><empty_output><jupyter_text>Next, define the acquisition function<jupyter_code>def expected_improvement(mean_x, sigma_squarred_x, y_val_for_best_hyperparameters, normal_dist=None, point_est = False):
    """Finds the expected improvement of a point give the current best point.
    If point_est = False, then computes the expected value on a vector from the posterior distribution.
    """
    from scipy.stats import norm
    
    with np.errstate(divide='ignore'): # in case sigma equals zero
        # Expected val for one point
        if point_est ==True:
            sigma_x = np.sqrt(sigma_squarred_x) # get the standard deviation from the variance

            Z = (mean_x - y_val_for_best_hyperparameters) / sigma_x

            if round(sigma_x,8) == 0:
                return 0
            else:
                return (mean_x - y_val_for_best_hyperparameters)*normal_dist.cdf(Z) + sigma_x*normal_dist.pdf(Z)

        else:
            # Sample from the posterior functions
            for _ in range(len(mean_x)):
                list_of_improvements = []
                m_s = []
                for m,z,s in zip(mean_x, ((mean_x -y_val_for_best_hyperparameters)\
                                     / np.std(sigma_squarred_x)),np.sqrt(sigma_squarred_x) ):
                    
                    list_of_improvements.append(((m-y_val_for_best_hyperparameters)*norm().cdf(z)\
                                                 +s*norm().pdf(z)))
                    m_s.append(m)

                return list_of_improvements<jupyter_output><empty_output><jupyter_text>Then, create the function to predict the next points to sample<jupyter_code>def sample_from_function( y_function, 
                        train_x, train_y_numbers, kernel, y_var, testing_domain = range(-100,100,10), n = 20,
                        verbose=True):
    """Sample N times from the unknown function and for each time find the 
    point that will have the highest expected improvement (find the maxima of the function).
    Verbose signifies if the function should print out the points where it is sampling
    
    Returns a generator of x_sampled_points, y_sampled_points, vars_, best_x, best_y, \
                list_of_expected_improvements, testing_domain, mus
          for improvements. Mus and Vars are the mean and var for each sampled point in the gaussian process. 
    
    Starts off the search for expected improvement with a coarse search and then hones in on 
    the domain the the highest expected improvement.
    
    Note - the y-function can EITHER by the actual y-function (for evaluation purposes, or an objective function 
    (i.e. - RMSE))"""

    # STARTING POINTS
    if train_x.shape[1]>1: # Two dimensional 
        # Random starting point
        try:
            start_point = np.array([testing_domain[np.random.choice(len(testing_domain)),:]])
        except:
            start_point = np.array([testing_domain[np.random.choice(len(testing_domain))]])
        #start at random point in the domain
        best_x = start_point
        try: # in case we are passing the actual function here
            best_y = y_function(start_point[0][0] , start_point[0][1])
        except: # we are passing the objective function (min RMSE)
            best_y = y_function(start_point[0][0],dimensions = 'two',hyperparameter_value_two = start_point[0][1]) 
        # for plotting the points sampled
        x_sampled_points = []
        y_sampled_points = []
        sampled_variance = [] 
        
    else: # One dimensional case
        try:
            start_point = np.random.choice(testing_domain)#start at random point in the domain
        except:
            start_point = np.random.choice(list(testing_domain.ravel()))
        best_x = start_point
        best_y = y_function(start_point) # initial best y value based on start point
        start_y = y_function(start_point) 
        # for plotting the points sampled
        x_sampled_points = []
        y_sampled_points = []
        sampled_variance =[] 
    

    for i in range(n):
        if i == 0:
            if train_x.shape[1]==1: ## one dimensional case
                testing_domain = np.array([testing_domain]).reshape(-1,1)
                        
            # find the next x-point to sample
            mus, vars_,_,post = test_gaussian_process(testing_domain, train_x, train_y_numbers,
                                               y_var, kernel, return_sample = True, n_posteriors=30)
            
           
            sigmas_post = np.var(post,axis=1)
            mus_post = np.mean(post,axis=1)
            # get the expected values from the posterior distribution
            list_of_expected_improvements = expected_improvement(mus_post, sigmas_post ,best_y)

            max_improv_x_idx = np.argmax(list_of_expected_improvements)
            max_improv_x = testing_domain[max_improv_x_idx]
            
            if train_x.shape[1]==1:
                max_improv_y = y_function(max_improv_x)
            else: # Two D
                try: # see if we are passing in the actual function
                    max_improv_y = y_function(max_improv_x[0], max_improv_x[1])
                except: # we are passing the objective function in     
                    max_improv_y = y_function(max_improv_x[0],dimensions = 'two', 
                                              hyperparameter_value_two = max_improv_x[1])
                
            if max_improv_y > best_y: ## use to find out where to search next
                best_y = max_improv_y 
                best_x = max_improv_x
            else:
                best_x = start_point
            
            if verbose:
                print(f"Bayesian Optimization just sampled point = {best_x}")
                print(f"Best x (Bayesian Optimization) = {best_x}, Best y = {best_y}")
                sampled_variance.append(sigmas_post[max_improv_x_idx])

                # append the point to sample
                x_sampled_points.append(max_improv_x)
                y_sampled_points.append(max_improv_y)

                # append our new the newly sampled point to the training data
                train_x = np.vstack((train_x, max_improv_x))
                train_y_numbers = np.vstack((train_y_numbers, max_improv_y))

                yield x_sampled_points, y_sampled_points, vars_, best_x, best_y, \
                    list_of_expected_improvements, testing_domain, mus
                    
            else:

                sampled_variance.append(sigmas_post[max_improv_x_idx])

                # append the point to sample
                x_sampled_points.append(max_improv_x)
                y_sampled_points.append(max_improv_y)

                # append our new the newly sampled point to the training data
                train_x = np.vstack((train_x, max_improv_x))
                train_y_numbers = np.vstack((train_y_numbers, max_improv_y))

                yield x_sampled_points, y_sampled_points, vars_, best_x, best_y, \
                    list_of_expected_improvements, testing_domain, mus
            

        else:
            # reformat testing domain to include more possible points
            if train_x.shape[1]==1:
                testing_domain = np.array([testing_domain]).reshape(-1,1)
            
            
            mus, vars_,_,post = test_gaussian_process(testing_domain, train_x, train_y_numbers,
                                               y_var, kernel, return_sample = True, n_posteriors=30)

            sigmas_post = np.var(post,axis=1)
            mus_post = np.mean(post,axis=1)
            
            # get the expected values from the posterior distribution
            list_of_expected_improvements = expected_improvement(mus_post, sigmas_post, best_y)
            
            
            # find the best x value
            max_improv_x_idx = np.argmax(list_of_expected_improvements)
            max_improv_x = testing_domain[max_improv_x_idx]
            
            if train_x.shape[1]==1:
                max_improv_y = y_function(max_improv_x)
            else: # Two D
                try: # see if we are passing in the actual function
                    max_improv_y = y_function(max_improv_x[0], max_improv_x[1])
                except: # we are passing the objective function in     
                    max_improv_y = y_function(max_improv_x[0],dimensions = 'two', 
                                              hyperparameter_value_two = max_improv_x[1])
                
  
            if max_improv_y > best_y: ## use to find out where to search next
                best_y = max_improv_y 
                best_x = max_improv_x
            else:
                pass
            if verbose:
                print(f"Bayesian Optimization just sampled point = {max_improv_x}")
                print(f"Best x (Bayesian Optimization) = {best_x}, Best y = {best_y}")
                sampled_variance.append(sigmas_post[max_improv_x_idx])

                # append the point to sample
                x_sampled_points.append(max_improv_x)
                y_sampled_points.append(max_improv_y)

                # append our new the newly sampled point to the training data
                train_x = np.vstack((train_x, max_improv_x))
                train_y_numbers = np.vstack((train_y_numbers, max_improv_y))

                yield x_sampled_points, y_sampled_points, vars_, best_x, best_y, \
                    list_of_expected_improvements, testing_domain, mus
                    
            else:

                sampled_variance.append(sigmas_post[max_improv_x_idx])

                # append the point to sample
                x_sampled_points.append(max_improv_x)
                y_sampled_points.append(max_improv_y)

                # append our new the newly sampled point to the training data
                train_x = np.vstack((train_x, max_improv_x))
                train_y_numbers = np.vstack((train_y_numbers, max_improv_y))

                yield x_sampled_points, y_sampled_points, vars_, best_x, best_y, \
                    list_of_expected_improvements, testing_domain, mus<jupyter_output><empty_output><jupyter_text>Compare my implementation to the BayesianOptimization package- One dimension first<jupyter_code>bo.maximize(init_points=1, n_iter=10, acq='ei')
# my implementation
# Create the generator, ONE DIMENSIONAL
samples = 10
gen = sample_from_function(oneD_function, train_x, train_y_numbers,squarred_kernel, y_var, n=samples,
                           testing_domain=np.linspace(-10,10,1000) )

for i in range(samples):
    print('STEP Number = ',i+1)
    # sample from generator
    x_sampled_points_1d, y_sampled_points_1d, sampled_var_1d, best_x_1d, best_y_1d,\
    improvements_1d, domain_1d, mus_1d = next(gen)
# Finally, we take a look at the final results.
print("Bayesian Optimization package results",bo.res['max'])

print(f"The best y-value for the 1D case that my implementation found was {best_y_1d} with x_val = {best_x_1d}")<jupyter_output>The best y-value for the 1D case that my implementation found was [ 317.97323004] with x_val = [ 4.03403403]<jupyter_text>Next, explore a two-dimensional comparison $f(x,y)=100{\sqrt {\left|y-0.01x^{2}\right|}}+0.01\left|x+10\right|.\quad $<jupyter_code>twoD_function = lambda x,y: 100*np.sqrt(abs(y-.01*x**2))+.01*abs(x+10)
twoD_domain = np.linspace(-20,20)
combo_domain = list(itertools.product(*[twoD_domain,twoD_domain]))
xs = [combo_domain[i][0] for i in range(len(combo_domain))]
ys = [combo_domain[i][1] for i in range(len(combo_domain))]
twoD_function_y = [twoD_function(combo_domain[i][0],combo_domain[i][1]) for i in range(len(combo_domain))]
fig = plt.figure(figsize=(13,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xs, ys, twoD_function_y,cmap='hot',c=twoD_function_y);
%autoreload 2
bo_implementation = IBO()
test_x = np.array(np.linspace(-10,10,1000)).reshape(-1,1)
train_x = np.array(np.random.choice(test_x.ravel())).reshape(-1,1)
train_y_numbers = np.array([oneD_function(_) for _ in train_x]).reshape(-1,1)<jupyter_output><empty_output><jupyter_text>One Dimensional<jupyter_code>bo_implementation.fit(train_x,train_y_numbers, test_x, oneD_function, y_func_type='real' ,samples = 10)
# one sample
x_sampled_points, y_sampled_points, best_x, best_y = bo_implementation.predict()
# find the max
bo_implementation.maximize(n_steps=10)<jupyter_output>Bayesian Optimization just sampled point = [-3.05305305]
Best x (Bayesian Optimization) = [-3.05305305],                         Best y = [ 279.90363718]
<class 'numpy.ndarray'> type testing domain
Bayesian Optimization just sampled point = [-4.09409409]
Best x (Bayesian Optimization) = [-4.09409409], Best y = [ 320.07987715]
<class 'numpy.ndarray'> type testing domain
Bayesian Optimization just sampled point = [-3.91391391]
Best x (Bayesian Optimization) = [-4.09409409], Best y = [ 320.07987715]
<class 'numpy.ndarray'> type testing domain
Bayesian Optimization just sampled point = [-5.45545546]
Best x (Bayesian Optimization) = [-4.09409409], Best y = [ 320.07987715]
<class 'numpy.ndarray'> type testing domain
Bayesian Optimization just sampled point = [-4.43443443]
Best x (Bayesian Optimization) = [-4.09409409], Best y = [ 320.07987715]
<class 'numpy.ndarray'> type testing domain
Bayesian Optimization just sampled point = [-4.31431431]
Best x (Bayesian Optimization) = [-4.31431431], B[...]<jupyter_text>Two dimensional real function<jupyter_code>%autoreload 2
bo_implementation_2d = IBO()
bo_implementation_2d .fit(train_x_2d, train_y_2d, combo_domain, twoD_function, y_func_type='real' ,samples = 10)
x_sampled_points, y_sampled_points, sampled_var, \
            best_x, best_y, improvements, domain, mus  = next(bo_implementation_2d.bo_gen)
plt.plot(improvements)
x_sampled_points, y_sampled_points, best_x, best_y = bo_implementation_2d.predict()
 bo_implementation_2d.maximize(n_steps=10)
twoD_domain = np.linspace(-20,20)
combo_domain = np.squeeze(np.array([list(itertools.product(*[twoD_domain,twoD_domain]))]))
xs = [combo_domain[i][0] for i in range(len(combo_domain))]
ys = [combo_domain[i][1] for i in range(len(combo_domain))]
twoD_function_y = [twoD_function(combo_domain[i][0],combo_domain[i][1]) for i in range(len(combo_domain))]
train_x_2d = np.array([[1,4],
                      [3,5]])
train_y_2d = np.array([[twoD_function(i[0],i[1]) for i in train_x_2d]]).reshape(-1,1)
train_y_2d
train_x_2d
np.argmax(train_y_2d)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09576363241364327, 'fraction_numerical': 0.02337605909189659, 'mean_word_length': 2.6843284776692813, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 53, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '374430', 'license': 'MIT', 'n_tokens_mistral': 7871, 'n_tokens_neox': 7268, 'n_words': 1723, 'path': 'advanced_stats_bayesian_optimization/.ipynb_checkpoints/Compare_IBO_to_Implementation-checkpoint.ipynb', 'repo_name': 'jonhilgart22/glavanize-projects', 'size': '23027'}