more restructuring

fvgp/gp.py

@@ -30,6 +30,7 @@
 #   neither minres nor random logdet are doing a good job, cg is better but we need a preconditioner, maybe a large LU?
 #   when using gp2Scale but the solution is dense we should go with dense linalg [Done]
 #   init hps and bounds should only be initiated for the default kernel
+#   the mcmc in default mode should not need proposal distributions explicitly
 
 class GP:
     """
@@ -1221,7 +1222,8 @@ def _compute_gp_linalg(self, vec, KV, calc_inv=False):
             elif mode == "dense_Chol":
                 KV = KV.toarray()
                 KVinvY, KVlogdet, factorization_obj = self._Chol(KV, vec)
-            else: KVinvY, KVlogdet, factorization_obj = self._gp2Scale_linalg(KV, vec)
+            elif mode == "gp2Scale": KVinvY, KVlogdet, factorization_obj = self._gp2Scale_linalg(KV, vec)
+            else: raise Exception("No linear algebra mode applicable in '_compute_gp_linalg'")
             KVinv = None
         else:
             if calc_inv:
@@ -1262,11 +1264,13 @@ def _LU(self, KV, vec):
         return KVinvY, KVlogdet, factorization_obj
 
     def _Chol(self, KV, vec):
+        if self.info: print("Dense Cholesky in progress ...")
         c, l = cho_factor(KV)
         factorization_obj = ("Chol", c, l)
         KVinvY = cho_solve((c, l), vec)
         upper_diag = abs(c.diagonal())
         KVlogdet = 2.0 * np.sum(np.log(upper_diag))
+        if self.info: print("Dense Cholesky compute time: ", time.time() - st, "seconds.")
         return KVinvY, KVlogdet, factorization_obj
 
     def _gp2Scale_linalg(self, KV, vec):

fvgp/gp_kernels.py

@@ -341,33 +341,6 @@ def polynomial_kernel(x1, x2, p):
     return p
 
 
-def default_kernel(x1, x2, hyperparameters, obj):
-    """
-    Function for the default kernel, a Matern kernel of first-order differentiability.
-
-    Parameters
-    ----------
-    x1 : np.ndarray
-        Numpy array of shape (U x D).
-    x2 : np.ndarray
-        Numpy array of shape (V x D).
-    hyperparameters : np.ndarray
-        Array of hyperparameters. For this kernel we need D + 1 hyperparameters.
-    obj : object instance
-        GP object instance.
-
-    Return
-    ------
-    Covariance matrix : np.ndarray
-    """
-    hps = hyperparameters
-    distance_matrix = np.zeros((len(x1), len(x2)))
-    for i in range(len(x1[0])):
-        distance_matrix += abs(np.subtract.outer(x1[:, i], x2[:, i]) / hps[1 + i]) ** 2
-    distance_matrix = np.sqrt(distance_matrix)
-    return hps[0] * obj.matern_kernel_diff1(distance_matrix, 1)
-
-
 def wendland_anisotropic(x1, x2, hyperparameters, obj):
     """
     Function for the Wendland kernel.
@@ -506,6 +479,14 @@ def wendland_anisotropic_gp2Scale_cpu(x1, x2, hps, obj):
     return kernel
 
 
+def _get_distance_matrix_gpu(x1, x2, device, hps):  # pragma: no cover
+    import torch
+    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]) / hps[1 + i]) ** 2
+    return torch.sqrt(d)
+
+
 def wendland_anisotropic_gp2Scale_gpu(x1, x2, hps, obj):  # pragma: no cover
     import torch
     cuda_device = torch.device("cuda:0")

fvgp/gp_likelihood.py

@@ -2,8 +2,8 @@
 
 
 class GPlikelihood:  # pragma: no cover
-    def __init__(self, K, noise_variances, gp_noise_function):
-        assert isinstance(K, np.ndarray) and np.ndim(K) == 2
+    def __init__(self, prior_obj, noise_variances, gp_noise_function):
+        self.prior_obj = prior_obj
 
         if noise_variances is not None and callable(gp_noise_function):
             warnings.warn("Noise function and measurement noise provided. noise_variances set to None.", stacklevel=2)
@@ -31,22 +31,6 @@ def __init__(self, K, noise_variances, gp_noise_function):
                 self.noise_function_grad = self._default_dnoise_dh
 
     ##################################################################################
-    def _KVsolve(self, b):
-        if self.factorization_obj[0] == "LU":
-            LU = self.factorization_obj[1]
-            return LU.solve(b)
-        elif self.factorization_obj[0] == "Chol":
-            c, l = self.factorization_obj[1], self.factorization_obj[2]
-            return cho_solve((c, l), b)
-        else:
-            res = np.empty((len(b), b.shape[1]))
-            if b.shape[1] > 100: warnings.warn(
-                "You want to predict at >100 points. \n When using gp2Scale, this takes a while. \n \
-                Better predict at only a handful of points.")
-            for i in range(b.shape[1]):
-                res[:, i], exit_status = minres(self.KV, b[:, i])
-                # if exit_status != 0: res[:,i], exit_status = minres(self.KV,b[:,i])
-            return res
 
     def log_likelihood(self, hyperparameters=None):
         """
@@ -211,3 +195,35 @@ def _default_noise_function(self, x, hyperparameters, gp_obj):
             return self.gp2Scale_obj.calculate_sparse_noise_covariance(noise)
         else:
             return np.diag(noise)
+
+    def _default_dnoise_dh(self, x, hps, gp_obj):
+        gr = np.zeros((len(hps), len(x), len(x)))
+        return gr
+
+    def _default_dnoise_dh_econ(self, x, i, hps, gp_obj):
+        gr = np.zeros((len(x), len(x)))
+        return gr
+
+    ##########################
+    def _finitediff_dnoise_dh(self, x, hps, gp_obj):
+        gr = np.zeros((len(hps), len(x), len(x)))
+        for i in range(len(hps)):
+            temp_hps1 = np.array(hps)
+            temp_hps1[i] = temp_hps1[i] + 1e-6
+            temp_hps2 = np.array(hps)
+            temp_hps2[i] = temp_hps2[i] - 1e-6
+            a = self.noise_function(x, temp_hps1, self)
+            b = self.noise_function(x, temp_hps2, self)
+            gr[i] = (a - b) / 2e-6
+        return gr
+
+    ##########################
+    def _finitediff_dnoise_dh_econ(self, x, i, hps, gp_obj):
+        temp_hps1 = np.array(hps)
+        temp_hps1[i] = temp_hps1[i] + 1e-6
+        temp_hps2 = np.array(hps)
+        temp_hps2[i] = temp_hps2[i] - 1e-6
+        a = self.noise_function(x, temp_hps1, self)
+        b = self.noise_function(x, temp_hps2, self)
+        gr = (a - b) / 2e-6
+        return gr

fvgp/gp_marginal_density.py

@@ -0,0 +1,98 @@
+import numpy as np
+
+
+class GPMarginalDensity:
+    def __init__(self, data_obj, prior_obj, likelihood_obj):
+        self.data_obj = data_obj
+        self.prior_obj = prior_obj
+        self.likelihood_obj = likelihood_obj
+
+    def _compute_GPpriorV(self, x_data, y_data, hyperparameters, calc_inv=False):
+        # get the prior mean
+        prior_mean_vec = self.mean_function(x_data, hyperparameters, self)
+        assert np.ndim(prior_mean_vec) == 1
+        # get the latest noise
+        V = self.noise_function(x_data, hyperparameters, self)
+        assert np.ndim(V) == 2
+        K = self._compute_K(x_data, x_data, hyperparameters)
+        # check if shapes are correct
+        if K.shape != V.shape: raise Exception("Noise covariance and prior covariance not of the same shape.")
+        # get K + V
+        KV = K + V
+
+        # get Kinv/KVinvY, LU, Chol, logdet(KV)
+        KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(y_data-prior_mean_vec, KV,
+                                                calc_inv=calc_inv)
+        return K, KV, KVinvY, KVlogdet, factorization_obj, KVinv, prior_mean_vec, V
+
+    ##################################################################################
+
+    def _update_GPpriorV(self, x_data_old, x_new, y_data, hyperparameters, calc_inv=False):
+        # get the prior mean
+        prior_mean_vec = np.append(self.prior_mean_vec, self.mean_function(x_new, hyperparameters, self))
+        assert np.ndim(prior_mean_vec) == 1
+        # get the latest noise
+        V = self.noise_function(self.data.x_data, hyperparameters, self) #can be avoided by update
+        assert np.ndim(V) == 2
+        # get K
+        K = self.update_K(x_data_old, x_new, hyperparameters)
+
+        # check if shapes are correct
+        if K.shape != V.shape: raise Exception("Noise covariance and prior covariance not of the same shape.")
+
+        # get K + V
+        KV = K + V
+        # get Kinv/KVinvY, LU, Chol, logdet(KV)
+
+        if self.online_mode is True: KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(
+            y_data - prior_mean_vec, k, kk)
+        else: KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(y_data-prior_mean_vec, KV,
+                                                                             calc_inv=calc_inv)
+        return K, KV, KVinvY, KVlogdet, factorization_obj, KVinv, prior_mean_vec, V
+
+    ##################################################################################
+    def _compute_gp_linalg(self, vec, KV, calc_inv=False):
+        if calc_inv:
+            KVinv = self._inv(KV)
+            factorization_obj = ("Inv", None)
+            KVinvY = KVinv @ vec
+            KVlogdet = self._logdet(KV)
+        else:
+            KVinv = None
+            KVinvY, KVlogdet, factorization_obj = self._Chol(KV, vec)
+        return KVinvY, KVlogdet, factorization_obj, KVinv
+
+    def _update_gp_linalg(self, vec, k, kk):
+        X = self._inv(kk - C @ self.KVinv @ B)
+        F = -self.KVinv @ k @ X
+        KVinv = np.block([[self.KVinv + self.KVinv @ B @ X @ C @ self.KVinv, F],
+                          [F.T,                                              X]])
+        factorization_obj = ("Inv", None)
+        KVinvY = KVinv @ vec
+        KVlogdet = self.KVlogdet + self._logdet(kk - k.T @ self.KVinv @ k)
+        return KVinvY, KVlogdet, factorization_obj, KVinv
+
+    def _LU(self, KV, vec):
+        st = time.time()
+        if self.info: print("LU in progress ...")
+        LU = splu(KV.tocsc())
+        factorization_obj = ("LU", LU)
+        KVinvY = LU.solve(vec)
+        upper_diag = abs(LU.U.diagonal())
+        KVlogdet = np.sum(np.log(upper_diag))
+        if self.info: print("LU compute time: ", time.time() - st, "seconds.")
+        return KVinvY, KVlogdet, factorization_obj
+
+    def _Chol(self, KV, vec):
+        if self.info: print("Dense Cholesky in progress ...")
+        c, l = cho_factor(KV)
+        factorization_obj = ("Chol", c, l)
+        KVinvY = cho_solve((c, l), vec)
+        upper_diag = abs(c.diagonal())
+        KVlogdet = 2.0 * np.sum(np.log(upper_diag))
+        if self.info: print("Dense Cholesky compute time: ", time.time() - st, "seconds.")
+        return KVinvY, KVlogdet, factorization_obj
+
+    ##################################################################################
+
+

fvgp/gp_posterior.py

@@ -1,7 +1,30 @@
-
-class GPosterior: # pragma: no cover
-    def __init__(self, KVinvY=None):
+class GPosterior:  # pragma: no cover
+    def __init__(self, prior_obj, data_obj):
         assert isinstance(KVinvY, np.ndarray)
+        self.KV = prior_obj.KV
+        self.factorization_obj = prior_obj.factorization_obj
+        self.prior_obj = prior_obj
+        self.data_obj = data_obj
+
+    def _KVsolve(self, b):
+        if self.factorization_obj[0] == "LU":
+            LU = self.factorization_obj[1]
+            return LU.solve(b)
+        elif self.factorization_obj[0] == "Chol":
+            c, l = self.factorization_obj[1], self.factorization_obj[2]
+            return cho_solve((c, l), b)
+        elif self.factorization_obj[0] == "gp2Scale":
+            res = np.empty((len(b), b.shape[1]))
+            if b.shape[1] > 100: warnings.warn(
+                "You want to predict at >100 points. \n When using gp2Scale, this takes a while. \n \
+                Better predict at only a handful of points.")
+            for i in range(b.shape[1]):
+                res[:, i], exit_status = cg(self.KV, b[:, i])
+            return res
+        elif self.factorization_obj[0] == "Inv":
+            return self.KVinv @ b
+        else:
+            raise Exception("Non-permitted factorization object encountered.")
 
     def posterior_mean(self, x_pred, hyperparameters=None, x_out=None):
         """