more restructuring

fvgp/__init__.py

@@ -7,10 +7,11 @@
 from loguru import logger
 import sys
 from .gp import GP
+from .gp2 import GP as GP2
 from .fvgp import fvGP
 from .gpMCMC import gpMCMC
 
-__all__ = ['GP', 'fvGP']
+__all__ = ['GP','GP2', 'fvGP']
 __version__ = _version.get_versions()['version']
 
 logger.disable('fvgp')

fvgp/gp_data.py

@@ -32,24 +32,36 @@ def update(self, x_data_new, y_data_new, noise_variances_new=None, append=True):
         assert ((isinstance(noise_variances_new, np.ndarray) and np.ndim(noise_variances_new) == 1)
                 or noise_variances_new is None)
         if self.Euclidean: assert isinstance(x_data_new, np.ndarray) and np.ndim(x_data_new) == 2
-        else: assert isinstance(x_data_new, list)
-        if self.noise_variances and noise_variances_new is None:
-            raise Exception("Please provide noise_variances in the data update.")
+        else: assert (isinstance(x_data_new, list) and
+                      np.ndim(x_data_new) == 2 and
+                      self.input_space_dim == x_data_new.shape[1])
+
+        if self.noise_variances is not None and noise_variances_new is None:
+            print(self.noise_variances, noise_variances_new)
+            raise Exception("Please provide noise_variances in the data update because you did at initialization"
+                            "or during a previous update.")
+        if self.noise_variances is None and noise_variances_new is not None:
+            raise Exception("You did not initialize noise and but included noise in the update"
+                            "but this changes the settings. Please reinitialize in this case.")
+        if callable(noise_variances_new): raise Exception("The update noise_variances cannot be a callable.")
+        if noise_variances_new is not None:
+            assert isinstance(noise_variances_new, np.ndarray) and np.ndim(noise_variances_new) == 1
 
         if append is False:
             self.x_data = x_data_new
             self.y_data = y_data_new
+            self.noise_variances = noise_variances_new
         else:
-            if self.Euclidean:
-                np.row_stack([self.x_data, x_data_new])
-            else:
-                self.x_data = self.x_data + x_data_new
+            if self.Euclidean: self.x_data = np.row_stack([self.x_data, x_data_new])
+            else: self.x_data = self.x_data + x_data_new
             self.y_data = np.append(self.y_data, y_data_new)
-            if noise_variances_new is not None:
+            if isinstance(noise_variances_new, np.ndarray):
                 self.noise_variances = np.append(self.noise_variances, noise_variances_new)
         self.point_number = len(self.x_data)
         self._check_for_nan()
 
     def _check_for_nan(self):
         if self.Euclidean:
             if np.isnan(np.sum(self.x_data) + np.sum(self.y_data)): raise Exception("NaNs encountered in dataset.")
+
+

fvgp/gp_kernels.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def squared_exponential_kernel(distance, length):
     """
     Function for the squared exponential kernel.
@@ -418,6 +419,29 @@ def get_distance_matrix(x1, x2):
     return np.sqrt(d)
 
 
+def get_anisotropic_distance_matrix(x1, x2, hps):
+    """
+    Function to calculate the pairwise axial-anisotropic distance matrix of
+    points in x1 and x2.
+
+    Parameters
+    ----------
+    x1 : np.ndarray
+        Numpy array of shape (U x D).
+    x2 : np.ndarray
+        Numpy array of shape (V x D).
+    hps : np.ndarray
+        1d array of values. The diagonal of the metric tensor describing the axial anisotropy.
+
+    Return
+    ------
+    distance matrix : np.ndarray
+    """
+    d = np.zeros((len(x1), len(x2)))
+    for i in range(len(x1[0])): d += abs(np.subtract.outer(x1[:, i], x2[:, i]) / hps[i]) ** 2
+    return np.sqrt(d)
+
+
 def _g(x, x0, w, l):
     d = get_distance_matrix(x, x0)
     e = np.exp(-(d ** 2) / l)
@@ -463,3 +487,40 @@ def wendland_anisotropic_gp2Scale_gpu(x1, x2, hps, obj):  # pragma: no cover
     d[d > 1.] = 1.
     kernel = hps[0] * (1. - d) ** 8 * (35. * d ** 3 + 25. * d ** 2 + 8. * d + 1.)
     return kernel.cpu().numpy()
+
+
+def example_deep_kernel(x1, x2, hps, obj):
+    hps_nn_weights_l1 = hps[0:49]
+    hps_nn_weights_l2 = hps[49:98]
+    hps_nn_weights_l3 = hps[98:147]
+    hps_nn_weights_l4 = hps[147:196]
+    hps_nn_weights_l5 = hps[196:245]
+
+    hps_nn_biases_l1 = hps[245:252]
+    hps_nn_biases_l2 = hps[252:259]
+    hps_nn_biases_l3 = hps[259:266]
+    hps_nn_biases_l4 = hps[266:273]
+    hps_nn_biases_l5 = hps[273:280]
+
+    gp_deep_kernel_layer_width = 7
+    iset_dim = 7
+
+    n.set_weights(hps_nn_weights_l1.reshape(gp_deep_kernel_layer_width, iset_dim),
+                  hps_nn_weights_l2.reshape(gp_deep_kernel_layer_width, gp_deep_kernel_layer_width),
+                  hps_nn_weights_l3.reshape(gp_deep_kernel_layer_width, gp_deep_kernel_layer_width),
+                  hps_nn_weights_l4.reshape(gp_deep_kernel_layer_width, gp_deep_kernel_layer_width),
+                  hps_nn_weights_l5.reshape(iset_dim, gp_deep_kernel_layer_width)
+                  )
+    n.set_biases(hps_nn_biases_l1.reshape(gp_deep_kernel_layer_width),
+                 hps_nn_biases_l2.reshape(gp_deep_kernel_layer_width),
+                 hps_nn_biases_l3.reshape(gp_deep_kernel_layer_width),
+                 hps_nn_biases_l4.reshape(gp_deep_kernel_layer_width),
+                 hps_nn_biases_l5.reshape(iset_dim)
+                 )
+    signal_var = hps[280]
+    length_scale = hps[281]
+    x1_nn = n.forward(x1)
+    x2_nn = n.forward(x2)
+    d = obj.get_distance_matrix(x1_nn, x2_nn)
+    k = signal_var * obj.matern_kernel_diff1(d, length_scale)
+    return k

fvgp/gp_likelihood.py

@@ -2,7 +2,7 @@
 import warnings
 
 
-class GPlikelihood:  # pragma: no cover
+class GPlikelihood:
     def __init__(self,
                  x_data,
                  y_data,
@@ -28,7 +28,7 @@ def __init__(self,
         self.gp2Scale = gp2Scale
 
         if self.noise_variances is not None and callable(gp_noise_function):
-            raise Exception("Noise function and measurement noise provided. Decide which one to use.", stacklevel=2)
+            raise Exception("Noise function and measurement noise provided. Decide which one to use.")
         if callable(gp_noise_function):
             self.noise_function = gp_noise_function
         elif self.noise_variances is not None:
@@ -59,7 +59,10 @@ def update(self, x_data, y_data, noise_variances, hyperparameters):
         self.x_data = x_data
         self.y_data = y_data
         self.noise_variances = noise_variances
-        self.V = self.noise_function(self.x_data, hyperparameters, self)
+        self.V = self.calculate_V(hyperparameters)
+
+    def calculate_V(self, hyperparameters):
+        return self.noise_function(self.x_data, hyperparameters, self)
 
     def _default_noise_function(self, x, hyperparameters, gp_obj):
         noise = np.ones((len(x))) * (np.mean(abs(self.y_data)) / 100.0)

fvgp/gp_marginal_density.py

@@ -2,6 +2,12 @@
 from scipy.sparse.linalg import splu
 from scipy.sparse.linalg import minres, cg
 from scipy.linalg import cho_factor, cho_solve
+import time
+from loguru import logger
+from .misc import solve
+from .misc import logdet
+from .misc import inv
+
 
 class GPMarginalDensity:
     def __init__(self,
@@ -11,43 +17,63 @@ def __init__(self,
                  online=False,
                  calc_inv=False,
                  info=False):
+
         self.data_obj = data_obj
         self.prior_obj = prior_obj
         self.likelihood_obj = likelihood_obj
         self.calc_inv = calc_inv
+        self.online = online
         self.info = info
         self.K = prior_obj.K
         self.V = likelihood_obj.V
         self.y_data = data_obj.y_data
         self.y_mean = data_obj.y_data - prior_obj.prior_mean_vector
-        if online: self.calc_inv = True
+        if self.online: self.calc_inv = True
 
-        self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, self.KVinv = self._compute_GPpriorV()
+        self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, self.KVinv, mean = \
+            self.compute_GPpriorV(self.prior_obj.hyperparameters)
 
-    def update(self):
-        self.K = prior_obj.K
-        self.V = likelihood_obj.V
-        self.y_data = data_obj.y_data
-        self.y_mean = data_obj.y_data - prior_obj.prior_mean_vector
-        self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, self.KVinv = self._update_GPpriorV()
+    def update_data(self):
+        """Update the marginal PDF when the data has changed in data likelihood or prior objects"""
+        self.K = self.prior_obj.K
+        self.V = self.likelihood_obj.V
+        self.y_data = self.data_obj.y_data
+        self.y_mean = self.data_obj.y_data - self.prior_obj.prior_mean_vector
+        self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, self.KVinv = self.update_GPpriorV()
+
+    def update_hyperparameters(self):
+        """Update the marginal PDF when if hyperparameters have changed"""
+        self.K = self.prior_obj.K
+        self.V = self.likelihood_obj.V
+
+        self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, self.KVinv, self.prior_obj.prior_mean_vector = \
+            self.compute_GPpriorV(self.prior_obj.hyperparameters)
 
-    def _compute_GPpriorV(self):
+        self.y_mean = self.y_data - self.prior_obj.prior_mean_vector
+
+    def compute_GPpriorV(self, hyperparameters, calc_inv=None):
+        """Recomputed the prior mean for new hyperparameters (e.g. during training)"""
         # get K + V
-        K = self.prior_obj.K
-        V = self.likelihood_obj.V
+        if calc_inv is None: calc_inv = self.calc_inv
+        K = self.prior_obj.compute_kernel(self.data_obj.x_data, self.data_obj.x_data, hyperparameters=hyperparameters)
+        m = self.prior_obj.compute_mean(self.data_obj.x_data, hyperparameters=hyperparameters)
+        V = self.likelihood_obj.calculate_V(hyperparameters)
+        y_mean = self.data_obj.y_data - m
         # check if shapes are correct
         assert K.shape == V.shape
         KV = K + V
         # get Kinv/KVinvY, LU, Chol, logdet(KV)
-        KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(self.y_mean, KV)
-        return KV, KVinvY, KVlogdet, factorization_obj, KVinv
+        KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(y_mean, KV, calc_inv)
+        return KV, KVinvY, KVlogdet, factorization_obj, KVinv, m
 
     ##################################################################################
 
-    def _update_GPpriorV(self):
+    def update_GPpriorV(self):
+        """This updates the prior after new data was communicated"""
         # get K
         K = self.prior_obj.K
         V = self.likelihood_obj.V
+
         # check if shapes are correct
         assert K.shape == V.shape
 
@@ -59,17 +85,17 @@ def _update_GPpriorV(self):
             KVinvY, KVlogdet, factorization_obj, KVinv = self._update_gp_linalg(
                 self.y_mean, KV)
         else:
-            KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(self.y_mean, KV)
+            KVinvY, KVlogdet, factorization_obj, KVinv = self._compute_gp_linalg(self.y_mean, KV, self.calc_inv)
         return KV, KVinvY, KVlogdet, factorization_obj, KVinv
 
     ##################################################################################
-    def _compute_gp_linalg(self, vec, KV):
-        if self.calc_inv:
-            KVinv = self._inv(KV)
+    def _compute_gp_linalg(self, vec, KV, calc_inv):
+        if calc_inv:
+            KVinv = inv(KV)
             factorization_obj = ("Inv", None)
             KVinvY = KVinv @ vec
-            KVlogdet = self._logdet(KV)
-        elif not self.calc_inv:
+            KVlogdet = logdet(KV)
+        elif not calc_inv:
             KVinv = None
             KVinvY, KVlogdet, factorization_obj = self._Chol(KV, vec)
         else: raise Exception("calc_inv unspecified")
@@ -79,13 +105,13 @@ def _update_gp_linalg(self, vec, KV):
         size_KVinv = np.len(self.KVinv)
         kk = KV[size_KVinv:, size_KVinv:]
         k = KV[size_KVinv:, 0:size_KVinv]
-        X = self._inv(kk - C @ self.KVinv @ B)
+        X = 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)
+        KVlogdet = self.KVlogdet + logdet(kk - k.T @ self.KVinv @ k)
         return KVinvY, KVlogdet, factorization_obj, KVinv
 
     def _LU(self, KV, vec):
@@ -100,6 +126,7 @@ def _LU(self, KV, vec):
         return KVinvY, KVlogdet, factorization_obj
 
     def _Chol(self, KV, vec):
+        st = time.time()
         if self.info: print("Dense Cholesky in progress ...")
         c, l = cho_factor(KV)
         factorization_obj = ("Chol", c, l)
@@ -126,12 +153,8 @@ def log_likelihood(self, hyperparameters=None):
         log marginal likelihood of the data : float
         """
         if hyperparameters is None:
-            K, KV, KVinvY, KVlogdet, FO, KVinv, mean, cov = \
-                (self.K, self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj,
-                 self.KVinv, self.prior_mean_vec, self.V)
-        else:
-            K, KV, KVinvY, KVlogdet, FO, KVinv, mean, cov = \
-                self._compute_GPpriorV(self.x_data, self.y_data, hyperparameters, calc_inv=False)
+            KVinvY, KVlogdet, mean = self.KVinvY, self.KVlogdet, self.prior_obj.prior_mean_vector
+        else: KV, KVinvY, KVlogdet, factorization_obj, KVinv, mean = self.compute_GPpriorV(hyperparameters)
         n = len(self.y_data)
         return -(0.5 * ((self.y_data - mean).T @ KVinvY)) - (0.5 * KVlogdet) - (0.5 * n * np.log(2.0 * np.pi))
 
@@ -169,48 +192,44 @@ def neg_log_likelihood_gradient(self, hyperparameters=None):
         """
         logger.debug("log-likelihood gradient is being evaluated...")
         if hyperparameters is None:
-            K, KV, KVinvY, KVlogdet, FO, KVinv, mean, cov = \
-                self.K, self.KV, self.KVinvY, self.KVlogdet, self.factorization_obj, \
-                    self.KVinv, self.prior_mean_vec, self.V
-        else:
-            K, KV, KVinvY, KVlogdet, FO, KVinv, mean, cov = \
-                self._compute_GPpriorV(self.x_data, self.y_data, hyperparameters, calc_inv=False)
+            KVinvY, KVlogdet, mean = self.KVinvY, self.KVlogdet, self.prior_obj.prior_mean_vector
+        else: KV, KVinvY, KVlogdet, factorization_obj, KVinv, mean = self.compute_GPpriorV(hyperparameters)
 
         b = KVinvY
         y = self.y_data - mean
-        if self.ram_economy is False:
+        if self.prior_obj.ram_economy is False:
             try:
-                dK_dH = self.dk_dh(self.x_data, self.x_data, hyperparameters, self) + self.noise_function_grad(
-                    self.x_data, hyperparameters, self)
+                dK_dH = self.prior_obj.dk_dh(self.data_obj.x_data, self.data_obj.x_data, hyperparameters, self) + \
+                        self.likelihood_obj.noise_function_grad(self.data_obj.x_data, hyperparameters, self)
             except Exception as e:
                 raise Exception(
                     "The gradient evaluation dK/dh + dNoise/dh was not successful. \n \
                     That normally means the combination of ram_economy and definition \
                     of the gradient function is wrong. ",
                     str(e))
             KV = np.array([KV, ] * len(hyperparameters))
-            a = self._solve(KV, dK_dH)
+            a = solve(KV, dK_dH)
         bbT = np.outer(b, b.T)
         dL_dH = np.zeros((len(hyperparameters)))
         dL_dHm = np.zeros((len(hyperparameters)))
-        dm_dh = self.dm_dh(self.x_data, hyperparameters, self)
+        dm_dh = self.prior_obj.dm_dh(self.data_obj.x_data, hyperparameters, self)
 
         for i in range(len(hyperparameters)):
             dL_dHm[i] = -dm_dh[i].T @ b
-            if self.ram_economy is False:
+            if self.prior_obj.ram_economy is False:
                 matr = a[i]
             else:
                 try:
-                    dK_dH = self.dk_dh(self.x_data, self.x_data, i, hyperparameters, self) + self.noise_function_grad(
-                        self.x_data, i, hyperparameters, self)
+                    dK_dH = self.prior_obj.dk_dh(self.data_obj.x_data, self.data_obj.x_data, i, hyperparameters, self)+\
+                            self.noise_function_grad(self.data_obj.x_data, i, hyperparameters, self)
                 except:
                     raise Exception(
                         "The gradient evaluation dK/dh + dNoise/dh was not successful. \n \
                         That normally means the combination of ram_economy and definition of \
                         the gradient function is wrong.")
-                matr = np.linalg.solve(KV, dK_dH)
+                matr = solve(KV, dK_dH)
             if dL_dHm[i] == 0.0:
-                if self.ram_economy is False:
+                if self.prior_obj.ram_economy is False:
                     mtrace = np.einsum('ij,ji->', bbT, dK_dH[i])
                 else:
                     mtrace = np.einsum('ij,ji->', bbT, dK_dH)