16945523497

Committed 13 Aug 2025 06:07PM UTC coverage: 97.987% (+6.4%) from 91.538%

Build # 16945523497

Build Type

Pull #12

github

Committed by

web-flow

Commit Message

Merge 84ca55fbe into 4a2aa04a5

Pull Request Pull Request #12: Replacing pytorch with tensorflow for using native L-BFGS with autograd

Run Details

303 of 309 new or added lines in 8 files covered. (98.06%)

1 existing line in 1 file now uncovered.

438 of 447 relevant lines covered (97.99%)

3.91 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.0

/lcgp/lcgp.py

import tensorflow as tf
import tensorflow_probability as tfp
import gpflow
from .covmat import Matern32
import numpy as np

# for Python 3.9 inclusion
from typing import Optional

# Display only code-breaking errors
tf.get_logger().setLevel('ERROR')
# Set default float type to float64
tf.keras.backend.set_floatx('float64')


class LCGP(gpflow.Module):
    def __init__(self,
                 y: Optional[np.ndarray] = tf.Tensor,
                 x: Optional[np.ndarray] = tf.Tensor,
                 q: int = None,
                 var_threshold: float = None,
                 diag_error_structure: list = None,
                 parameter_clamp_flag: bool = False,
                 robust_mean: bool = True,
                 penalty_const: dict = None,
                 submethod: str = 'full',
                 verbose: bool = False):
        """
        Constructor for LCGP class.
        """
        super().__init__()
        self.verbose = verbose
        self.robust_mean = robust_mean
        self.x = self._verify_data_types(x)
        self.y = self._verify_data_types(y)

        self.method = 'LCGP'
        self.submethod = submethod
        self.submethod_loss_map = {'full': self.neglpost,
                                   # 'elbo': self.negelbo,
                                   # 'proflik': self.negproflik
                                   }
        self.submethod_predict_map = {'full': self.predict_full,
                                      # 'elbo': self.predict_elbo,
                                      # 'proflik': self.predict_proflik
                                      }

        self.parameter_clamp_flag = parameter_clamp_flag

        if (q is not None) and (var_threshold is not None):
            raise ValueError('Include only q or var_threshold but not both.')
        self.q = q
        self.var_threshold = var_threshold

        # standardize x to unit hypercube
        self.x, self.x_min, self.x_max, self.x_orig, self.xnorm = \
            self.init_standard_x(self.x)
        # standardize y
        self.y, self.ymean, self.ystd, self.y_orig = self.init_standard_y(self.y)

        # placeholders for variables
        self.n, self.d, self.p = 0., 0., 0.
        # verify that input and output dimensions match
        # sets n, d, and p
        self.verify_dim(self.y, self.x)

        # reset q if none is provided
        self.g, self.phi, self.diag_D, self.q = \
            self.init_phi(var_threshold=var_threshold)

        if diag_error_structure is None:
            self.diag_error_structure = [1] * int(self.p)
        else:
            self.diag_error_structure = diag_error_structure

        self.verify_error_structure(self.diag_error_structure, self.y)

        # Initialize parameters
        self.lLmb = gpflow.Parameter(tf.ones([self.q, self.x.shape[1]], dtype=tf.float64),
                                     name='Latent GP log-scale',
                                     transform=tfp.bijectors.SoftClip(
                                         low=tf.constant(1e-6, dtype=tf.float64),
                                         high=tf.constant(1e4, dtype=tf.float64)
                                     ), dtype=tf.float64)
        self.lLmb0 = gpflow.Parameter(tf.ones([self.q], dtype=tf.float64),
                                      name='Latent GP log-lengthscale',
                                      transform=tfp.bijectors.SoftClip(
                                          low=tf.constant(1e-4, dtype=tf.float64),
                                          high=tf.constant(1e4, dtype=tf.float64)
                                      ), dtype=tf.float64)
        self.lsigma2s = gpflow.Parameter(tf.ones([len(self.diag_error_structure)], dtype=tf.float64),
                                         name='Diagonal error log-variance') #, transform=tfp.bijectors.Exp())
        self.lnugGPs = gpflow.Parameter(tf.ones([self.q], dtype=tf.float64) * 1e-6,
                                        name='Latent GP nugget scale',
                                        transform=tfp.bijectors.SoftClip(
                                            low=tf.math.exp(tf.constant(-16, dtype=tf.float64)),
                                            high=tf.math.exp(tf.constant(-2, dtype=tf.float64))
                                        ), dtype=tf.float64)

        if penalty_const is None:
            pc = {'lLmb': 40, 'lLmb0': 5}
        else:
            pc = penalty_const
            for k, v in pc.items():
                assert v >= 0, 'penalty constant should be nonnegative.'
        self.penalty_const = pc

        self.init_params()

        # placeholders for predictive quantities
        self.CinvMs = tf.fill([self.q, self.n], tf.constant(float('nan'), dtype=tf.float64))
        self.Ths = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))
        self.Th_hats = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))
        self.Cinvhs = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))

    @staticmethod
    def init_standard_x(x):
        """
        Standardizes training inputs and collects summary information.
        """
        x_max = tf.reduce_max(x, axis=0)
        x_min = tf.reduce_min(x, axis=0)
        xs = (x - x_min) / (x_max - x_min)

        xnorm = tf.zeros(x.shape[1], dtype=tf.float64)
        for j in range(x.shape[1]):
            xdist = tf.abs((tf.reshape(x[:, j], (-1, 1)) - x[:, j]))

            positive_xdist = tf.boolean_mask(xdist, xdist > 0)
            mean_val = tf.reduce_mean(positive_xdist)

            xnorm = tf.tensor_scatter_nd_update(xnorm, [[j]], [mean_val])
        return xs, x_min, x_max, x, xnorm

    def init_standard_y(self, y):
        """
        Standardizes outputs and collects summary information.
        """
        if self.robust_mean:
            ycenter = tfp.stats.percentile(y, 50.0, axis=1, keepdims=True)
            yspread = tfp.stats.percentile(tf.abs(y - ycenter), 50.0, axis=1, keepdims=True)
        else:
            ycenter = tf.reduce_mean(y, axis=1, keepdims=True)
            yspread = tf.math.reduce_std(y, axis=1, keepdims=True)

        ys = (y - ycenter) / yspread
        return ys, ycenter, yspread, y

    def __repr__(self):
        params = gpflow.utilities.tabulate_module_summary(self)
        desc = 'LCGP(\n' \
               '\tsubmethod:\t{:s}\n' \
               '\toutput dimension:\t{:d}\n' \
               '\tnumber of latent components:\t{:d}\n' \
               '\tparameter_clamping:\t{:s}\n' \
               '\trobust_standardization:\t{:s}\n' \
               '\tdiagonal_error structure:\t{:s}\n' \
               '\tparameters:\t\n{}\n)'.format(self.submethod, self.p,
                                         self.q, str(self.parameter_clamp_flag),
                                         str(self.robust_mean),
                                         str(self.diag_error_structure),
                                         params)
        return desc

    def init_phi(self, var_threshold: float = None):
        """
        Initialization of orthogonal basis, computed with singular value decomposition.
        """
        y, q = self.y, self.q
        n, p = self.n, self.p

        singvals, left_u, _ = tf.linalg.svd(y, full_matrices=False)

        if (q is None) and (var_threshold is None):
            q = p
        elif (q is None) and (var_threshold is not None):
            cumvar = tf.cumsum(singvals ** 2) / tf.reduce_sum(singvals ** 2)
            q = int(tf.argmax(cumvar > var_threshold) + 1)

        assert left_u.shape[1] == min(n, p)
        singvals = singvals[:q]

        # Compute phi and diag_D
        phi = left_u[:, :q] * tf.sqrt(tf.cast(n, tf.float64)) / singvals
        diag_D = tf.reduce_sum(phi ** 2, axis=0)

        g = tf.matmul(phi, y, transpose_a=True)
        return g, phi, diag_D, q

    def init_params(self):
        """
        Initializes parameters for LCGP.
        """
        x = self.x
        d = self.d

        llmb = np.exp(0.5 * np.log(d) + np.log(np.std(x, axis=0)))
        lLmb = np.tile(llmb, self.q).reshape((self.q, d))
        lLmb0 = np.ones(self.q, dtype=np.float64)
        lnugGPs = np.exp(-10.) * np.ones(self.q, dtype=np.float64)

        err_struct = self.diag_error_structure
        lsigma2_diag = np.zeros(len(err_struct), dtype=np.float64)
        col = 0
        for k in range(len(err_struct)):
            lsigma2_diag[k] = np.log(np.var(self.y[col:(col + err_struct[k])]))
            col += err_struct[k]

        self.lLmb.assign(lLmb)
        self.lLmb0.assign(lLmb0)
        self.lnugGPs.assign(lnugGPs)
        self.lsigma2s.assign(lsigma2_diag)
        return

    def verify_dim(self, y, x):
        """
        Verifies if input and output dimensions match. Sets class variables for
        dimensions. Throws error if the dimensions do not match.
        """
        p, ny = tf.shape(y)[0], tf.shape(y)[1]
        nx, d = tf.shape(x)[0], tf.shape(x)[1]

        assert ny == nx, 'Number of inputs (x) differs from number of outputs (y), y.shape[1] != x.shape[0]'

        self.n = tf.constant(nx, tf.int32)
        self.d = tf.constant(d, tf.int32)
        self.p = tf.constant(p, tf.int32)
        return

    def tx_x(self, xs):
        """
        Reverts standardization of inputs.
        """
        return xs * (self.x_max - self.x_min) + self.x_min

    def tx_y(self, ys):
        """
        Reverts output standardization.
        """
        return ys * self.ystd + self.ymean

    def fit(self, verbose=False):
        opt = gpflow.optimizers.Scipy()
        opt.minimize(self.loss, self.trainable_variables)
        return

    def loss(self):
        """
        Computes the loss based on the submethod.
        """
        if self.submethod == 'full':
            return self.neglpost()
        # elif self.submethod == 'elbo':
        #     return self.negelbo()
        # elif self.submethod == 'proflik':
        #     return self.negproflik()
        else:
            raise ValueError("Invalid submethod. Choices are 'full', 'elbo', or 'proflik'.")

    @tf.function
    def neglpost(self):
        lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()
        x = self.x
        y = self.y

        pc = self.penalty_const

        n = self.n
        q = self.q
        D = self.diag_D
        phi = self.phi
        psi_c = tf.transpose(phi) / tf.sqrt(tf.exp(lsigma2s))

        nlp = tf.constant(0., dtype=tf.float64)

        for k in range(q):
            Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])

            Wk, Uk = tf.linalg.eigh(Ck)

            Qk = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))
            Pk = tf.matmul(tf.expand_dims(psi_c[k], axis=1), tf.expand_dims(psi_c[k], axis=0))

            yQk = tf.matmul(y, Qk)
            yPk = tf.matmul(tf.transpose(y), tf.transpose(Pk))

            nlp += (0.5 * tf.reduce_sum(tf.math.log(1 + D[k] * Wk)))
            nlp += -(0.5 * tf.reduce_sum(yQk * tf.transpose(yPk)))

        nlp += (n / 2 * tf.reduce_sum(lsigma2s))
        nlp += (0.5 * tf.reduce_sum(tf.square(tf.transpose(y) / tf.sqrt(tf.exp(lsigma2s)))))

        # Regularization
        nlp += (pc['lLmb'] * tf.reduce_sum(tf.square(tf.math.log(lLmb))) +
                pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0.unconstrained_variable)))
        nlp += (-tf.reduce_sum(tf.math.log(tf.math.log(lnugGPs) + 100)))
        # nlp += (tf.reduce_sum(tf.math.log(lnugGPs - 100)))
        nlp /= tf.cast(n, tf.float64)
        return nlp

    # def negelbo(self):
    #     n = self.n
    #     x = self.x
    #     y = self.y
    #     pc = self.penalty_const
    #
    #     lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()
    #     B = tf.matmul(tf.transpose(y / tf.sqrt(tf.exp(lsigma2s))), self.phi)
    #     D = self.diag_D
    #     phi = self.phi
    #
    #     psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))
    #
    #     M = tf.zeros([self.q, n], dtype=tf.float64)
    #
    #     negelbo = tf.constant(0., dtype=tf.float64)
    #     for k in range(self.q):
    #         Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])
    #
    #         Wk, Uk = tf.linalg.eigh(Ck)
    #         dkInpCkinv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / Wk), tf.transpose(Uk))) + \
    #                      D[k] * tf.eye(n, dtype=tf.float64)
    #
    #         # (dk * I + Ck^{-1})^{-1}
    #         dkInpCkinv_inv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))
    #         Mk = tf.linalg.matvec(dkInpCkinv_inv, tf.transpose(B)[k])
    #         Vk = 1 / tf.linalg.diag_part(dkInpCkinv)
    #
    #         CkinvhMk = tf.linalg.matvec(tf.matmul(tf.matmul(Uk, tf.linalg.diag(1 / tf.sqrt(Wk))), tf.transpose(Uk)),
    #                                     Mk)
    #
    #         M = tf.tensor_scatter_nd_update(M, [[k]], tf.expand_dims(Mk, axis=0))
    #
    #         negelbo += 0.5 * tf.reduce_sum(tf.math.log(Wk))
    #         negelbo += 0.5 * tf.reduce_sum(tf.square(CkinvhMk))
    #         negelbo -= 0.5 * tf.reduce_sum(tf.math.log(Vk))
    #         negelbo += 0.5 * tf.reduce_sum(
    #             Vk * D[k] * tf.linalg.diag_part(tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / Wk), tf.transpose(Uk)))))
    #
    #     resid = (tf.transpose(y) - tf.matmul(tf.transpose(M), psi)) / tf.sqrt(tf.exp(lsigma2s))
    #
    #     negelbo += 0.5 * tf.reduce_sum(tf.square(resid))
    #     negelbo += n / 2 * tf.reduce_sum(lsigma2s)
    #
    #     # Regularization
    #     negelbo += pc['lLmb'] * tf.reduce_sum(tf.square(lLmb)) + \
    #                pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0))
    #     negelbo += -tf.reduce_sum(tf.math.log(lnugGPs + 100))
    #
    #     negelbo /= tf.cast(n, tf.float64)
    #
    #     return negelbo
    #
    # def negproflik(self):
    #     lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()
    #     x = self.x
    #     y = self.y
    #
    #     pc = self.penalty_const
    #
    #     n = self.n
    #     q = self.q
    #     D = self.diag_D
    #     phi = self.phi
    #     psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))
    #
    #     B = tf.matmul(tf.transpose(y / tf.sqrt(tf.exp(lsigma2s))), self.phi)
    #     G = tf.zeros([self.q, n], dtype=tf.float64)
    #
    #     negproflik = tf.constant(0., dtype=tf.float64)
    #
    #     for k in range(q):
    #         Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])
    #         Wk, Uk = tf.linalg.eigh(Ck)
    #
    #         dkInpCkinv_inv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))
    #         Gk = tf.matmul(dkInpCkinv_inv, tf.transpose(B)[k])
    #
    #         CkinvhGk = tf.matmul(tf.matmul(Uk, tf.linalg.diag(1 / Wk)), tf.transpose(Uk)) @ Gk
    #
    #         G = tf.tensor_scatter_nd_update(G, [[k]], tf.expand_dims(Gk, axis=0))
    #
    #         negproflik += 0.5 * tf.reduce_sum(tf.math.log(Wk))
    #         negproflik += 0.5 * tf.reduce_sum(tf.square(CkinvhGk))
    #
    #     resid = (tf.transpose(y) - tf.matmul(tf.transpose(G), psi)) / tf.sqrt(tf.exp(lsigma2s))
    #
    #     negproflik += 0.5 * tf.reduce_sum(tf.square(resid))
    #     negproflik += n / 2 * tf.reduce_sum(lsigma2s)
    #
    #     negproflik += pc['lLmb'] * tf.reduce_sum(tf.square(lLmb)) + \
    #                   pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0))
    #     negproflik += -tf.reduce_sum(tf.math.log(lnugGPs + 100))
    #
    #     negproflik /= tf.cast(n, tf.float64)
    #     return negproflik


    def predict(self, x0, return_fullcov=False):
        """
        Returns predictive quantities at new input `x0`.  Both outputs are of
        size (number of new input, output dimension).
        :param x0: New input of size (number of new input, dimension of input).
        :param return_fullcov: Returns (predictive mean, predictive variance,
        variance for the true mean, full predictive covariance) if True.  Otherwise,
        only return the first three quantities.
        """
        x0 = self._verify_data_types(x0)
        submethod = self.submethod
        predict_map = self.submethod_predict_map
        try:
            predict_call = predict_map[submethod]
        except KeyError as e:
            print(e)
            # print('Invalid submethod.  Choices are \'full\', \'elbo\', or \'proflik\'.')
            raise KeyError('Invalid submethod.  Choices are \'full\'.')
        return tf.stop_gradient(predict_call(x0=x0, return_fullcov=return_fullcov))


    def compute_aux_predictive_quantities(self):
        """
        Compute auxiliary quantities for predictions using full posterior approach.
        """
        x = self.x
        lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()

        D = self.diag_D
        # B := Y @ Sigma^{-1/2} @ Phi
        B = tf.matmul(tf.transpose(self.y) / tf.sqrt(tf.exp(lsigma2s)), self.phi)

        CinvM = tf.zeros([self.q, self.n], dtype=tf.float64)
        Th = tf.zeros([self.q, self.n, self.n], dtype=tf.float64)

        for k in range(self.q):
            Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])

            Wk, Uk = tf.linalg.eigh(Ck)

            # (I + D_k * C_k)^{-1}
            IpdkCkinv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1.0 / (1.0 + D[k] * Wk)), tf.transpose(Uk)))

            CkinvMk = tf.linalg.matvec(IpdkCkinv, tf.transpose(B)[k])
            Thk = tf.matmul(Uk, tf.matmul(tf.linalg.diag(tf.sqrt((D[k] * Wk ** 2) / (Wk ** 2 + D[k] * Wk ** 3))),
                                          tf.transpose(Uk)))

            CinvM = tf.tensor_scatter_nd_update(CinvM, [[k]], tf.expand_dims(CkinvMk, axis=0))
            Th = tf.tensor_scatter_nd_update(Th, [[k]], tf.expand_dims(Thk, axis=0))

        self.CinvMs = CinvM
        self.Ths = Th

    def predict_full(self, x0, return_fullcov=False):
        """
        Returns predictions using full posterior approach.
        """
        if tf.reduce_any(tf.math.is_nan(self.CinvMs)) or tf.reduce_any(tf.math.is_nan(self.Ths)):
            self.compute_aux_predictive_quantities()

        x = self.x
        lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()

        phi = self.phi

        CinvM = self.CinvMs
        Th = self.Ths

        x0 = self._verify_data_types(x0)
        x0 = (x0 - self.x_min) / (self.x_max - self.x_min)  # Standardize x0
        n0 = tf.shape(x0)[0]

        ghat = tf.zeros([self.q, n0], dtype=tf.float64)
        gvar = tf.zeros([self.q, n0], dtype=tf.float64)
        for k in range(self.q):
            c00k = Matern32(x0, x0, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k],
                            diag_only=True)  # Diagonal-only covariance
            c0k = Matern32(x0, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k],
                           diag_only=False)

            ghat_k = tf.linalg.matvec(c0k, CinvM[k])
            gvar_k = c00k - tf.reduce_sum(tf.square(tf.matmul(c0k, Th[k])), axis=1)

            ghat = tf.tensor_scatter_nd_update(ghat, [[k]], [ghat_k])
            gvar = tf.tensor_scatter_nd_update(gvar, [[k]], [gvar_k])

        self.ghat = ghat
        self.gvar = gvar

        psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))

        predmean = tf.matmul(psi, ghat, transpose_a=True)
        confvar = tf.matmul(tf.transpose(gvar), tf.square(psi))
        predvar = confvar + tf.exp(lsigma2s)

        ypred = self.tx_y(predmean)
        yconfvar = tf.transpose(confvar) * tf.square(self.ystd)
        ypredvar = tf.transpose(predvar) * tf.square(self.ystd)

        if return_fullcov:
            CH = tf.sqrt(gvar)[..., tf.newaxis] * psi[tf.newaxis, ...]
            yfullpredcov = (tf.einsum('nij,njk->nik', CH,
                                     tf.transpose(CH, perm=[0, 2, 1])) +
                            tf.linalg.diag(tf.exp(lsigma2s)))
            yfullpredcov *= tf.square(self.ystd)
            return ypred, ypredvar, yconfvar, yfullpredcov

        return ypred, ypredvar, yconfvar

    @staticmethod
    def _verify_data_types(t):
        """
        Verify if inputs are TensorFlow tensors, if not, cast into tensors.
        Verify if inputs are at least 2-dimensional, if not, expand dimensions to 2.
        """
        if not isinstance(t, tf.Tensor):
            t = tf.convert_to_tensor(t, dtype=tf.float64)
        if t.ndim < 2:
            t = tf.expand_dims(t, axis=1)
        return t

    # def predict_elbo(self, x0, return_fullcov=False):
    #     pass
    #
    # def predict_proflik(self, x0, return_fullcov=False):
    #     pass


    @staticmethod
    def verify_error_structure(diag_error_structure, y):
        """
        Verifies if diagonal error structure input, if any, is valid.
        """
        assert sum(diag_error_structure) == y.shape[0], \
            'Sum of error_structure should' \
            ' equal the output dimension.'

    def get_param(self):
        """
        Returns the parameters for LCGP instance.
        """
        # if self.parameter_clamp_flag:
        #     lLmb, lLmb0, lsigma2s, lnugGPs = \
        #         self.parameter_clamp(lLmb=self.lLmb, lLmb0=self.lLmb0,
        #                              lsigma2s=self.lsigma2s, lnugs=self.lnugGPs)
        # else:
        lLmb, lLmb0, lsigma2s, lnugGPs = \
            self.lLmb, self.lLmb0, self.lsigma2s, self.lnugGPs

        built_lsigma2s = tf.zeros(self.p, dtype=tf.float64)
        err_struct = self.diag_error_structure
        col = 0
        for k in range(len(err_struct)):
            built_lsigma2s = tf.tensor_scatter_nd_update(
                built_lsigma2s,
                tf.range(col, col + err_struct[k])[:, tf.newaxis],
                tf.fill([err_struct[k]], lsigma2s[k])
            )
            col += err_struct[k]

        return lLmb, lLmb0, built_lsigma2s, lnugGPs

1	import tensorflow as tf	4✔
2	import tensorflow_probability as tfp	4✔
3	import gpflow	4✔
4	from .covmat import Matern32	4✔
5	import numpy as np	4✔
6
7	# for Python 3.9 inclusion
8	from typing import Optional	4✔
9
10	# Display only code-breaking errors
11	tf.get_logger().setLevel('ERROR')	4✔
12	# Set default float type to float64
13	tf.keras.backend.set_floatx('float64')	4✔
14
15
16	class LCGP(gpflow.Module):	4✔
17	def __init__(self,	4✔
18	y: Optional[np.ndarray] = tf.Tensor,
19	x: Optional[np.ndarray] = tf.Tensor,
20	q: int = None,
21	var_threshold: float = None,
22	diag_error_structure: list = None,
23	parameter_clamp_flag: bool = False,
24	robust_mean: bool = True,
25	penalty_const: dict = None,
26	submethod: str = 'full',
27	verbose: bool = False):
28	"""
29	Constructor for LCGP class.
30	"""
31	super().__init__()	4✔
32	self.verbose = verbose	4✔
33	self.robust_mean = robust_mean	4✔
34	self.x = self._verify_data_types(x)	4✔
35	self.y = self._verify_data_types(y)	4✔
36
37	self.method = 'LCGP'	4✔
38	self.submethod = submethod	4✔
39	self.submethod_loss_map = {'full': self.neglpost,	4✔
40	# 'elbo': self.negelbo,
41	# 'proflik': self.negproflik
42	}
43	self.submethod_predict_map = {'full': self.predict_full,	4✔
44	# 'elbo': self.predict_elbo,
45	# 'proflik': self.predict_proflik
46	}
47
48	self.parameter_clamp_flag = parameter_clamp_flag	4✔
49
50	if (q is not None) and (var_threshold is not None):	4✔
51	raise ValueError('Include only q or var_threshold but not both.')	4✔
52	self.q = q	4✔
53	self.var_threshold = var_threshold	4✔
54
55	# standardize x to unit hypercube
56	self.x, self.x_min, self.x_max, self.x_orig, self.xnorm = \	4✔
57	self.init_standard_x(self.x)
58	# standardize y
59	self.y, self.ymean, self.ystd, self.y_orig = self.init_standard_y(self.y)	4✔
60
61	# placeholders for variables
62	self.n, self.d, self.p = 0., 0., 0.	4✔
63	# verify that input and output dimensions match
64	# sets n, d, and p
65	self.verify_dim(self.y, self.x)	4✔
66
67	# reset q if none is provided
68	self.g, self.phi, self.diag_D, self.q = \	4✔
69	self.init_phi(var_threshold=var_threshold)
70
71	if diag_error_structure is None:	4✔
72	self.diag_error_structure = [1] * int(self.p)	4✔
73	else:
74	self.diag_error_structure = diag_error_structure	4✔
75
76	self.verify_error_structure(self.diag_error_structure, self.y)	4✔
77
78	# Initialize parameters
79	self.lLmb = gpflow.Parameter(tf.ones([self.q, self.x.shape[1]], dtype=tf.float64),	4✔
80	name='Latent GP log-scale',
81	transform=tfp.bijectors.SoftClip(
82	low=tf.constant(1e-6, dtype=tf.float64),
83	high=tf.constant(1e4, dtype=tf.float64)
84	), dtype=tf.float64)
85	self.lLmb0 = gpflow.Parameter(tf.ones([self.q], dtype=tf.float64),	4✔
86	name='Latent GP log-lengthscale',
87	transform=tfp.bijectors.SoftClip(
88	low=tf.constant(1e-4, dtype=tf.float64),
89	high=tf.constant(1e4, dtype=tf.float64)
90	), dtype=tf.float64)
91	self.lsigma2s = gpflow.Parameter(tf.ones([len(self.diag_error_structure)], dtype=tf.float64),	4✔
92	name='Diagonal error log-variance') #, transform=tfp.bijectors.Exp())
93	self.lnugGPs = gpflow.Parameter(tf.ones([self.q], dtype=tf.float64) * 1e-6,	4✔
94	name='Latent GP nugget scale',
95	transform=tfp.bijectors.SoftClip(
96	low=tf.math.exp(tf.constant(-16, dtype=tf.float64)),
97	high=tf.math.exp(tf.constant(-2, dtype=tf.float64))
98	), dtype=tf.float64)
99
100	if penalty_const is None:	4✔
101	pc = {'lLmb': 40, 'lLmb0': 5}	4✔
102	else:
103	pc = penalty_const	4✔
104	for k, v in pc.items():	4✔
105	assert v >= 0, 'penalty constant should be nonnegative.'	4✔
106	self.penalty_const = pc	4✔
107
108	self.init_params()	4✔
109
110	# placeholders for predictive quantities
111	self.CinvMs = tf.fill([self.q, self.n], tf.constant(float('nan'), dtype=tf.float64))	4✔
112	self.Ths = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))	4✔
113	self.Th_hats = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))	4✔
114	self.Cinvhs = tf.fill([self.q, self.n, self.n], tf.constant(float('nan'), dtype=tf.float64))	4✔
115
116	@staticmethod	4✔
117	def init_standard_x(x):	4✔
118	"""
119	Standardizes training inputs and collects summary information.
120	"""
121	x_max = tf.reduce_max(x, axis=0)	4✔
122	x_min = tf.reduce_min(x, axis=0)	4✔
123	xs = (x - x_min) / (x_max - x_min)	4✔
124
125	xnorm = tf.zeros(x.shape[1], dtype=tf.float64)	4✔
126	for j in range(x.shape[1]):	4✔
127	xdist = tf.abs((tf.reshape(x[:, j], (-1, 1)) - x[:, j]))	4✔
128
129	positive_xdist = tf.boolean_mask(xdist, xdist > 0)	4✔
130	mean_val = tf.reduce_mean(positive_xdist)	4✔
131
132	xnorm = tf.tensor_scatter_nd_update(xnorm, [[j]], [mean_val])	4✔
133	return xs, x_min, x_max, x, xnorm	4✔
134
135	def init_standard_y(self, y):	4✔
136	"""
137	Standardizes outputs and collects summary information.
138	"""
139	if self.robust_mean:	4✔
140	ycenter = tfp.stats.percentile(y, 50.0, axis=1, keepdims=True)	4✔
141	yspread = tfp.stats.percentile(tf.abs(y - ycenter), 50.0, axis=1, keepdims=True)	4✔
142	else:
143	ycenter = tf.reduce_mean(y, axis=1, keepdims=True)	4✔
144	yspread = tf.math.reduce_std(y, axis=1, keepdims=True)	4✔
145
146	ys = (y - ycenter) / yspread	4✔
147	return ys, ycenter, yspread, y	4✔
148
149	def __repr__(self):	4✔
150	params = gpflow.utilities.tabulate_module_summary(self)	4✔
151	desc = 'LCGP(\n' \	4✔
152	'\tsubmethod:\t{:s}\n' \
153	'\toutput dimension:\t{:d}\n' \
154	'\tnumber of latent components:\t{:d}\n' \
155	'\tparameter_clamping:\t{:s}\n' \
156	'\trobust_standardization:\t{:s}\n' \
157	'\tdiagonal_error structure:\t{:s}\n' \
158	'\tparameters:\t\n{}\n)'.format(self.submethod, self.p,
159	self.q, str(self.parameter_clamp_flag),
160	str(self.robust_mean),
161	str(self.diag_error_structure),
162	params)
163	return desc	4✔
164
165	def init_phi(self, var_threshold: float = None):	4✔
166	"""
167	Initialization of orthogonal basis, computed with singular value decomposition.
168	"""
169	y, q = self.y, self.q	4✔
170	n, p = self.n, self.p	4✔
171
172	singvals, left_u, _ = tf.linalg.svd(y, full_matrices=False)	4✔
173
174	if (q is None) and (var_threshold is None):	4✔
175	q = p	4✔
176	elif (q is None) and (var_threshold is not None):	4✔
177	cumvar = tf.cumsum(singvals 2) / tf.reduce_sum(singvals 2)	4✔
178	q = int(tf.argmax(cumvar > var_threshold) + 1)	4✔
179
180	assert left_u.shape[1] == min(n, p)	4✔
181	singvals = singvals[:q]	4✔
182
183	# Compute phi and diag_D
184	phi = left_u[:, :q] * tf.sqrt(tf.cast(n, tf.float64)) / singvals	4✔
185	diag_D = tf.reduce_sum(phi ** 2, axis=0)	4✔
186
187	g = tf.matmul(phi, y, transpose_a=True)	4✔
188	return g, phi, diag_D, q	4✔
189
190	def init_params(self):	4✔
191	"""
192	Initializes parameters for LCGP.
193	"""
194	x = self.x	4✔
195	d = self.d	4✔
196
197	llmb = np.exp(0.5 * np.log(d) + np.log(np.std(x, axis=0)))	4✔
198	lLmb = np.tile(llmb, self.q).reshape((self.q, d))	4✔
199	lLmb0 = np.ones(self.q, dtype=np.float64)	4✔
200	lnugGPs = np.exp(-10.) * np.ones(self.q, dtype=np.float64)	4✔
201
202	err_struct = self.diag_error_structure	4✔
203	lsigma2_diag = np.zeros(len(err_struct), dtype=np.float64)	4✔
204	col = 0	4✔
205	for k in range(len(err_struct)):	4✔
206	lsigma2_diag[k] = np.log(np.var(self.y[col:(col + err_struct[k])]))	4✔
207	col += err_struct[k]	4✔
208
209	self.lLmb.assign(lLmb)	4✔
210	self.lLmb0.assign(lLmb0)	4✔
211	self.lnugGPs.assign(lnugGPs)	4✔
212	self.lsigma2s.assign(lsigma2_diag)	4✔
213	return	4✔
214
215	def verify_dim(self, y, x):	4✔
216	"""
217	Verifies if input and output dimensions match. Sets class variables for
218	dimensions. Throws error if the dimensions do not match.
219	"""
220	p, ny = tf.shape(y)[0], tf.shape(y)[1]	4✔
221	nx, d = tf.shape(x)[0], tf.shape(x)[1]	4✔
222
223	assert ny == nx, 'Number of inputs (x) differs from number of outputs (y), y.shape[1] != x.shape[0]'	4✔
224
225	self.n = tf.constant(nx, tf.int32)	4✔
226	self.d = tf.constant(d, tf.int32)	4✔
227	self.p = tf.constant(p, tf.int32)	4✔
228	return	4✔
229
230	def tx_x(self, xs):	4✔
231	"""
232	Reverts standardization of inputs.
233	"""
234	return xs * (self.x_max - self.x_min) + self.x_min	4✔
235
236	def tx_y(self, ys):	4✔
237	"""
238	Reverts output standardization.
239	"""
240	return ys * self.ystd + self.ymean	4✔
241
242	def fit(self, verbose=False):	4✔
243	opt = gpflow.optimizers.Scipy()	4✔
244	opt.minimize(self.loss, self.trainable_variables)	4✔
245	return	4✔
246
247	def loss(self):	4✔
248	"""
249	Computes the loss based on the submethod.
250	"""
251	if self.submethod == 'full':	4✔
252	return self.neglpost()	4✔
253	# elif self.submethod == 'elbo':
254	# return self.negelbo()
255	# elif self.submethod == 'proflik':
256	# return self.negproflik()
257	else:
258	raise ValueError("Invalid submethod. Choices are 'full', 'elbo', or 'proflik'.")	4✔
259
260	@tf.function	4✔
261	def neglpost(self):	4✔
262	lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()	4✔
263	x = self.x	4✔
264	y = self.y	4✔
265
266	pc = self.penalty_const	4✔
267
268	n = self.n	4✔
269	q = self.q	4✔
270	D = self.diag_D	4✔
271	phi = self.phi	4✔
272	psi_c = tf.transpose(phi) / tf.sqrt(tf.exp(lsigma2s))	4✔
273
274	nlp = tf.constant(0., dtype=tf.float64)	4✔
275
276	for k in range(q):	4✔
277	Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])	4✔
278
279	Wk, Uk = tf.linalg.eigh(Ck)	4✔
280
281	Qk = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))	4✔
282	Pk = tf.matmul(tf.expand_dims(psi_c[k], axis=1), tf.expand_dims(psi_c[k], axis=0))	4✔
283
284	yQk = tf.matmul(y, Qk)	4✔
285	yPk = tf.matmul(tf.transpose(y), tf.transpose(Pk))	4✔
286
287	nlp += (0.5 * tf.reduce_sum(tf.math.log(1 + D[k] * Wk)))	4✔
288	nlp += -(0.5 * tf.reduce_sum(yQk * tf.transpose(yPk)))	4✔
289
290	nlp += (n / 2 * tf.reduce_sum(lsigma2s))	4✔
291	nlp += (0.5 * tf.reduce_sum(tf.square(tf.transpose(y) / tf.sqrt(tf.exp(lsigma2s)))))	4✔
292
293	# Regularization
294	nlp += (pc['lLmb'] * tf.reduce_sum(tf.square(tf.math.log(lLmb))) +	4✔
295	pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0.unconstrained_variable)))
296	nlp += (-tf.reduce_sum(tf.math.log(tf.math.log(lnugGPs) + 100)))	4✔
297	# nlp += (tf.reduce_sum(tf.math.log(lnugGPs - 100)))
298	nlp /= tf.cast(n, tf.float64)	4✔
299	return nlp	4✔
300
301	# def negelbo(self):
302	# n = self.n
303	# x = self.x
304	# y = self.y
305	# pc = self.penalty_const
306	#
307	# lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()
308	# B = tf.matmul(tf.transpose(y / tf.sqrt(tf.exp(lsigma2s))), self.phi)
309	# D = self.diag_D
310	# phi = self.phi
311	#
312	# psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))
313	#
314	# M = tf.zeros([self.q, n], dtype=tf.float64)
315	#
316	# negelbo = tf.constant(0., dtype=tf.float64)
317	# for k in range(self.q):
318	# Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])
319	#
320	# Wk, Uk = tf.linalg.eigh(Ck)
321	# dkInpCkinv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / Wk), tf.transpose(Uk))) + \
322	# D[k] * tf.eye(n, dtype=tf.float64)
323	#
324	# # (dk * I + Ck^{-1})^{-1}
325	# dkInpCkinv_inv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))
326	# Mk = tf.linalg.matvec(dkInpCkinv_inv, tf.transpose(B)[k])
327	# Vk = 1 / tf.linalg.diag_part(dkInpCkinv)
328	#
329	# CkinvhMk = tf.linalg.matvec(tf.matmul(tf.matmul(Uk, tf.linalg.diag(1 / tf.sqrt(Wk))), tf.transpose(Uk)),
330	# Mk)
331	#
332	# M = tf.tensor_scatter_nd_update(M, [[k]], tf.expand_dims(Mk, axis=0))
333	#
334	# negelbo += 0.5 * tf.reduce_sum(tf.math.log(Wk))
335	# negelbo += 0.5 * tf.reduce_sum(tf.square(CkinvhMk))
336	# negelbo -= 0.5 * tf.reduce_sum(tf.math.log(Vk))
337	# negelbo += 0.5 * tf.reduce_sum(
338	# Vk * D[k] * tf.linalg.diag_part(tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / Wk), tf.transpose(Uk)))))
339	#
340	# resid = (tf.transpose(y) - tf.matmul(tf.transpose(M), psi)) / tf.sqrt(tf.exp(lsigma2s))
341	#
342	# negelbo += 0.5 * tf.reduce_sum(tf.square(resid))
343	# negelbo += n / 2 * tf.reduce_sum(lsigma2s)
344	#
345	# # Regularization
346	# negelbo += pc['lLmb'] * tf.reduce_sum(tf.square(lLmb)) + \
347	# pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0))
348	# negelbo += -tf.reduce_sum(tf.math.log(lnugGPs + 100))
349	#
350	# negelbo /= tf.cast(n, tf.float64)
351	#
352	# return negelbo
353	#
354	# def negproflik(self):
355	# lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()
356	# x = self.x
357	# y = self.y
358	#
359	# pc = self.penalty_const
360	#
361	# n = self.n
362	# q = self.q
363	# D = self.diag_D
364	# phi = self.phi
365	# psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))
366	#
367	# B = tf.matmul(tf.transpose(y / tf.sqrt(tf.exp(lsigma2s))), self.phi)
368	# G = tf.zeros([self.q, n], dtype=tf.float64)
369	#
370	# negproflik = tf.constant(0., dtype=tf.float64)
371	#
372	# for k in range(q):
373	# Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])
374	# Wk, Uk = tf.linalg.eigh(Ck)
375	#
376	# dkInpCkinv_inv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1 / (D[k] + 1 / Wk)), tf.transpose(Uk)))
377	# Gk = tf.matmul(dkInpCkinv_inv, tf.transpose(B)[k])
378	#
379	# CkinvhGk = tf.matmul(tf.matmul(Uk, tf.linalg.diag(1 / Wk)), tf.transpose(Uk)) @ Gk
380	#
381	# G = tf.tensor_scatter_nd_update(G, [[k]], tf.expand_dims(Gk, axis=0))
382	#
383	# negproflik += 0.5 * tf.reduce_sum(tf.math.log(Wk))
384	# negproflik += 0.5 * tf.reduce_sum(tf.square(CkinvhGk))
385	#
386	# resid = (tf.transpose(y) - tf.matmul(tf.transpose(G), psi)) / tf.sqrt(tf.exp(lsigma2s))
387	#
388	# negproflik += 0.5 * tf.reduce_sum(tf.square(resid))
389	# negproflik += n / 2 * tf.reduce_sum(lsigma2s)
390	#
391	# negproflik += pc['lLmb'] * tf.reduce_sum(tf.square(lLmb)) + \
392	# pc['lLmb0'] * (2 / n) * tf.reduce_sum(tf.square(lLmb0))
393	# negproflik += -tf.reduce_sum(tf.math.log(lnugGPs + 100))
394	#
395	# negproflik /= tf.cast(n, tf.float64)
396	# return negproflik
397
398
399	def predict(self, x0, return_fullcov=False):	4✔
400	"""
401	Returns predictive quantities at new input `x0`. Both outputs are of
402	size (number of new input, output dimension).
403	:param x0: New input of size (number of new input, dimension of input).
404	:param return_fullcov: Returns (predictive mean, predictive variance,
405	variance for the true mean, full predictive covariance) if True. Otherwise,
406	only return the first three quantities.
407	"""
408	x0 = self._verify_data_types(x0)	4✔
409	submethod = self.submethod	4✔
410	predict_map = self.submethod_predict_map	4✔
411	try:	4✔
412	predict_call = predict_map[submethod]	4✔
413	except KeyError as e:	×
414	print(e)	×
415	# print('Invalid submethod. Choices are \'full\', \'elbo\', or \'proflik\'.')
NEW 416	raise KeyError('Invalid submethod. Choices are \'full\'.')	×
417	return tf.stop_gradient(predict_call(x0=x0, return_fullcov=return_fullcov))	4✔
418
419
420	def compute_aux_predictive_quantities(self):	4✔
421	"""
422	Compute auxiliary quantities for predictions using full posterior approach.
423	"""
424	x = self.x	4✔
425	lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()	4✔
426
427	D = self.diag_D	4✔
428	# B := Y @ Sigma^{-1/2} @ Phi
429	B = tf.matmul(tf.transpose(self.y) / tf.sqrt(tf.exp(lsigma2s)), self.phi)	4✔
430
431	CinvM = tf.zeros([self.q, self.n], dtype=tf.float64)	4✔
432	Th = tf.zeros([self.q, self.n, self.n], dtype=tf.float64)	4✔
433
434	for k in range(self.q):	4✔
435	Ck = Matern32(x, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k])	4✔
436
437	Wk, Uk = tf.linalg.eigh(Ck)	4✔
438
439	# (I + D_k * C_k)^{-1}
440	IpdkCkinv = tf.matmul(Uk, tf.matmul(tf.linalg.diag(1.0 / (1.0 + D[k] * Wk)), tf.transpose(Uk)))	4✔
441
442	CkinvMk = tf.linalg.matvec(IpdkCkinv, tf.transpose(B)[k])	4✔
443	Thk = tf.matmul(Uk, tf.matmul(tf.linalg.diag(tf.sqrt((D[k] * Wk 2) / (Wk 2 + D[k] * Wk ** 3))),	4✔
444	tf.transpose(Uk)))
445
446	CinvM = tf.tensor_scatter_nd_update(CinvM, [[k]], tf.expand_dims(CkinvMk, axis=0))	4✔
447	Th = tf.tensor_scatter_nd_update(Th, [[k]], tf.expand_dims(Thk, axis=0))	4✔
448
449	self.CinvMs = CinvM	4✔
450	self.Ths = Th	4✔
451
452	def predict_full(self, x0, return_fullcov=False):	4✔
453	"""
454	Returns predictions using full posterior approach.
455	"""
456	if tf.reduce_any(tf.math.is_nan(self.CinvMs)) or tf.reduce_any(tf.math.is_nan(self.Ths)):	4✔
457	self.compute_aux_predictive_quantities()	4✔
458
459	x = self.x	4✔
460	lLmb, lLmb0, lsigma2s, lnugGPs = self.get_param()	4✔
461
462	phi = self.phi	4✔
463
464	CinvM = self.CinvMs	4✔
465	Th = self.Ths	4✔
466
467	x0 = self._verify_data_types(x0)	4✔
468	x0 = (x0 - self.x_min) / (self.x_max - self.x_min) # Standardize x0	4✔
469	n0 = tf.shape(x0)[0]	4✔
470
471	ghat = tf.zeros([self.q, n0], dtype=tf.float64)	4✔
472	gvar = tf.zeros([self.q, n0], dtype=tf.float64)	4✔
473	for k in range(self.q):	4✔
474	c00k = Matern32(x0, x0, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k],	4✔
475	diag_only=True) # Diagonal-only covariance
476	c0k = Matern32(x0, x, llmb=lLmb[k], llmb0=lLmb0[k], lnug=lnugGPs[k],	4✔
477	diag_only=False)
478
479	ghat_k = tf.linalg.matvec(c0k, CinvM[k])	4✔
480	gvar_k = c00k - tf.reduce_sum(tf.square(tf.matmul(c0k, Th[k])), axis=1)	4✔
481
482	ghat = tf.tensor_scatter_nd_update(ghat, [[k]], [ghat_k])	4✔
483	gvar = tf.tensor_scatter_nd_update(gvar, [[k]], [gvar_k])	4✔
484
485	self.ghat = ghat	4✔
486	self.gvar = gvar	4✔
487
488	psi = tf.transpose(phi) * tf.sqrt(tf.exp(lsigma2s))	4✔
489
490	predmean = tf.matmul(psi, ghat, transpose_a=True)	4✔
491	confvar = tf.matmul(tf.transpose(gvar), tf.square(psi))	4✔
492	predvar = confvar + tf.exp(lsigma2s)	4✔
493
494	ypred = self.tx_y(predmean)	4✔
495	yconfvar = tf.transpose(confvar) * tf.square(self.ystd)	4✔
496	ypredvar = tf.transpose(predvar) * tf.square(self.ystd)	4✔
497
498	if return_fullcov:	4✔
NEW 499	CH = tf.sqrt(gvar)[..., tf.newaxis] * psi[tf.newaxis, ...]	×
NEW 500	yfullpredcov = (tf.einsum('nij,njk->nik', CH,	×
501	tf.transpose(CH, perm=[0, 2, 1])) +
502	tf.linalg.diag(tf.exp(lsigma2s)))
NEW 503	yfullpredcov *= tf.square(self.ystd)	×
UNCOV 504	return ypred, ypredvar, yconfvar, yfullpredcov	×
505
506	return ypred, ypredvar, yconfvar	4✔
507
508	@staticmethod	4✔
509	def _verify_data_types(t):	4✔
510	"""
511	Verify if inputs are TensorFlow tensors, if not, cast into tensors.
512	Verify if inputs are at least 2-dimensional, if not, expand dimensions to 2.
513	"""
514	if not isinstance(t, tf.Tensor):	4✔
515	t = tf.convert_to_tensor(t, dtype=tf.float64)	4✔
516	if t.ndim < 2:	4✔
517	t = tf.expand_dims(t, axis=1)	4✔
518	return t	4✔
519
520	# def predict_elbo(self, x0, return_fullcov=False):
521	# pass
522	#
523	# def predict_proflik(self, x0, return_fullcov=False):
524	# pass
525
526
527	@staticmethod	4✔
528	def verify_error_structure(diag_error_structure, y):	4✔
529	"""
530	Verifies if diagonal error structure input, if any, is valid.
531	"""
532	assert sum(diag_error_structure) == y.shape[0], \	4✔
533	'Sum of error_structure should' \
534	' equal the output dimension.'
535
536	def get_param(self):	4✔
537	"""
538	Returns the parameters for LCGP instance.
539	"""
540	# if self.parameter_clamp_flag:
541	# lLmb, lLmb0, lsigma2s, lnugGPs = \
542	# self.parameter_clamp(lLmb=self.lLmb, lLmb0=self.lLmb0,
543	# lsigma2s=self.lsigma2s, lnugs=self.lnugGPs)
544	# else:
545	lLmb, lLmb0, lsigma2s, lnugGPs = \	4✔
546	self.lLmb, self.lLmb0, self.lsigma2s, self.lnugGPs
547
548	built_lsigma2s = tf.zeros(self.p, dtype=tf.float64)	4✔
549	err_struct = self.diag_error_structure	4✔
550	col = 0	4✔
551	for k in range(len(err_struct)):	4✔
552	built_lsigma2s = tf.tensor_scatter_nd_update(	4✔
553	built_lsigma2s,
554	tf.range(col, col + err_struct[k])[:, tf.newaxis],
555	tf.fill([err_struct[k]], lsigma2s[k])
556	)
557	col += err_struct[k]	4✔
558
559	return lLmb, lLmb0, built_lsigma2s, lnugGPs	4✔

mosesyhc / LCGP / 16945523497

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous