code: def joint_min(mu, var, with_derivatives=False, **kwargs):
      """
      Computes the probability of every given point to be the minimum
      based on the EPMGP[1] algorithm.
      [1] Cunningham, P. Hennig, and S. Lacoste-Julien.
      Gaussian probabilities and expectation propagation.
      under review. Preprint at arXiv, November 2011.
      Parameters
      ----------
      M: np.ndarray(N,)
          Mean value of each of the N points.
      V: np.ndarray(N, N)
          Covariance matrix for all points
      with_derivatives: bool
          If true than also the gradients are computed
      Returns
      -------
      np.ndarray(N,1)
          pmin distribution
      """
  
      logP = np.zeros(mu.shape)
      D = mu.shape[0]
      if with_derivatives:
          dlogPdMu = np.zeros((D, D))
          dlogPdSigma = np.zeros((D, int(0.5 * D * (D + 1))))
          dlogPdMudMu = np.zeros((D, D, D))
      for i in range(mu.shape[0]):
          # logP[k] ) self._min_faktor(mu, var, 0)
          a = min_faktor(mu, var, i)
  
          logP[i] = next(a)
          if with_derivatives:
              dlogPdMu[i, :] = next(a).T
              dlogPdMudMu[i, :, :] = next(a)
              dlogPdSigma[i, :] = next(a).T
  
      logP[np.isinf(logP)] = -500
      # re-normalize at the end, to smooth out numerical imbalances:
      logPold = logP
      Z = np.sum(np.exp(logPold))
      maxLogP = np.max(logP)
      s = maxLogP + np.log(np.sum(np.exp(logP - maxLogP)))
      s = maxLogP if np.isinf(s) else s
  
      logP = logP - s
      if not with_derivatives:
          return logP
  
      dlogPdMuold = dlogPdMu
      dlogPdSigmaold = dlogPdSigma
      dlogPdMudMuold = dlogPdMudMu
      # adjust derivatives, too. This is a bit tedious.
      Zm = (
          sum(np.rot90((np.exp(logPold) * np.rot90(dlogPdMuold, 1)), 3))
          / Z
      )
      Zs = (
          sum(
              np.rot90(
                  (np.exp(logPold) * np.rot90(dlogPdSigmaold, 1)), 3
              )
          )
          / Z
      )
  
      dlogPdMu = dlogPdMuold - Zm
      dlogPdSigma = dlogPdSigmaold - Zs
  
      ff = np.einsum("ki,kj->kij", dlogPdMuold, dlogPdMuold)
      gg = (
          np.einsum("kij,k->ij", dlogPdMudMuold + ff, np.exp(logPold))
          / Z
      )
      Zij = Zm.T * Zm
      adds = np.reshape(-gg + Zij, (1, D, D))
      dlogPdMudMu = dlogPdMudMuold + adds
      return logP, dlogPdMu, dlogPdSigma, dlogPdMudMu
  
imports: import numpy
  from scipy import special
label_desc: Use image node_id 1 for calling a global function with example usage: joint_min(mu, var, with_derivatives) and returns: logP,  dlogPdMu,  dlogPdSigma,  dlogPdMudMu, logP