Compute the expectation:
    expectation[n] = <Ka_{Z1, x_n} Kb_{x_n, Z2}>_p(x_n)
        - Ka_{.,.}, Kb_{.,.} :: RBF kernels
    Ka and Kb as well as Z1 and Z2 can differ from each other.

    :return: N x dim(Z1) x dim(Z2)
12095
def _expectation(p, kern1, feat1, kern2, feat2, nghp=None):
    

    if kern1.on_separate_dims(kern2) and isinstance(p, DiagonalGaussian):
        # no joint expectations required
        eKxz1 = expectation(p, (kern1, feat1))
        eKxz2 = expectation(p, (kern2, feat2))
        return eKxz1[:, :, None] * eKxz2[:, None, :]

    Ka, Kb = kern1, kern2

    with params_as_tensors_for(Ka, feat1, Kb, feat2):

        # use only active dimensions
        Xcov = Ka._slice_cov(tf.matrix_diag(p.cov) if isinstance(p, DiagonalGaussian) else p.cov)

        Z1, Xmu = Ka._slice(feat1.Z, p.mu)

        N = tf.shape(Xmu)[0]
        D = tf.shape(Xmu)[1]

        def get_squared_length_scales(kern):
            squared_lengthscales = kern.lengthscales ** 2. if kern.ARD \
                        else tf.zeros((D,), dtype=settings.float_type) + kern.lengthscales ** 2.
            return squared_lengthscales

        if Ka == Kb:
            La = get_squared_length_scales(Ka)
            Lb = La
            half_mean_L = La * 0.5  # average length scale
        else:
            La, Lb = map(get_squared_length_scales, (Ka, Kb))
            half_mean_L = La * Lb / (La + Lb)  # average length scale

        sqrt_det_L = tf.reduce_prod(half_mean_L) ** 0.5
        C = tf.cholesky(tf.matrix_diag(half_mean_L) + Xcov)  # [N, D, D]
        dets = sqrt_det_L / tf.exp(tf.reduce_sum(tf.log(tf.matrix_diag_part(C)), axis=1))  # N

        # for mahalanobis computation we need Zᵀ (CCᵀ)⁻¹ Z  as well as C⁻¹ Z
        # with Z = Z₁, Z₂  for two rbf kernels
        def get_cholesky_solve_terms(Z, C=C):
            C_inv_z = tf.matrix_triangular_solve(
                C, tf.tile(tf.expand_dims(tf.transpose(Z), 0),
                           [N, 1, 1]), lower=True)  # [N, D, M]
            z_CC_inv_z = tf.reduce_sum(tf.square(C_inv_z), 1)  # [N, M]

            return C_inv_z, z_CC_inv_z

        C_inv_mu = tf.matrix_triangular_solve(C, tf.expand_dims(Xmu, 2), lower=True)  # [N, D, 1]
        mu_CC_inv_mu = tf.expand_dims(tf.reduce_sum(tf.square(C_inv_mu), 1), 2)  # [N, 1, 1]

        C_inv_z1, z1_CC_inv_z1 = get_cholesky_solve_terms(Z1 / La * half_mean_L)
        z1_CC_inv_mu = 2 * tf.matmul(C_inv_z1, C_inv_mu, transpose_a=True)[:, :, 0]  # [N, M1]

        if feat1 == feat2 and Ka == Kb:
            # in this case Z2==Z1 so we can reuse the Z1 terms
            C_inv_z2, z2_CC_inv_z2  = C_inv_z1, z1_CC_inv_z1
            z2_CC_inv_mu = z1_CC_inv_mu  # [N, M]
            Z2 = Z1
        else:
            # compute terms related to Z2
            Z2, _ = Kb._slice(feat2.Z, p.mu)
            C_inv_z2, z2_CC_inv_z2 = get_cholesky_solve_terms(Z2 / Lb * half_mean_L)
            z2_CC_inv_mu = 2 * tf.matmul(C_inv_z2, C_inv_mu, transpose_a=True)[:, :, 0]  # [N, M2]

        z1_CC_inv_z2 = tf.matmul(C_inv_z1, C_inv_z2, transpose_a=True)  # [N, M1, M2]

        # expand dims for broadcasting
        # along M1
        z2_CC_inv_mu = tf.expand_dims(z2_CC_inv_mu, 1)  # [N, 1, M2]
        z2_CC_inv_z2 = tf.expand_dims(z2_CC_inv_z2, 1)

        # along M2
        z1_CC_inv_mu = tf.expand_dims(z1_CC_inv_mu, 2)  # [N, M1, 1]
        z1_CC_inv_z1 = tf.expand_dims(z1_CC_inv_z1, 2)

        # expanded version of ((Z1 + Z2)-mu) (CCT)-1 ((Z1 + Z2)-mu)
        mahalanobis = mu_CC_inv_mu + z2_CC_inv_z2 + \
                      z1_CC_inv_z1 + 2 * z1_CC_inv_z2 - \
                      z1_CC_inv_mu - z2_CC_inv_mu  # [N, M1, M2]

        exp_mahalanobis = tf.exp(-0.5 * mahalanobis)  # [N, M1, M2]

        if Z1 == Z2:
            # CAVEAT : Compute sqrt(self.K(Z)) explicitly
            # to prevent automatic gradient from
            # being NaN sometimes, see pull request #615
            sqrt_exp_dist = tf.exp(-0.25 * Ka.scaled_square_dist(Z1, None))
        else:
            # Compute exp( -.5 (Z-Z')^top (L_1+L_2)^{-1} (Z-Z') )
            lengthscales_rms = tf.sqrt(La + Lb)
            Z1 = Z1 / lengthscales_rms
            Z1sqr = tf.reduce_sum(tf.square(Z1), axis=1)
            Z2 = Z2 / lengthscales_rms
            Z2sqr = tf.reduce_sum(tf.square(Z2), axis=1)
            dist = -2 * tf.matmul(Z1, Z2, transpose_b=True) \
                   + tf.reshape(Z1sqr, (-1, 1)) + tf.reshape(Z2sqr, (1, -1))
            sqrt_exp_dist = tf.exp(-0.5 * dist)  # M1 x M2

        return Ka.variance * Kb.variance * sqrt_exp_dist * \
               tf.reshape(dets, [N, 1, 1]) * exp_mahalanobis
394987
def _M(self, k, t, tips=None, gaps=None):
        """"""
        if self._distributionmodel:
            return self.model.M(k, t, tips, gaps)
        else:
            return self.model.M(t, tips, gaps)
557383
def _build_likelihood(self):
        
        pX = DiagonalGaussian(self.X_mean, self.X_var)

        num_inducing = len(self.feature)
        psi0 = tf.reduce_sum(expectation(pX, self.kern))
        psi1 = expectation(pX, (self.kern, self.feature))
        psi2 = tf.reduce_sum(expectation(pX, (self.kern, self.feature), (self.kern, self.feature)), axis=0)
        Kuu = features.Kuu(self.feature, self.kern, jitter=settings.jitter)
        L = tf.cholesky(Kuu)
        sigma2 = self.likelihood.variance
        sigma = tf.sqrt(sigma2)

        # Compute intermediate matrices
        A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
        tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
        AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
        B = AAT + tf.eye(num_inducing, dtype=settings.float_type)
        LB = tf.cholesky(B)
        log_det_B = 2. * tf.reduce_sum(tf.log(tf.matrix_diag_part(LB)))
        c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma

        # KL[q(x) || p(x)]
        dX_var = self.X_var if len(self.X_var.get_shape()) == 2 else tf.matrix_diag_part(self.X_var)
        NQ = tf.cast(tf.size(self.X_mean), settings.float_type)
        D = tf.cast(tf.shape(self.Y)[1], settings.float_type)
        KL = -0.5 * tf.reduce_sum(tf.log(dX_var)) \
             + 0.5 * tf.reduce_sum(tf.log(self.X_prior_var)) \
             - 0.5 * NQ \
             + 0.5 * tf.reduce_sum((tf.square(self.X_mean - self.X_prior_mean) + dX_var) / self.X_prior_var)

        # compute log marginal bound
        ND = tf.cast(tf.size(self.Y), settings.float_type)
        bound = -0.5 * ND * tf.log(2 * np.pi * sigma2)
        bound += -0.5 * D * log_det_B
        bound += -0.5 * tf.reduce_sum(tf.square(self.Y)) / sigma2
        bound += 0.5 * tf.reduce_sum(tf.square(c))
        bound += -0.5 * D * (tf.reduce_sum(psi0) / sigma2 -
                             tf.reduce_sum(tf.matrix_diag_part(AAT)))
        bound -= KL
        return bound
935018
def _nbinom_ztrunc_p(mu, k_agg):
        

        p_eq = lambda p, mu, k_agg: (k_agg * p) / (1 - (1 + p)**-k_agg) - mu

        # The upper bound needs to be large. p will increase with increasing mu
        # and decreasing k_agg
        p = optim.brentq(p_eq, 1e-10, 1e10, args=(mu, k_agg))
        return p
990307
def kernelize(self, kernel):
    '''
    '''
    if kernel == 'none':
      return self
    if kernel == 'binary':
      if self.is_weighted():
        return self._update_edges(1, copy=True)
      return self
    if kernel == 'rbf':
      w = self.edge_weights()
      r = np.exp(-w / w.std())
      return self._update_edges(r, copy=True)
    raise ValueError('Invalid kernel type: %r' % kernel)
1314859
def _expectation(p, kern, feat, none2, none3, nghp=None):
    
    """
    return functools.reduce(tf.add, [
        expectation(p, (k, feat), nghp=nghp) for k in kern.kernels])
378778
def log_K(Z, X, t):
    
    m, L = update_params(X, t)
    z = Z['x']
    return np.log(det(L))-la.norm(np.dot(L.T,z-m))**2/2.
1333466
def KGE(s, o):
    
    # s,o = filter_nan(s, o)
    cc = correlation(s, o)
    alpha = np.std(s)/np.std(o)
    beta = np.sum(s)/np.sum(o)
    kge = 1 - np.sqrt((cc-1)**2 + (alpha-1)**2 + (beta-1)**2)
    return kge, cc, alpha, beta
228759
def est_kl_divergence(self, other, kernel=None, delta=1e-2):
        
        return self._kl_divergence(
            other.particle_locations,
            other.particle_weights,
            kernel, delta
        )
670785
def _build_predict(self, Xnew, full_cov=False):
        
        num_inducing = len(self.feature)
        err = self.Y - self.mean_function(self.X)
        Kuf = features.Kuf(self.feature, self.kern, self.X)
        Kuu = features.Kuu(self.feature, self.kern, jitter=settings.numerics.jitter_level)
        Kus = features.Kuf(self.feature, self.kern, Xnew)
        sigma = tf.sqrt(self.likelihood.variance)
        L = tf.cholesky(Kuu)
        A = tf.matrix_triangular_solve(L, Kuf, lower=True) / sigma
        B = tf.matmul(A, A, transpose_b=True) + tf.eye(num_inducing, dtype=settings.float_type)
        LB = tf.cholesky(B)
        Aerr = tf.matmul(A, err)
        c = tf.matrix_triangular_solve(LB, Aerr, lower=True) / sigma
        tmp1 = tf.matrix_triangular_solve(L, Kus, lower=True)
        tmp2 = tf.matrix_triangular_solve(LB, tmp1, lower=True)
        mean = tf.matmul(tmp2, c, transpose_a=True)
        if full_cov:
            var = self.kern.K(Xnew) + tf.matmul(tmp2, tmp2, transpose_a=True) \
                  - tf.matmul(tmp1, tmp1, transpose_a=True)
            var = tf.tile(var[None, ...], [self.num_latent, 1, 1])  # P x N x N
        else:
            var = self.kern.Kdiag(Xnew) + tf.reduce_sum(tf.square(tmp2), 0) \
                  - tf.reduce_sum(tf.square(tmp1), 0)
            var = tf.tile(var[:, None], [1, self.num_latent])
        return mean + self.mean_function(Xnew), var
1328166
def _expectation(p, linear_mean, none, kern, feat, nghp=None):
    
    with params_as_tensors_for(linear_mean):
        N = p.mu.shape[0].value
        D = p.mu.shape[1].value
        exKxz = expectation(p, mean_functions.Identity(D), (kern, feat), nghp=nghp)
        eKxz = expectation(p, (kern, feat), nghp=nghp)
        eAxKxz = tf.matmul(tf.tile(linear_mean.A[None, :, :], (N, 1, 1)),
                           exKxz, transpose_a=True)
        ebKxz = linear_mean.b[None, :, None] * eKxz[:, None, :]
        return eAxKxz + ebKxz
518376
def inference_likelihood(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape
        input_dim = Z.shape[0]
        if self.mpi_comm != None:
            from mpi4py import MPI
            num_data_all = np.array(num_data,dtype=np.int32)
            self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])
            num_data = num_data_all

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        if het_noise:
            self.batchsize = 1

        psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kmm = kern.K(Z).copy()
        diag.add(Kmm, self.const_jitter)
        if not np.isfinite(Kmm).all():
            print(Kmm)
        Lm = jitchol(Kmm)
        LmInv = dtrtri(Lm)

        LmInvPsi2LmInvT = LmInv.dot(psi2_full.dot(LmInv.T))
        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LLInv = dtrtri(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        LmLLInv = LLInv.dot(LmInv)
        
        b  = psi1Y_full.dot(LmLLInv.T)
        bbt = np.square(b).sum()
        v = b.dot(LmLLInv).T
        LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)
        
        tmp = -LLInv.T.dot(LLinvPsi1TYYTPsi1LLinvT+output_dim*np.eye(input_dim)).dot(LLInv)
        dL_dpsi2R = LmInv.T.dot(tmp+output_dim*np.eye(input_dim)).dot(LmInv)/2.
        
        # Cache intermediate results
        self.midRes['dL_dpsi2R'] = dL_dpsi2R
        self.midRes['v'] = v

        #======================================================================
        # Compute log-likelihood
        #======================================================================
        if het_noise:
            logL_R = -np.log(beta).sum()
        else:
            logL_R = -num_data*np.log(beta)
        logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        dL_dKmm =  dL_dpsi2R - output_dim*LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        if not self.Y_speedup or het_noise:
            wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
            post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
        else:
            post = None

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        if not het_noise:
            dL_dthetaL = (YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta)/2. - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()
            self.midRes['dL_dthetaL'] = dL_dthetaL

        return logL, dL_dKmm, post
101609
def _build_likelihood(self):
        

        # Get prior KL.
        KL = self.build_prior_KL()

        # Get conditionals
        fmean, fvar = self._build_predict(self.X, full_cov=False, full_output_cov=False)

        # Get variational expectations.
        var_exp = self.likelihood.variational_expectations(fmean, fvar, self.Y)

        # re-scale for minibatch size
        scale = tf.cast(self.num_data, settings.float_type) / tf.cast(tf.shape(self.X)[0], settings.float_type)

        return tf.reduce_sum(var_exp) * scale - KL
309657
def _expectation(p, kern, feat, none1, none2, nghp=None):
    
    with params_as_tensors_for(kern, feat):
        # use only active dimensions
        Xcov = kern._slice_cov(p.cov)
        Z, Xmu = kern._slice(feat.Z, p.mu)
        D = tf.shape(Xmu)[1]
        if kern.ARD:
            lengthscales = kern.lengthscales
        else:
            lengthscales = tf.zeros((D,), dtype=settings.float_type) + kern.lengthscales

        chol_L_plus_Xcov = tf.cholesky(tf.matrix_diag(lengthscales ** 2) + Xcov)  # NxDxD

        all_diffs = tf.transpose(Z) - tf.expand_dims(Xmu, 2)  # NxDxM
        exponent_mahalanobis = tf.matrix_triangular_solve(chol_L_plus_Xcov, all_diffs, lower=True)  # NxDxM
        exponent_mahalanobis = tf.reduce_sum(tf.square(exponent_mahalanobis), 1)  # NxM
        exponent_mahalanobis = tf.exp(-0.5 * exponent_mahalanobis)  # NxM

        sqrt_det_L = tf.reduce_prod(lengthscales)
        sqrt_det_L_plus_Xcov = tf.exp(tf.reduce_sum(tf.log(tf.matrix_diag_part(chol_L_plus_Xcov)), axis=1))
        determinants = sqrt_det_L / sqrt_det_L_plus_Xcov  # N

        return kern.variance * (determinants[:, None] * exponent_mahalanobis)
203797
def _kp(a, b):
    
    if a.shape != b.shape or a.shape[-1] != 1:
        raise(ValueError)
    N = a.shape[0]
    # take the outer product over the last two axes, then reshape:
    return np.einsum('ijk,ilk->ijkl', a, b).reshape(N, -1, 1)
370801
def _expectation(self, X):
        """"""

        # Calculate probabilities of X belonging to the different clusters
        weighted_likelihoods = self._get_likelihoods(X) * self.priors
        sum_likelihoods = np.expand_dims(np.sum(weighted_likelihoods, axis=1), axis=1)
        # Determine responsibility as P(X|y)*P(y)/P(X)
        self.responsibility = weighted_likelihoods / sum_likelihoods
        # Assign samples to cluster that has largest probability
        self.sample_assignments = self.responsibility.argmax(axis=1)
        # Save value for convergence check
        self.responsibilities.append(np.max(self.responsibility, axis=1))
442685
def predict_covariance(self, X, with_noise=True):
        
        _, v = self._predict(X, True, with_noise)
        return v