def _adjusted_mutual_info_score(reference_indices, estimated_indices):
    """Compute the mutual information between two sequence labelings, adjusted for
    chance.

    Parameters
    ----------
    reference_indices : np.ndarray
        Array of reference indices

    estimated_indices : np.ndarray
        Array of estimated indices

    Returns
    -------
    ami : float <= 1.0
        Mutual information

    .. note:: Based on sklearn.metrics.cluster.adjusted_mutual_info_score
        and sklearn.metrics.cluster.expected_mutual_info_score

    """
    n_samples = len(reference_indices)
    ref_classes = np.unique(reference_indices)
    est_classes = np.unique(estimated_indices)
    # Special limit cases: no clustering since the data is not split.
    # This is a perfect match hence return 1.0.
    if (ref_classes.shape[0] == est_classes.shape[0] == 1 or
            ref_classes.shape[0] == est_classes.shape[0] == 0):
        return 1.0
    contingency = _contingency_matrix(reference_indices,
                                      estimated_indices).astype(float)
    # Calculate the MI for the two clusterings
    mi = _mutual_info_score(reference_indices, estimated_indices,
                            contingency=contingency)
    # The following code is based on
    # sklearn.metrics.cluster.expected_mutual_information
    R, C = contingency.shape
    N = float(n_samples)
    a = np.sum(contingency, axis=1).astype(np.int32)
    b = np.sum(contingency, axis=0).astype(np.int32)
    # There are three major terms to the EMI equation, which are multiplied to
    # and then summed over varying nij values.
    # While nijs[0] will never be used, having it simplifies the indexing.
    nijs = np.arange(0, max(np.max(a), np.max(b)) + 1, dtype='float')
    # Stops divide by zero warnings. As its not used, no issue.
    nijs[0] = 1
    # term1 is nij / N
    term1 = nijs / N
    # term2 is log((N*nij) / (a * b)) == log(N * nij) - log(a * b)
    # term2 uses the outer product
    log_ab_outer = np.log(np.outer(a, b))
    # term2 uses N * nij
    log_Nnij = np.log(N * nijs)
    # term3 is large, and involved many factorials. Calculate these in log
    # space to stop overflows.
    gln_a = scipy.special.gammaln(a + 1)
    gln_b = scipy.special.gammaln(b + 1)
    gln_Na = scipy.special.gammaln(N - a + 1)
    gln_Nb = scipy.special.gammaln(N - b + 1)
    gln_N = scipy.special.gammaln(N + 1)
    gln_nij = scipy.special.gammaln(nijs + 1)
    # start and end values for nij terms for each summation.
    start = np.array([[v - N + w for w in b] for v in a], dtype='int')
    start = np.maximum(start, 1)
    end = np.minimum(np.resize(a, (C, R)).T, np.resize(b, (R, C))) + 1
    # emi itself is a summation over the various values.
    emi = 0
    for i in range(R):
        for j in range(C):
            for nij in range(start[i, j], end[i, j]):
                term2 = log_Nnij[nij] - log_ab_outer[i, j]
                # Numerators are positive, denominators are negative.
                gln = (gln_a[i] + gln_b[j] + gln_Na[i] + gln_Nb[j] -
                       gln_N - gln_nij[nij] -
                       scipy.special.gammaln(a[i] - nij + 1) -
                       scipy.special.gammaln(b[j] - nij + 1) -
                       scipy.special.gammaln(N - a[i] - b[j] + nij + 1))
                term3 = np.exp(gln)
                emi += (term1[nij] * term2 * term3)
    # Calculate entropy for each labeling
    h_true, h_pred = _entropy(reference_indices), _entropy(estimated_indices)
    ami = (mi - emi) / (max(h_true, h_pred) - emi)
    return ami
Compute the mutual information between two sequence labelings, adjusted for
    chance.

    Parameters
    ----------
    reference_indices : np.ndarray
        Array of reference indices

    estimated_indices : np.ndarray
        Array of estimated indices

    Returns
    -------
    ami : float <= 1.0
        Mutual information

    .. note:: Based on sklearn.metrics.cluster.adjusted_mutual_info_score
        and sklearn.metrics.cluster.expected_mutual_info_score