US Topography

The data for this script can be downloaded here

This script uses version 6 of gpCAM

import numpy as np

from gpcam.gp_optimizer import GPOptimizer

import matplotlib.pyplot as plt

from numpy.random import default_rng


a = np.load("us_topo.npy")

rng = default_rng()

ind = rng.choice(len(a)-1, size=3000, replace=False)

points = a[ind,0:2]

values = a[ind,2:3]

print("x_min ", np.min(points[:,0])," x_max ",np.max(points[:,0]))

print("y_min ", np.min(points[:,1])," y_max ",np.max(points[:,1]))

print("length of data set: ", len(points))


index_set_bounds = np.array([[0,99],[0,248]])

hyperparameter_bounds = np.array([[0.001,1e9],[1,1000],[1,1000]])

hps_guess = np.array([4.71907062e+06, 4.07439017e+02, 3.59068120e+02])


###################################################################################

gp = GPOptimizer(2,1,1, index_set_bounds)

gp.tell(points,values)

gp.init_gp(hps_guess)

gp.train_gp(hyperparameter_bounds,likelihood_optimization_pop_size = 20,

likelihood_optimization_tolerance = 1e-6, likelihood_optimization_max_iter = 2)


x_pred = np.empty((10000,2))

counter = 0

x = np.linspace(0,99,100)

y = np.linspace(0,248,100)


for i in x:

for j in y:

x_pred[counter] = np.array([i,j])

counter += 1


res1 = gp.gp.posterior_mean(x_pred)

res2 = gp.gp.posterior_covariance(x_pred)

#res3 = gp.gp.shannon_information_gain(x_pred)

X,Y = np.meshgrid(x,y)


PM = np.reshape(res1["f(x)"],(100,100))

PV = np.reshape(res2["v(x)"],(100,100))

plt.figure(figsize= (10,10))

plt.pcolormesh(X,Y,PM)

plt.figure(figsize= (10,10))

plt.pcolormesh(X,Y,PV)

plt.show()

next = gp.ask(position = None, n = 1, objective_function = "covariance", optimization_bounds = None,

optimization_method = "global", optimization_pop_size = 50, optimization_max_iter = 20,

optimization_tol = 10e-6, dask_client = False)

print(next)