r"""
        This method computes the variational lower bound on the likelihood,
        which is:

            E_{q(F)} [ \log p(Y|F) ] - KL[ q(F) || p(F)]

        with

            q(\mathbf f) = N(\mathbf f \,|\, \boldsymbol \mu, \boldsymbol \Sigma)
774414
def _build_likelihood(self):
        

        """

        # Get prior KL.
        KL = gauss_kl(self.q_mu, self.q_sqrt)

        # Get conditionals
        K = self.kern.K(self.X) + tf.eye(self.num_data, dtype=settings.float_type) * \
            settings.numerics.jitter_level
        L = tf.cholesky(K)

        fmean = tf.matmul(L, self.q_mu) + self.mean_function(self.X)  # NN,ND->ND

        q_sqrt_dnn = tf.matrix_band_part(self.q_sqrt, -1, 0)  # D x N x N

        L_tiled = tf.tile(tf.expand_dims(L, 0), tf.stack([self.num_latent, 1, 1]))

        LTA = tf.matmul(L_tiled, q_sqrt_dnn)  # D x N x N
        fvar = tf.reduce_sum(tf.square(LTA), 2)

        fvar = tf.transpose(fvar)

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

        return tf.reduce_sum(var_exp) - KL
450270
def grad_local_log_likelihood(self, x):
        
        # Observation likelihoods
        lmbda = np.exp(np.dot(x, self.C.T) + np.dot(self.inputs, self.D.T))
        return (self.data - lmbda).dot(self.C)
610643
def log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # prior log p(x) -- quadratic terms
        J_diag, J_upper_diag = self.sparse_J_prior
        lp = -0.5 * np.sum(x * symm_block_tridiag_matmul(J_diag, J_upper_diag, x))

        # prior log p(x) -- linear terms
        _, h_init, log_Z_init = self.info_init_params
        _, _, _, h1, h2, log_Z_dyn = self.info_dynamics_params
        lp += x[0].dot(h_init)
        lp += np.sum(x[:-1] * h1)
        lp += np.sum(x[1:] * h2)

        # prior log p(x) -- normalization constants
        lp += log_Z_init
        lp += np.sum(log_Z_dyn)

        # likelihood log p(y | x)
        lp += self.log_conditional_likelihood(x)

        return lp
1187803
def hessian_local_log_likelihood(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        hfun = hessian(self.local_log_likelihood)
        H_diag = np.zeros((T, D, D))
        for t in range(T):
            H_diag[t] = hfun(x[t], self.data[t], self.inputs[t])
        return H_diag
791801
def log_likelihood_pairwise(data, params):
    """"""
    loglik = 0
    for winner, loser in data:
        loglik -= np.logaddexp(0, -(params[winner] - params[loser]))
    return loglik
1151312
def sparse_hessian_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        J_diag, J_upper_diag = self.sparse_J_prior
        H_diag, H_upper_diag = -J_diag, -J_upper_diag

        # Collect the likelihood terms
        H_diag += self.hessian_local_log_likelihood(x)

        # Subtract a little bit to ensure negative definiteness
        H_diag -= 1e-8 * np.eye(D)

        return H_diag, H_upper_diag
917828
def inference_minibatch(self, kern, X, Z, likelihood, Y, Zp):
        

        num_data, output_dim = Y.shape

        #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
        # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
        #self.YYTfactor = beta*self.get_YYTfactor(Y)
        if self.Y_speedup and not het_noise:
            YYT_factor = self.get_YYTfactor(Y)
        else:
            YYT_factor = Y

        n_start = self.batch_pos
        batchsize = num_data if self.batchsize is None else self.batchsize
        n_end = min(batchsize+n_start, num_data)
        if n_end == num_data:
            isEnd = True
            self.batch_pos = 0
        else:
            isEnd = False
            self.batch_pos = n_end

        if batchsize == num_data:
            Y_slice = YYT_factor
            X_slice = X
        else:
            Y_slice = YYT_factor[n_start:n_end]
            X_slice = X[n_start:n_end]

        psi0 = kern._Kdiag(X_slice) #Kffdiag
        psi1 = kern.K(X_slice, Z) #Kfu
        betapsi1 = np.einsum('n,nm->nm', beta, psi1)

        X_slice = X_slice.values
        Z = Z.values

        Zp = Zp.gamma.values
        indX = np.int_(X_slice[:, -1])
        indZ = np.int_(Z[:, -1]) - Zp.shape[0]

        betaY = beta*Y_slice

        #======================================================================
        # Load Intermediate Results
        #======================================================================

        dL_dpsi2R = self.midRes['dL_dpsi2R']
        v = self.midRes['v']

        #======================================================================
        # Compute dL_dpsi
        #======================================================================

        dL_dpsi0 = -.5*output_dim*(beta * Zp[indX, :]) #XxQ #TODO: Check this gradient

        dL_dpsi1 = np.dot(betaY, v.T)
        dL_dEZp  = psi1*dL_dpsi1
        dL_dpsi1 = Zp[np.ix_(indX, indZ)]*dL_dpsi1
        dL_dgamma = np.zeros(Zp.shape)
        for d in np.unique(indX):
            indd = indX == d
            betapsi1d = betapsi1[indd, :]
            psi1d = psi1[indd, :]
            Zpd = Zp[d, :]
            Zp2 = Zpd[:, None]*Zpd[None, :] - np.diag(np.power(Zpd, 2)) + np.diag(Zpd)
            dL_dpsi1[indd, :] += np.dot(betapsi1d, Zp2[np.ix_(indZ, indZ)] * dL_dpsi2R)*2.

            dL_EZp2 = dL_dpsi2R * (np.dot(psi1d.T, psi1d) * beta)*2.  # Zpd*Kufd*Kfud*beta
            #Gradient of Likelihood wrt gamma is calculated here
            EZ = Zp[d, indZ]
            for q in range(Zp.shape[1]):
                EZt = EZ.copy()
                indq = indZ == q
                EZt[indq] = .5
                dL_dgamma[d, q] = np.sum(dL_dEZp[np.ix_(indd, indq)]) + np.sum(dL_EZp2[:, indq]*EZt[:, None]) -\
                    .5*beta*(np.sum(psi0[indd, q]))

        #======================================================================
        # Compute dL_dthetaL
        #======================================================================
        if isEnd:
            dL_dthetaL = self.midRes['dL_dthetaL']
        else:
            dL_dthetaL = 0.

        grad_dict = {'dL_dKdiag': dL_dpsi0,
                     'dL_dKnm': dL_dpsi1,
                     'dL_dthetaL': dL_dthetaL,
                     'dL_dgamma': dL_dgamma}

        return isEnd, (n_start, n_end), grad_dict
956087
def gradient_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        _, h_init, _ = self.info_init_params
        _, _, _, h1, h2, _ = self.info_dynamics_params
        H_diag, H_upper_diag = self.sparse_J_prior

        # Compute the gradient from the prior
        g = -1 * symm_block_tridiag_matmul(H_diag, H_upper_diag, x)
        g[0] += h_init
        g[:-1] += h1
        g[1:] += h2

        # Compute gradient from the likelihood terms
        g += self.grad_local_log_likelihood(x)

        return g
383666
def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None):
        
        assert mean_function is None, "inference with a mean function not implemented"

        # Compute K
        K = kern.K(X)

        #Find mode
        if self.bad_fhat or self.first_run:
            Ki_f_init = np.zeros_like(Y)
            self.first_run = False
        else:
            Ki_f_init = self._previous_Ki_fhat

        Ki_f_init = np.zeros_like(Y)# FIXME: take this out

        f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)

        #Compute hessian and other variables at mode
        log_marginal, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata)

        self._previous_Ki_fhat = Ki_fhat.copy()
        return Posterior(woodbury_vector=Ki_fhat, woodbury_inv=woodbury_inv, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL}
138527
def log_prior(self):
        """"""
        for p, b in zip(self.parameter_vector, self.parameter_bounds):
            if b[0] is not None and p < b[0]:
                return -np.inf
            if b[1] is not None and p > b[1]:
                return -np.inf
        return 0.0
397700
def likelihood(self, x, cl):
        

        # sha = x.shape
        # xr = x.reshape(-1, sha[-1])
        # outsha = sha[:-1]
        # from PyQt4.QtCore import pyqtRemoveInputHook
        # pyqtRemoveInputHook()
        logger.debug("likel " + str(x.shape))
        if self.modelparams["type"] == "gmmsame":

            px = self.mdl[cl].score_samples(x)

        # todo ošetřit více dimenzionální fv
        # px = px.reshape(outsha)
        elif self.modelparams["type"] == "kernel":
            px = self.mdl[cl].score_samples(x)
        elif self.modelparams["type"] == "gaussian_kde":
            # print x
            # np.log because it is likelihood
            # @TODO Zde je patrně problém s reshape
            # old
            # px = np.log(self.mdl[cl](x.reshape(-1)))
            # new
            px = np.log(self.mdl[cl](x))
            # px = px.reshape(outsha)
            # from PyQt4.QtCore import pyqtRemoveInputHook
            # pyqtRemoveInputHook()
        elif self.modelparams["type"] == "dpgmm":
            # todo here is a hack
            # dpgmm z nějakého důvodu nefunguje pro naše data
            # vždy natrénuje jednu složku v blízkosti nuly
            # patrně to bude mít něco společného s parametrem alpha
            # přenásobí-li se to malým číslem, zázračně to chodí
            logger.warning(".score() replaced with .score_samples() . Check it.")
            # px = self.mdl[cl].score(x * 0.01)
            px = self.mdl[cl].score_samples(x * 0.01)
        elif self.modelparams["type"] == "stored":
            px = self.mdl[cl].score(x)
        return px
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
1202326
def log_gaussian(x, y, x_sigma, y_sigma, mu):
    
    if x_sigma==0.0 or y_sigma==0.0:
        return x * 0.0

    with float_error_ignore():
        x_w = np.divide(np.log(x)-mu, x_sigma*x_sigma)
        y_h = np.divide(np.log(y)-mu, y_sigma*y_sigma)

        return np.exp(-0.5*x_w*x_w + -0.5*y_h*y_h)
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.
382572
function ecp_proj(R, N, p, l) {
	var x = N[0], y = N[1], z = N[2];
	var xr = R[0], yr = R[1], zr = R[2];
	var t0 = _ecp_t0, t1 = _ecp_t1, t2 = _ecp_t2;
	fp_mul(t2, t0, fp_sqr(t1, fp_inv(t0, z, p, l), p, l), p, l);
	fp_mul(xr, x, t1, p, l);
	fp_mul(yr, y, t2, p, l);
	mpn_zero(zr, l), zr[0] = 1;
	return R;
}
645707
def log_likelihood(self):
        
        ll = GP.log_likelihood(self)
        jacobian = self.warping_function.fgrad_y(self.Y_untransformed)
        return ll + np.log(jacobian).sum()
1316092
def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
    '''
    '''
    if not numpy.all(numpy.array(deltas) > epsilon):
        raise ArgumentError('all deltas have to be greater than epsilon')

    bound = []
    for delta in deltas:
        # get boundary paths
        paths = pseudo.contour_paths(delta)

        # get vertices on boundary
        vertices = paths.vertices()

        # evaluate polynomial
        supremum = numpy.max(numpy.abs(p(vertices)))

        # compute bound
        bound.append(epsilon/(delta-epsilon)
                     * paths.length()/(2*numpy.pi*delta)
                     * supremum)
    return bound