def rasm_mode(self, K, Y, likelihood, Ki_f_init, Y_metadata=None, *args, **kwargs):
Ki_f = Ki_f_init.copy()<EOL>f = np.dot(K, Ki_f)<EOL>def obj(Ki_f, f):<EOL><INDENT>ll = -<NUM_LIT:0.5>*np.sum(np.dot(Ki_f.T, f)) + np.sum(likelihood.logpdf(f, Y, Y_metadata=Y_metadata))<EOL>if np.isnan(ll):<EOL><INDENT>import ipdb; ipdb.set_trace()  <EOL>return -np.inf<EOL><DEDENT>else:<EOL><INDENT>return ll<EOL><DEDENT><DEDENT>difference = np.inf<EOL>iteration = <NUM_LIT:0><EOL>while difference > self._mode_finding_tolerance and iteration < self._mode_finding_max_iter:<EOL><INDENT>W = -likelihood.d2logpdf_df2(f, Y, Y_metadata=Y_metadata)<EOL>if np.any(np.isnan(W)):<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>grad = likelihood.dlogpdf_df(f, Y, Y_metadata=Y_metadata)<EOL>if np.any(np.isnan(grad)):<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>W_f = W*f<EOL>b = W_f + grad <EOL>W12BiW12, _, _, _ = self._compute_B_statistics(K, W, likelihood.log_concave, *args, **kwargs)<EOL>W12BiW12Kb = np.dot(W12BiW12, np.dot(K, b))<EOL>full_step_Ki_f = b - W12BiW12Kb <EOL>dKi_f = full_step_Ki_f - Ki_f<EOL>def inner_obj(step_size):<EOL><INDENT>Ki_f_trial = Ki_f + step_size*dKi_f<EOL>f_trial = np.dot(K, Ki_f_trial)<EOL>return -obj(Ki_f_trial, f_trial)<EOL><DEDENT>step = optimize.brent(inner_obj, tol=<NUM_LIT>, maxiter=<NUM_LIT:12>)<EOL>Ki_f_new = Ki_f + step*dKi_f<EOL>f_new = np.dot(K, Ki_f_new)<EOL>old_obj = obj(Ki_f, f)<EOL>new_obj = obj(Ki_f_new, f_new)<EOL>if new_obj < old_obj:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>difference = np.abs(new_obj - old_obj)<EOL>Ki_f = Ki_f_new<EOL>f = f_new<EOL>iteration += <NUM_LIT:1><EOL><DEDENT>if difference > self._mode_finding_tolerance:<EOL><INDENT>if not self.bad_fhat:<EOL><INDENT>warnings.warn("<STR_LIT>".format(difference, iteration, self._mode_finding_max_iter))<EOL><DEDENT>self.bad_fhat = True<EOL><DEDENT>elif self.bad_fhat:<EOL><INDENT>self.bad_fhat = False<EOL>warnings.warn("<STR_LIT>".format(difference, iteration, self._mode_finding_max_iter))<EOL><DEDENT>return f, Ki_f<EOL>
Rasmussen's numerically stable mode finding
For nomenclature see Rasmussen & Williams 2006
Influenced by GPML (BSD) code, all errors are our own

:param K: Covariance matrix evaluated at locations X
:type K: NxD matrix
:param Y: The data
:type Y: np.ndarray
:param likelihood: the likelihood of the latent function value for the given data
:type likelihood: a GPy.likelihood object
:param Ki_f_init: the initial guess at the mode
:type Ki_f_init: np.ndarray
:param Y_metadata: information about the data, e.g. which likelihood to take from a multi-likelihood object
:type Y_metadata: np.ndarray | None
:returns: f_hat, mode on which to make laplace approxmiation
:rtype: np.ndarray