def inference_likelihood(self, kern, X, Z, likelihood, Y, Zp):
num_data, output_dim = Y.shape<EOL>input_dim = Z.shape[<NUM_LIT:0>]<EOL>if self.mpi_comm is not None:<EOL><INDENT>from mpi4py import MPI<EOL>num_data_all = np.array(num_data,dtype=np.int32)<EOL>self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])<EOL>num_data = num_data_all<EOL><DEDENT>beta = <NUM_LIT:1.>/np.fmax(likelihood.variance, <NUM_LIT>)<EOL>het_noise = beta.size > <NUM_LIT:1><EOL>if het_noise:<EOL><INDENT>self.batchsize = <NUM_LIT:1><EOL><DEDENT>psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, Zp)<EOL>Kmm = kern.K(Z).copy()<EOL>diag.add(Kmm, self.const_jitter)<EOL>if not np.isfinite(Kmm).all():<EOL><INDENT>print(Kmm)<EOL><DEDENT>Lm = jitchol(Kmm)<EOL>LmInv = dtrtri(Lm)<EOL>LmInvPsi2LmInvT = np.dot(LmInv, np.dot(psi2_full, LmInv.T))<EOL>Lambda = np.eye(Kmm.shape[<NUM_LIT:0>])+LmInvPsi2LmInvT<EOL>LL = jitchol(Lambda)<EOL>LLInv = dtrtri(LL)<EOL>logdet_L = <NUM_LIT>*np.sum(np.log(np.diag(LL)))<EOL>LmLLInv = np.dot(LLInv, LmInv)<EOL>b = np.dot(psi1Y_full, LmLLInv.T)<EOL>bbt = np.sum(np.square(b))<EOL>v = np.dot(b, LmLLInv).T<EOL>LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)<EOL>tmp = -np.dot(np.dot(LLInv.T, LLinvPsi1TYYTPsi1LLinvT + output_dim*np.eye(input_dim)), LLInv)<EOL>dL_dpsi2R = <NUM_LIT>*np.dot(np.dot(LmInv.T, tmp + output_dim*np.eye(input_dim)), LmInv)<EOL>self.midRes['<STR_LIT>'] = dL_dpsi2R<EOL>self.midRes['<STR_LIT:v>'] = v<EOL>if het_noise:<EOL><INDENT>logL_R = -np.sum(np.log(beta))<EOL><DEDENT>else:<EOL><INDENT>logL_R = -num_data*np.log(beta)<EOL><DEDENT>logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)*<NUM_LIT> - output_dim*logdet_L*<NUM_LIT><EOL>dL_dKmm = dL_dpsi2R - <NUM_LIT>*output_dim*np.dot(np.dot(LmInv.T, LmInvPsi2LmInvT), LmInv)<EOL>if not self.Y_speedup or het_noise:<EOL><INDENT>wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='<STR_LIT:left>'), transpose='<STR_LIT:left>')<EOL>post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)<EOL><DEDENT>else:<EOL><INDENT>post = None<EOL><DEDENT>if not het_noise:<EOL><INDENT>dL_dthetaL = <NUM_LIT>*(YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta) - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()<EOL>self.midRes['<STR_LIT>'] = dL_dthetaL<EOL><DEDENT>return logL, dL_dKmm, post<EOL>
The first phase of inference:
Compute: log-likelihood, dL_dKmm

Cached intermediate results: Kmm, KmmInv,