def _psicomputations(variance, lengthscale, Z, variational_posterior):
        # here are the "statistics" for psi0, psi1 and psi2
        # Produced intermediate results:
        # _psi1                NxM
        mu = variational_posterior.mean
        S = variational_posterior.variance
        gamma = variational_posterior.binary_prob

        N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
        l2 = np.square(lengthscale)
        log_denom1 = np.log(S/l2+1)
        log_denom2 = np.log(2*S/l2+1)
        log_gamma = np.log(gamma)
        log_gamma1 = np.log(1.-gamma)
        variance = float(variance)
        psi0 = np.empty(N)
        psi0[:] = variance
        psi1 = np.empty((N,M))
        psi2n = np.empty((N,M,M))

        from ....util.misc import param_to_array
        S = param_to_array(S)
        mu = param_to_array(mu)
        gamma = param_to_array(gamma)
        Z = param_to_array(Z)

        support_code = """
        #include <math.h>
        """
        code = """
        for(int n=0; n<N; n++) {
            for(int m1=0;m1<M;m1++) {
                double log_psi1=0;
                for(int m2=0;m2<=m1;m2++) {
                    double log_psi2_n=0;
                    for(int q=0;q<Q;q++) {
                        double Snq = S(n,q);
                        double lq = l2(q);
                        double Zm1q = Z(m1,q);
                        double Zm2q = Z(m2,q);

                        if(m2==0) {
                            // Compute Psi_1
                            double muZ = mu(n,q)-Z(m1,q);

                            double psi1_exp1 = log_gamma(n,q) - (muZ*muZ/(Snq+lq) +log_denom1(n,q))/2.;
                            double psi1_exp2 = log_gamma1(n,q) -Zm1q*Zm1q/(2.*lq);
                            log_psi1 += (psi1_exp1>psi1_exp2)?psi1_exp1+log1p(exp(psi1_exp2-psi1_exp1)):psi1_exp2+log1p(exp(psi1_exp1-psi1_exp2));
                        }
                        // Compute Psi_2
                        double muZhat = mu(n,q) - (Zm1q+Zm2q)/2.;
                        double Z2 = Zm1q*Zm1q+ Zm2q*Zm2q;
                        double dZ = Zm1q - Zm2q;

                        double psi2_exp1 = dZ*dZ/(-4.*lq)-muZhat*muZhat/(2.*Snq+lq) - log_denom2(n,q)/2. + log_gamma(n,q);
                        double psi2_exp2 = log_gamma1(n,q) - Z2/(2.*lq);
                        log_psi2_n += (psi2_exp1>psi2_exp2)?psi2_exp1+log1p(exp(psi2_exp2-psi2_exp1)):psi2_exp2+log1p(exp(psi2_exp1-psi2_exp2));
                    }
                    double exp_psi2_n = exp(log_psi2_n);
                    psi2n(n,m1,m2) = variance*variance*exp_psi2_n;
                    if(m1!=m2) { psi2n(n,m2,m1) = variance*variance*exp_psi2_n;}
                }
                psi1(n,m1) = variance*exp(log_psi1);
            }
        }
        """
        weave.inline(code, support_code=support_code, arg_names=['psi1','psi2n','N','M','Q','variance','l2','Z','mu','S','gamma','log_denom1','log_denom2','log_gamma','log_gamma1'], type_converters=weave.converters.blitz)

        psi2 = psi2n.sum(axis=0)
        return psi0,psi1,psi2,psi2n
Z - MxQ
dataset/ETHPy150Open SheffieldML/GPy/GPy/kern/src/psi_comp/ssrbf_psi_comp.py/_psicomputations