work on gp2Scale

examples/run_gp2ScaleCPU.py

@@ -33,55 +33,6 @@ def client_is_ready(ip, n_workers):
         if counter > 20: print("getting the client is taking a long time: ", counter * 5, "seconds", flush = True)
 
 
-def kernel(x1,x2, hps,obj):
-    def b(x,x0,r = 0.1, ampl = 1.0):
-        """
-        evaluates the bump function
-        x ... a point (1d numpy array)
-        x0 ... 1d numpy array of location of bump function
-        returns the bump function b(x,x0) with radius r
-        """
-        x_new = x - x0
-        d = np.linalg.norm(x_new, axis = 1)
-        a = np.zeros(d.shape)
-        a = 1.0 - (d**2/r**2)
-        i = np.where(a > 0.0)
-        bump = np.zeros(a.shape)
-        bump[i] = ampl * np.exp((-1.0/a[i])+1)
-        return bump
-
-
-    def f(x,x0, radii, amplts):
-        b1 = b(x, x0[0:3],r = radii[0], ampl = amplts[0])  ###x0[0] ... D-dim location of bump func 1
-        b2 = b(x, x0[3:6],r = radii[1], ampl = amplts[1])  ###x0[1] ... D-dim location of bump func 2
-        b3 = b(x, x0[6:9],r = radii[2], ampl = amplts[2])  ###x0[1] ... D-dim location of bump func 2
-        b4 = b(x, x0[9:12],r = radii[3], ampl = amplts[3])  ###x0[1] ... D-dim location of bump func 2
-        return b1 + b2 + b3 + b4
-
-    def g(x,x0, radii, amplts):
-        b1 = b(x, x0[0:3],r = radii[0], ampl = amplts[0])  ###x0[0] ... D-dim location of bump func 1
-        b2 = b(x, x0[3:6],r = radii[1], ampl = amplts[1])  ###x0[1] ... D-dim location of bump func 2
-        b3 = b(x, x0[6:9],r = radii[2], ampl = amplts[2])  ###x0[1] ... D-dim location of bump func 2
-        b4 = b(x, x0[9:12],r = radii[3], ampl = amplts[3])  ###x0[1] ... D-dim location of bump func 2
-        return b1 + b2 + b3 + b4
-    def sparse_stat_kernel(x1,x2, hps):
-        d = 0
-        for i in range(len(x1[0])):
-            d += np.subtract.outer(x1[:,i],x2[:,i])**2
-        d = np.sqrt(d)
-        d[d == 0.0] = 1e-6
-        d[d > hps] = hps
-        kernel = (np.sqrt(2.0)/(3.0*np.sqrt(np.pi)))*\
-        ((3.0*(d/hps)**2*np.log((d/hps)/(1+np.sqrt(1.0 - (d/hps)**2))))+\
-        ((2.0*(d/hps)**2+1.0)*np.sqrt(1.0-(d/hps)**2)))
-
-        return kernel
-
-    k = np.outer(f(x1,hps[0:12],hps[12:16],hps[16:20]),
-                 f(x2,hps[0:12],hps[12:16],hps[16:20])) + \
-        np.outer(g(x1,hps[20:32],hps[32:36],hps[36:40]),
-                 g(x2,hps[20:32],hps[32:36],hps[36:40]))
-    return k + hps[40] * sparse_stat_kernel(x1,x2, hps[41])
 
 def normalize(v):
     v = v - np.min(v)

fvgp/advanced_kernels.py

@@ -0,0 +1,166 @@
+import time
+import numpy as np
+import torch
+
+
+###################################
+######gp2Scale GPU kernels#########
+###################################
+def sparse_stat_kernel(x1,x2, hps):
+    d = 0
+    for i in range(len(x1[0])): d += abs(np.subtract.outer(x1[:,i],x2[:,i]))**2
+    d = np.sqrt(d)
+    d[d == 0.0] = 1e-16
+    d[d > hps[1]] = hps[1]
+    kernel = (np.sqrt(2.0)/(3.0*np.sqrt(np.pi)))*\
+    ((3.0*(d/hps[1])**2*np.log((d/hps[1])/(1+np.sqrt(1.0 - (d/hps[1])**2))))+\
+    ((2.0*(d/hps[1])**2 + 1.0) * np.sqrt(1.0-(d/hps[1])**2)))
+    return hps[0] * kernel
+
+def sparse_stat_kernel_robust(x1,x2, hps):
+    d = 0
+    for i in range(len(x1[0])): d += abs(np.subtract.outer(x1[:,i],x2[:,i]))**2
+    d = np.sqrt(d)
+    d[d == 0.0] = 1e-16
+    d[d > 1./hps[1]**2] = 1./hps[1]**2
+    kernel = (np.sqrt(2.0)/(3.0*np.sqrt(np.pi)))*\
+    ((3.0*(d*hps[1]**2)**2*np.log((d*hps[1]**2)/(1+np.sqrt(1.0 - (d*hps[1]**2)**2))))+\
+    ((2.0*(d*hps[1]**2)**2 + 1.0) * np.sqrt(1.0-(d*hps[1]**2)**2)))
+    return (hps[0]**2) * kernel
+
+
+def f_gpu(x,x0, radii, amplts, device):
+    b1 = b_gpu(x, x0[0:3], radii[0], amplts[0], device)  ###x0[0] ... D-dim location of bump func 1
+    b2 = b_gpu(x, x0[3:6], radii[1], amplts[1], device)  ###x0[1] ... D-dim location of bump func 2
+    b3 = b_gpu(x, x0[6:9], radii[2], amplts[2], device)  ###x0[1] ... D-dim location of bump func 2
+    b4 = b_gpu(x, x0[9:12],radii[3], amplts[3], device)  ###x0[1] ... D-dim location of bump func 2
+    return b1 + b2 + b3 + b4
+
+def g_gpu(x,x0, radii, amplts,device):
+    b1 = b_gpu(x, x0[0:3], radii[0], amplts[0], device)  ###x0[0] ... D-dim location of bump func 1
+    b2 = b_gpu(x, x0[3:6], radii[1], amplts[1], device)  ###x0[1] ... D-dim location of bump func 2
+    b3 = b_gpu(x, x0[6:9], radii[2], amplts[2], device)  ###x0[1] ... D-dim location of bump func 2
+    b4 = b_gpu(x, x0[9:12],radii[3], amplts[3], device)  ###x0[1] ... D-dim location of bump func 2
+    return b1 + b2 + b3 + b4
+
+
+def b_gpu(x,x0, r, ampl, device):
+    """
+    evaluates the bump function
+    x ... a point (1d numpy array)
+    x0 ... 1d numpy array of location of bump function
+    returns the bump function b(x,x0) with radius r
+    """
+    x_new = x - x0
+    d = torch.linalg.norm(x_new, axis = 1)
+    a = torch.zeros(d.shape).to(device, dtype = torch.float32)
+    a = 1.0 - (d**2/r**2)
+    i = torch.where(a > 0.0)
+    bump = torch.zeros(a.shape).to(device, dtype = torch.float32)
+    e = torch.exp((-1.0/a[i])+1).to(device, dtype = torch.float32)
+    bump[i] = ampl * e
+    return bump
+
+
+def get_distance_matrix_gpu(x1,x2,device):
+    d = torch.zeros((len(x1),len(x2))).to(device, dtype = torch.float32)
+    for i in range(x1.shape[1]):
+        d += ((x1[:,i].reshape(-1, 1) - x2[:,i]))**2
+    return torch.sqrt(d)
+
+def sparse_stat_kernel_gpu(x1,x2, hps,device):
+    d = get_distance_matrix_gpu(x1,x2,device)
+    d[d == 0.0] = 1e-6
+    d[d > hps] = hps
+
+    d_hps = d/hps
+    d_hpss= d_hps**2
+    sq = torch.sqrt(1.0 - d_hpss)
+
+
+    kernel = (torch.sqrt(torch.tensor(2.0))/(3.0*torch.sqrt(torch.tensor(3.141592653))))*\
+    ((3.0*d_hpss * torch.log((d_hps)/(1+sq)))+\
+    ((2.0*d_hpss + 1.0)*sq))
+    return kernel
+
+
+def ks_gpu(x1,x2,hps,cuda_device):
+    k1 = torch.outer(f(x1,hps[0:12],hps[12:16],hps[16:20],cuda_device),
+                        f(x2,hps[0:12],hps[12:16],hps[16:20],cuda_device)) + \
+            torch.outer(g(x1,hps[20:32],hps[32:36],hps[36:40],cuda_device),
+                        g(x2,hps[20:32],hps[32:36],hps[36:40],cuda_device))
+    k2 = sparse_stat_kernel_gpu(x1,x2, hps[41],cuda_device)
+    return k1 + hps[40]*k2
+
+def kernel_gpu(x1,x2, hps, info = None):
+    cuda_device = torch.device("cuda:0")
+    x1_dev = torch.from_numpy(x1).to(cuda_device, dtype = torch.float32)
+    x2_dev = torch.from_numpy(x2).to(cuda_device, dtype = torch.float32)
+    hps_dev = torch.from_numpy(hps).to(cuda_device, dtype = torch.float32)
+    ksparse = ks_gpu(x1_dev,x2_dev,hps_dev,cuda_device).cpu().numpy()
+    return ksparse
+
+
+###################################
+######gp2Scale CPU kernels#########
+###################################
+
+def b_cpu(x,x0,r = 0.1, ampl = 1.0):
+    """
+    evaluates the bump function
+    x ... a point (1d numpy array)
+    x0 ... 1d numpy array of location of bump function
+    returns the bump function b(x,x0) with radius r
+    """
+    x_new = x - x0
+    d = np.linalg.norm(x_new, axis = 1)
+    a = np.zeros(d.shape)
+    a = 1.0 - (d**2/r**2)
+    i = np.where(a > 0.0)
+    bump = np.zeros(a.shape)
+    bump[i] = ampl * np.exp((-1.0/a[i])+1)
+    return bump
+
+
+def f_cpu(x,x0, radii, amplts):
+    b1 = b_cpu(x, x0[0:3],r = radii[0], ampl = amplts[0])  ###x0[0] ... D-dim location of bump func 1
+    b2 = b_cpu(x, x0[3:6],r = radii[1], ampl = amplts[1])  ###x0[1] ... D-dim location of bump func 2
+    b3 = b_cpu(x, x0[6:9],r = radii[2], ampl = amplts[2])  ###x0[1] ... D-dim location of bump func 2
+    b4 = b_cpu(x, x0[9:12],r = radii[3], ampl = amplts[3])  ###x0[1] ... D-dim location of bump func 2
+    return b1 + b2 + b3 + b4
+
+def g_cpu(x,x0, radii, amplts):
+    b1 = b_cpu(x, x0[0:3],r = radii[0], ampl = amplts[0])  ###x0[0] ... D-dim location of bump func 1
+    b2 = b_cpu(x, x0[3:6],r = radii[1], ampl = amplts[1])  ###x0[1] ... D-dim location of bump func 2
+    b3 = b_cpu(x, x0[6:9],r = radii[2], ampl = amplts[2])  ###x0[1] ... D-dim location of bump func 2
+    b4 = b_cpu(x, x0[9:12],r = radii[3], ampl = amplts[3])  ###x0[1] ... D-dim location of bump func 2
+    return b1 + b2 + b3 + b4
+def sparse_stat_kernel_cpu(x1,x2, hps):
+    d = 0
+    for i in range(len(x1[0])): d += np.subtract.outer(x1[:,i],x2[:,i])**2
+    d = np.sqrt(d)
+    d[d == 0.0] = 1e-6
+    d[d > hps] = hps
+    kernel = (np.sqrt(2.0)/(3.0*np.sqrt(np.pi)))*\
+    ((3.0*(d/hps)**2*np.log((d/hps)/(1+np.sqrt(1.0 - (d/hps)**2))))+\
+    ((2.0*(d/hps)**2+1.0)*np.sqrt(1.0-(d/hps)**2)))
+    return kernel
+
+
+def kernel_cpu(x1,x2, hps,info = None):
+        k = np.outer(f_cpu(x1,hps[0:12],hps[12:16],hps[16:20]),
+                 f_cpu(x2,hps[0:12],hps[12:16],hps[16:20])) + \
+        np.outer(g_cpu(x1,hps[20:32],hps[32:36],hps[36:40]),
+                 g_cpu(x2,hps[20:32],hps[32:36],hps[36:40]))
+    return k + hps[40] * sparse_stat_kernel(x1,x2, hps[41])
+
+
+############################################################
+############################################################
+############################################################
+############################################################
+############################################################
+############################################################
+############################################################
+
+

fvgp/gp.py

@@ -652,7 +652,10 @@ def test_log_likelihood_gradient(self,hyperparameters):
             thps_aux[i] = thps_aux[i] + eps
             grad[i] = (self.log_likelihood(thps_aux) - self.log_likelihood(thps))/eps
         analytical = self.log_likelihood_gradient(thps)
-        if np.linalg.norm(grad-analytical) > 1e-4: print("Gradient possibly wrong")
+        if np.linalg.norm(grad-analytical) > 1e-1: 
+            print("Gradient possibly wrong")
+            print(grad)
+            print(analytical)
         return grad, analytical
     ##################################################################################
     ##################################################################################
@@ -1724,16 +1727,6 @@ def non_stat_kernel(self,x1,x2,x0,w,l):
         non_stat = np.outer(self._g(x1,x0,w,l),self._g(x2,x0,w,l))
         return non_stat
 
-    def log_likelihood_gradient_test(self,hps):
-        eps = 1e-6
-        print("finite difference gradient:")
-        for i in range(len(hps)):
-            hps_a = np.array(hps)
-            hps_a[i] = hps_a[i] + eps
-            print((self.log_likelihood(hps_a) - self.log_likelihood(hps))/eps)
-        print("analytical gradient: ")
-        print(self.log_likelihood_gradient(hps))
-
 
     def non_stat_kernel_gradient(self,x1,x2,x0,w,l):
         dkdw = np.einsum('ij,k->ijk', self._dgdw(x1,x0,w,l), self._g(x2,x0,w,l)) + np.einsum('ij,k->ikj', self._dgdw(x2,x0,w,l), self._g(x1,x0,w,l))
@@ -1751,7 +1744,7 @@ def _get_distance_matrix(self,x1,x2):
     def _g(self,x,x0,w,l):
         d = self._get_distance_matrix(x,x0)
         e = np.exp( -(d**2) / l)
-        return  np.sum(w * e,axis = 1)
+        return np.sum(w * e,axis = 1)
 
     def _dgdw(self,x,x0,w,l):
         d = self._get_distance_matrix(x,x0)

fvgp/gp2Scale.py

@@ -42,6 +42,77 @@ def sparse_stat_kernel_robust(x1,x2, hps):
     ((2.0*(d*hps[1]**2)**2 + 1.0) * np.sqrt(1.0-(d*hps[1]**2)**2)))
     return (hps[0]**2) * kernel
 
+def kernel_gpu(x1,x2, hps):
+
+    def b(x,x0, r, ampl, device):
+        """
+        evaluates the bump function
+        x ... a point (1d numpy array)
+        x0 ... 1d numpy array of location of bump function
+        returns the bump function b(x,x0) with radius r
+        """
+        x_new = x - x0
+        d = torch.linalg.norm(x_new, axis = 1)
+        a = torch.zeros(d.shape).to(device, dtype = torch.float32)
+        a = 1.0 - (d**2/r**2)
+        i = torch.where(a > 0.0)
+        bump = torch.zeros(a.shape).to(device, dtype = torch.float32)
+        e = torch.exp((-1.0/a[i])+1).to(device, dtype = torch.float32)
+        bump[i] = ampl * e
+        return bump
+
+
+    def f(x,x0, radii, amplts, device):
+        b1 = b(x, x0[0:3], radii[0], amplts[0], device)  ###x0[0] ... D-dim location of bump func 1
+        b2 = b(x, x0[3:6], radii[1], amplts[1], device)  ###x0[1] ... D-dim location of bump func 2
+        b3 = b(x, x0[6:9], radii[2], amplts[2], device)  ###x0[1] ... D-dim location of bump func 2
+        b4 = b(x, x0[9:12],radii[3], amplts[3], device)  ###x0[1] ... D-dim location of bump func 2
+        return b1 + b2 + b3 + b4
+
+    def g(x,x0, radii, amplts,device):
+        b1 = b(x, x0[0:3], radii[0], amplts[0], device)  ###x0[0] ... D-dim location of bump func 1
+        b2 = b(x, x0[3:6], radii[1], amplts[1], device)  ###x0[1] ... D-dim location of bump func 2
+        b3 = b(x, x0[6:9], radii[2], amplts[2], device)  ###x0[1] ... D-dim location of bump func 2
+        b4 = b(x, x0[9:12],radii[3], amplts[3], device)  ###x0[1] ... D-dim location of bump func 2
+        return b1 + b2 + b3 + b4
+
+    def get_distance_matrix(x1,x2,device):
+        d = torch.zeros((len(x1),len(x2))).to(device, dtype = torch.float32)
+        for i in range(x1.shape[1]):
+            d += ((x1[:,i].reshape(-1, 1) - x2[:,i]))**2
+        return torch.sqrt(d)
+
+    def sparse_stat_kernel(x1,x2, hps,device):
+        d = get_distance_matrix(x1,x2,device)
+        d[d == 0.0] = 1e-6
+        d[d > hps] = hps
+
+        d_hps = d/hps
+        d_hpss= d_hps**2
+        sq = torch.sqrt(1.0 - d_hpss)
+
+
+        kernel = (torch.sqrt(torch.tensor(2.0))/(3.0*torch.sqrt(torch.tensor(3.141592653))))*\
+        ((3.0*d_hpss * torch.log((d_hps)/(1+sq)))+\
+        ((2.0*d_hpss + 1.0)*sq))
+        return kernel
+
+
+    def ks_gpu(x1,x2,hps,cuda_device):
+        k1 = torch.outer(f(x1,hps[0:12],hps[12:16],hps[16:20],cuda_device),
+                         f(x2,hps[0:12],hps[12:16],hps[16:20],cuda_device)) + \
+             torch.outer(g(x1,hps[20:32],hps[32:36],hps[36:40],cuda_device),
+                         g(x2,hps[20:32],hps[32:36],hps[36:40],cuda_device))
+        k2 = sparse_stat_kernel(x1,x2, hps[41],cuda_device)
+        return k1 + hps[40]*k2
+
+    cuda_device = torch.device("cuda:0")
+    x1_dev = torch.from_numpy(x1).to(cuda_device, dtype = torch.float32)
+    x2_dev = torch.from_numpy(x2).to(cuda_device, dtype = torch.float32)
+    hps_dev = torch.from_numpy(hps).to(cuda_device, dtype = torch.float32)
+    ksparse = ks_gpu(x1_dev,x2_dev,hps_dev,cuda_device).cpu().numpy()
+    return ksparse
+
 
 ############################################################
 ############################################################
@@ -202,7 +273,6 @@ def compute_prior_fvGP_pdf(self,client):
         self.covariance_value_prod = cov_y
 
     def _compute_covariance_value_product(self, hyperparameters,values, variances, mean, client):
-        np.save("latest_hps", hyperparameters)
         self.compute_covariance(self.x_data,self.x_data,hyperparameters, variances,client)
         y = values - mean
         #try: success = self.SparsePriorCovariance.compute_LU().result(timeout=self.LUtimeout)
@@ -391,6 +461,7 @@ def train(self,
             )
         print("computing the prior")
         self.compute_prior_fvGP_pdf(self.covariance_dask_client)
+        np.save("latest_hps", hyperparameters)
         ######################
         ######################
         ######################