def cross_entropy_neighbors_in_rep(adata, use_rep, n_points=3):
    """Compare neighborhood graph of representation based on cross entropy.

    `n_points` denotes the number of points to add as highlight annotation.

    Returns
    -------
    The cross entropy and the geodesic-distance-weighted cross entropy as
    ``entropy, geo_entropy_d, geo_entropy_o``.

    Adds the most overlapping or disconnected points as annotation to `adata`.
    """
    # see below why we need this
    if 'X_diffmap' not in adata.obsm.keys():
        raise ValueError('Run `tl.diffmap` on `adata`, first.')

    adata_ref = adata  # simple renaming, don't need copy here
    adata_cmp = adata.copy()
    n_neighbors = adata_ref.uns['neighbors']['params']['n_neighbors']
    from .neighbors import neighbors
    neighbors(adata_cmp, n_neighbors=n_neighbors, use_rep=use_rep)
    from .tools.diffmap import diffmap
    diffmap(adata_cmp)

    graph_ref = adata_ref.uns['neighbors']['connectivities']
    graph_cmp = adata_cmp.uns['neighbors']['connectivities']

    graph_ref = graph_ref.tocoo()  # makes a copy
    graph_cmp = graph_cmp.tocoo()

    edgeset_ref = {e for e in zip(graph_ref.row, graph_ref.col)}
    edgeset_cmp = {e for e in zip(graph_cmp.row, graph_cmp.col)}
    edgeset_union = list(edgeset_ref.union(edgeset_cmp))

    edgeset_union_indices = tuple(zip(*edgeset_union))
    edgeset_union_indices = (np.array(edgeset_union_indices[0]), np.array(edgeset_union_indices[1]))

    n_edges_ref = len(graph_ref.nonzero()[0])
    n_edges_cmp = len(graph_cmp.nonzero()[0])
    n_edges_union = len(edgeset_union)
    logg.msg(
        '... n_edges_ref', n_edges_ref,
        'n_edges_cmp', n_edges_cmp,
        'n_edges_union', n_edges_union)

    graph_ref = graph_ref.tocsr()  # need a copy of the csr graph anyways
    graph_cmp = graph_cmp.tocsr()

    p_ref = graph_ref[edgeset_union_indices].A1
    p_cmp = graph_cmp[edgeset_union_indices].A1

    # the following is how one compares it to log_loss form sklearn
    # p_ref[p_ref.nonzero()] = 1
    # from sklearn.metrics import log_loss
    # print(log_loss(p_ref, p_cmp))
    p_cmp = np.clip(p_cmp, EPS, 1-EPS)
    ratio = np.clip(p_ref / p_cmp, EPS, None)
    ratio_1m = np.clip((1 - p_ref) / (1 - p_cmp), EPS, None)

    entropy = np.sum(p_ref * np.log(ratio) + (1-p_ref) * np.log(ratio_1m))

    n_edges_fully_connected = (graph_ref.shape[0]**2 - graph_ref.shape[0])
    entropy /= n_edges_fully_connected

    fraction_edges = n_edges_ref / n_edges_fully_connected
    naive_entropy = (fraction_edges * np.log(1./fraction_edges)
                     + (1-fraction_edges) * np.log(1./(1-fraction_edges)))
    logg.msg('cross entropy of naive sparse prediction {:.3e}'.format(naive_entropy))
    logg.msg('cross entropy of random prediction {:.3e}'.format(-np.log(0.5)))
    logg.info('cross entropy {:.3e}'.format(entropy))

    # for manifold analysis, restrict to largest connected component in
    # reference
    # now that we clip at a quite high value below, this might not even be
    # necessary
    n_components, labels = scipy.sparse.csgraph.connected_components(graph_ref)
    largest_component = np.arange(graph_ref.shape[0], dtype=int)
    if n_components > 1:
        component_sizes = np.bincount(labels)
        logg.msg('largest component has size', component_sizes.max())
        largest_component = np.where(
            component_sizes == component_sizes.max())[0][0]
        graph_ref_red = graph_ref.tocsr()[labels == largest_component, :]
        graph_ref_red = graph_ref_red.tocsc()[:, labels == largest_component]
        graph_ref_red = graph_ref_red.tocoo()
        graph_cmp_red = graph_cmp.tocsr()[labels == largest_component, :]
        graph_cmp_red = graph_cmp_red.tocsc()[:, labels == largest_component]
        graph_cmp_red = graph_cmp_red.tocoo()
        edgeset_ref_red = {e for e in zip(graph_ref_red.row, graph_ref_red.col)}
        edgeset_cmp_red = {e for e in zip(graph_cmp_red.row, graph_cmp_red.col)}
        edgeset_union_red = edgeset_ref_red.union(edgeset_cmp_red)
        map_indices = np.where(labels == largest_component)[0]
        edgeset_union_red = {
            (map_indices[i], map_indices[j]) for (i, j) in edgeset_union_red}

    from .neighbors import Neighbors
    neigh_ref = Neighbors(adata_ref)
    dist_ref = neigh_ref.distances_dpt  # we expect 'X_diffmap' to be already present

    neigh_cmp = Neighbors(adata_cmp)
    dist_cmp = neigh_cmp.distances_dpt

    d_cmp = np.zeros_like(p_ref)
    d_ref = np.zeros_like(p_ref)
    for i, e in enumerate(edgeset_union):
        # skip contributions that are not in the largest component
        if n_components > 1 and e not in edgeset_union_red:
            continue
        d_cmp[i] = dist_cmp[e]
        d_ref[i] = dist_ref[e]

    MAX_DIST = 1000
    d_cmp = np.clip(d_cmp, 0.1, MAX_DIST)  # we don't want to measure collapsing clusters
    d_ref = np.clip(d_ref, 0.1, MAX_DIST)

    weights = np.array(d_cmp / d_ref)            # disconnected regions
    weights_overlap = np.array(d_ref / d_cmp)    # overlapping regions

    # the following is just for annotation of figures
    if 'highlights' not in adata_ref.uns:
        adata_ref.uns['highlights'] = {}
    else:
        # remove old disconnected and overlapping points
        new_highlights = {}
        for k, v in adata_ref.uns['highlights'].items():
            if v != 'O' and v not in {'D0', 'D1', 'D2', 'D3', 'D4'}:
                new_highlights[k] = v
        adata_ref.uns['highlights'] = new_highlights

    # points that are maximally disconnected
    max_weights = np.argpartition(weights, kth=-n_points)[-n_points:]
    points = list(edgeset_union_indices[0][max_weights])
    points2 = list(edgeset_union_indices[1][max_weights])
    found_disconnected_points = False
    for ip, p in enumerate(points):
        if d_cmp[max_weights][ip] == MAX_DIST:
            adata_ref.uns['highlights'][p] = 'D' + str(ip)
            adata_ref.uns['highlights'][points2[ip]] = 'D' + str(ip)
            found_disconnected_points = True
    if found_disconnected_points:
        logg.msg('most disconnected points', points)
        logg.msg('    with weights', weights[max_weights].round(1))

    max_weights = np.argpartition(
        weights_overlap, kth=-n_points)[-n_points:]
    points = list(edgeset_union_indices[0][max_weights])
    for p in points:
        adata_ref.uns['highlights'][p] = 'O'
    logg.msg('most overlapping points', points)
    logg.msg('    with weights', weights_overlap[max_weights].round(1))
    logg.msg('    with d_rep', d_cmp[max_weights].round(1))
    logg.msg('    with d_ref', d_ref[max_weights].round(1))

    geo_entropy_d = np.sum(weights * p_ref * np.log(ratio))
    geo_entropy_o = np.sum(weights_overlap * (1-p_ref) * np.log(ratio_1m))

    geo_entropy_d /= n_edges_fully_connected
    geo_entropy_o /= n_edges_fully_connected

    logg.info('geodesic cross entropy {:.3e}'.format(geo_entropy_d + geo_entropy_o))
    return entropy, geo_entropy_d, geo_entropy_o
Compare neighborhood graph of representation based on cross entropy.

    `n_points` denotes the number of points to add as highlight annotation.

    Returns
    -------
    The cross entropy and the geodesic-distance-weighted cross entropy as
    ``entropy, geo_entropy_d, geo_entropy_o``.

    Adds the most overlapping or disconnected points as annotation to `adata`.