def inference(self, kern, X, Z, likelihood, Y, indexD, output_dim, Y_metadata=None, Lm=None, dL_dKmm=None, Kuu_sigma=None):
input_dim = Z.shape[<NUM_LIT:0>]<EOL>uncertain_inputs = isinstance(X, VariationalPosterior)<EOL>beta = <NUM_LIT:1.>/likelihood.variance<EOL>if len(beta)==<NUM_LIT:1>:<EOL><INDENT>beta = np.zeros(output_dim)+beta<EOL><DEDENT>beta_exp = np.zeros(indexD.shape[<NUM_LIT:0>])<EOL>for d in range(output_dim):<EOL><INDENT>beta_exp[indexD==d] = beta[d]<EOL><DEDENT>psi0, psi1, psi2 = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)<EOL>psi2_sum = (beta_exp[:,None,None]*psi2).sum(<NUM_LIT:0>)/output_dim<EOL>Kmm = kern.K(Z).copy()<EOL>if Kuu_sigma is not None:<EOL><INDENT>diag.add(Kmm, Kuu_sigma)<EOL><DEDENT>else:<EOL><INDENT>diag.add(Kmm, self.const_jitter)<EOL><DEDENT>Lm = jitchol(Kmm)<EOL>logL = <NUM_LIT:0.><EOL>dL_dthetaL = np.zeros(output_dim)<EOL>dL_dKmm = np.zeros_like(Kmm)<EOL>dL_dpsi0 = np.zeros_like(psi0)<EOL>dL_dpsi1 = np.zeros_like(psi1)<EOL>dL_dpsi2 = np.zeros_like(psi2)<EOL>wv = np.empty((Kmm.shape[<NUM_LIT:0>],output_dim))<EOL>for d in range(output_dim):<EOL><INDENT>idx_d = indexD==d<EOL>Y_d = Y[idx_d]<EOL>N_d = Y_d.shape[<NUM_LIT:0>]<EOL>beta_d = beta[d]<EOL>psi2_d = psi2[idx_d].sum(<NUM_LIT:0>)*beta_d<EOL>psi1Y = Y_d.T.dot(psi1[idx_d])*beta_d<EOL>psi0_d = psi0[idx_d].sum()*beta_d<EOL>YRY_d = np.square(Y_d).sum()*beta_d<EOL>LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_d, '<STR_LIT:right>')<EOL>Lambda = np.eye(Kmm.shape[<NUM_LIT:0>])+LmInvPsi2LmInvT<EOL>LL = jitchol(Lambda)<EOL>LmLL = Lm.dot(LL)<EOL>b  = dtrtrs(LmLL, psi1Y.T)[<NUM_LIT:0>].T<EOL>bbt = np.square(b).sum()<EOL>v = dtrtrs(LmLL, b.T, trans=<NUM_LIT:1>)[<NUM_LIT:0>].T<EOL>LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)<EOL>tmp = -backsub_both_sides(LL, LLinvPsi1TYYTPsi1LLinvT)<EOL>dL_dpsi2R = backsub_both_sides(Lm, tmp+np.eye(input_dim))/<NUM_LIT:2><EOL>logL_R = -N_d*np.log(beta_d)<EOL>logL += -((N_d*log_2_pi+logL_R+psi0_d-np.trace(LmInvPsi2LmInvT))+YRY_d- bbt)/<NUM_LIT><EOL>dL_dKmm +=  dL_dpsi2R - backsub_both_sides(Lm, LmInvPsi2LmInvT)/<NUM_LIT:2><EOL>dL_dthetaL[d:d+<NUM_LIT:1>] = (YRY_d*beta_d + beta_d*psi0_d - N_d*beta_d)/<NUM_LIT> - beta_d*(dL_dpsi2R*psi2_d).sum() - beta_d*np.trace(LLinvPsi1TYYTPsi1LLinvT)<EOL>dL_dpsi0[idx_d] = -beta_d/<NUM_LIT><EOL>dL_dpsi1[idx_d] = beta_d*np.dot(Y_d,v)<EOL>dL_dpsi2[idx_d] = beta_d*dL_dpsi2R<EOL>wv[:,d] = v<EOL><DEDENT>LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_sum, '<STR_LIT:right>')<EOL>Lambda = np.eye(Kmm.shape[<NUM_LIT:0>])+LmInvPsi2LmInvT<EOL>LL = jitchol(Lambda)<EOL>LmLL = Lm.dot(LL)<EOL>logdet_L = <NUM_LIT>*np.sum(np.log(np.diag(LL)))<EOL>dL_dpsi2R_common = dpotri(LmLL)[<NUM_LIT:0>]/-<NUM_LIT><EOL>dL_dpsi2 += dL_dpsi2R_common[None,:,:]*beta_exp[:,None,None]<EOL>for d in range(output_dim):<EOL><INDENT>dL_dthetaL[d] += (dL_dpsi2R_common*psi2[indexD==d].sum(<NUM_LIT:0>)).sum()*-beta[d]*beta[d]<EOL><DEDENT>dL_dKmm += dL_dpsi2R_common*output_dim<EOL>logL += -output_dim*logdet_L/<NUM_LIT><EOL>LLInvLmT = dtrtrs(LL, Lm.T)[<NUM_LIT:0>]<EOL>cov = tdot(LLInvLmT.T)<EOL>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=wv, K=Kmm, mean=None, cov=cov, K_chol=Lm)<EOL>if not uncertain_inputs:<EOL><INDENT>dL_dpsi1 += (psi1[:,None,:]*dL_dpsi2).sum(<NUM_LIT:2>)*<NUM_LIT><EOL><DEDENT>if uncertain_inputs:<EOL><INDENT>grad_dict = {'<STR_LIT>': dL_dKmm,<EOL>'<STR_LIT>':dL_dpsi0,<EOL>'<STR_LIT>':dL_dpsi1,<EOL>'<STR_LIT>':dL_dpsi2,<EOL>'<STR_LIT>':dL_dthetaL}<EOL><DEDENT>else:<EOL><INDENT>grad_dict = {'<STR_LIT>': dL_dKmm,<EOL>'<STR_LIT>':dL_dpsi0,<EOL>'<STR_LIT>':dL_dpsi1,<EOL>'<STR_LIT>':dL_dthetaL}<EOL><DEDENT>return post, logL, grad_dict<EOL>
The first phase of inference:
Compute: log-likelihood, dL_dKmm

Cached intermediate results: Kmm, KmmInv,