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

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