code: def lt_factor(s, l, M, V, mp, p, gamma):
      cVc = (V[l, l] - 2 * V[s, l] + V[s, s]) / 2.0
      Vc = (V[:, l] - V[:, s]) / sq2
      cM = (M[l] - M[s]) / sq2
      cVnic = np.max([cVc / (1 - p * cVc), 0])
      cmni = cM + cVnic * (p * cM - mp)
      z = cmni / np.sqrt(cVnic + 1e-25)
      if np.isnan(z):
          z = -np.inf
      e, lP, exit_flag = log_relative_gauss(z)
      if exit_flag == 0:
          alpha = e / np.sqrt(cVnic)
          # beta  = alpha * (alpha + cmni / cVnic);
          # r     = beta * cVnic / (1 - cVnic * beta);
          beta = alpha * (alpha * cVnic + cmni)
          r = beta / (1 - beta)
          # new message
          pnew = r / cVnic
          mpnew = r * (alpha + cmni / cVnic) + alpha
  
          # update terms
          dp = np.max(
              [-p + eps, gamma * (pnew - p)]
          )  # at worst, remove message
          dmp = np.max([-mp + eps, gamma * (mpnew - mp)])
          d = np.max([dmp, dp])  # for convergence measures
  
          pnew = p + dp
          mpnew = mp + dmp
          # project out to marginal
          Vnew = V - dp / (1 + dp * cVc) * np.outer(Vc, Vc)
  
          Mnew = M + (dmp - cM * dp) / (1 + dp * cVc) * Vc
          if np.any(np.isnan(Vnew)):
              raise Exception(
                  "an error occurs while running expectation "
                  "propagation in entropy search. "
                  "Resulting variance contains NaN"
              )
          # % there is a problem here, when z is very large
          logS = (
              lP
              - 0.5 * (np.log(beta) - np.log(pnew) - np.log(cVnic))
              + (alpha * alpha) / (2 * beta) * cVnic
          )
  
      elif exit_flag == -1:
          d = np.NAN
          Mnew = 0
          Vnew = 0
          pnew = 0
          mpnew = 0
          logS = -np.Infinity
      elif exit_flag == 1:
          d = 0
          # remove message from marginal:
          # new message
          pnew = 0
          mpnew = 0
          # update terms
          dp = -p  # at worst, remove message
          dmp = -mp
          d = max([dmp, dp])  # for convergence measures
          # project out to marginal
          Vnew = V - dp / (1 + dp * cVc) * (np.outer(Vc, Vc))
          Mnew = M + (dmp - cM * dp) / (1 + dp * cVc) * Vc
          logS = 0
      return Mnew, Vnew, pnew, mpnew, logS, d
  
imports: import numpy
  from scipy import special
label_desc: Use image node_id 3 for calling a global function with example usage: lt_factor(s, l, M, V, mp, p, gamma) and returns: Mnew,  Vnew,  pnew,  mpnew,  logS,  d