new examples including gp2Scale posterior

examples/run_gp2ScaleCPU_local.py

@@ -92,7 +92,9 @@ def main():
     print("Log Likelihood: ", my_gp.log_likelihood(my_gp.hyperparameters))
     print("all done after: ",time.time() - st," seconds")
 
-    my_gp.train(hps_bounds, max_iter = 5, init_hyperparameters = init_hps)
+    my_gp.train(hps_bounds, max_iter = 2, init_hyperparameters = init_hps)
+    my_gp.posterior_mean(np.random.rand(2,3))
+    my_gp.posterior_covariance(np.random.rand(2,3))
 
 
 if __name__ == '__main__':

fvgp/gp2Scale.py

@@ -42,76 +42,6 @@ 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
 
 
 ############################################################
@@ -298,15 +228,12 @@ def compute_covariance(self, x1,x2,hyperparameters, variances,client):
         futures = []           ### a list of futures
         actor_futures = []
         finished_futures = set()
-        #scatter_future = []
         ###get workers
         compute_workers = set(self.compute_worker_set)
         idle_workers = set(compute_workers)
         ###future_worker_assignments
         self.future_worker_assignments = {}
         ###scatter data
-        #print("scatter future: ", scatter_future["hps"].result(), flush = True)
-        ###############
         start_time = time.time()
         count = 0
         num_batches_i = len(x1) // self.batch_size
@@ -339,10 +266,10 @@ def compute_covariance(self, x1,x2,hyperparameters, variances,client):
                 data = {"scattered_data": self.scatter_future,"hps": hyperparameters, "kernel" :self.kernel,  "range_i": (beg_i,end_i), "range_j": (beg_j,end_j), "mode": "prior","gpu": 0}
                 futures.append(client.submit(kernel_function, data, workers = current_worker))
                 self.assign_future_2_worker(futures[-1].key,current_worker)
-                #if self.info:
-                #    print("    submitted batch. i:", beg_i,end_i,"   j:",beg_j,end_j, "to worker ",current_worker, " Future: ", futures[-1].key)
-                #    print("    current time stamp: ", time.time() - start_time," percent finished: ",float(count)/self.total_number_of_batches())
-                #    print("")
+                if self.info:
+                    print("    submitted batch. i:", beg_i,end_i,"   j:",beg_j,end_j, "to worker ",current_worker, " Future: ", futures[-1].key)
+                    print("    current time stamp: ", time.time() - start_time," percent finished: ",float(count)/self.total_number_of_batches())
+                    print("")
                 count += 1
 
         if self.info:
@@ -567,6 +494,49 @@ def posterior_mean(self, x_iset):
         return {"x": p,
                 "f(x)": posterior_mean}
 
+
+    def posterior_covariance(self, x_iset, umfpack = True):
+        """
+        Function to compute the posterior covariance.
+        Parameters
+        ----------
+        x_iset : np.ndarray
+            A numpy array of shape (V x D), interpreted as  an array of input point positions.
+        variance_only : bool, optional
+            If True the compuation of the posterior covariance matrix is avoided which can save compute time.
+            In that case the return will only provide the variance at the input points.
+            Default = False.
+        Return
+        ------
+        solution dictionary : dict
+        """
+        p = np.array(x_iset)
+        if p.ndim == 1: p = np.array([p])
+        if len(p[0]) != len(self.x_data[0]): p = np.column_stack([p,np.zeros((len(p)))])
+        if len(p) == len(self.x_data): raise Warning("The number of prediction points and data points coincide which can lead to instabilities when user-defined kernels are used.")
+
+        k = self.kernel(self.x_data,p,self.hyperparameters)
+        kk = self.kernel(p, p,self.hyperparameters)
+        k_cov_prod = np.empty(k.shape)
+        for i in range(len(k[0])): k_cov_prod[:,i] = spsolve(self.SparsePriorCovariance.get_result().result(),k[:,i], use_umfpack = umfpack)
+
+        S = kk - (k_cov_prod.T @ k)
+        v = np.array(np.diag(S))
+        if np.any(v < -0.001):
+            logger.warning(inspect.cleandoc("""#
+            Negative variances encountered. That normally means that the model is unstable.
+            Rethink the kernel definitions, add more noise to the data,
+            or double check the hyperparameter optimization bounds. This will not
+            terminate the algorithm, but expect anomalies."""))
+            v[v<0.0] = 0.0
+            if not variance_only:
+                np.fill_diagonal(S, v)
+
+        return {"x": p,
+                "v(x)": v,
+                "S(x)": S}
+
+
 def kernel_function(data):
     st = time.time()
     hps= data["hps"]