Compute stats on two class vectors

Question

Below is a full working code example of code which is used to compute stats on two class vectors (partitions).

There are two functions:

• pairwise_indication and
• to_clust_set.

The first one is the function in question and the second is just for completeness, since it's apparently not so time consuming. pairwise_indication returns the indicators n11, n00, n10, n01 which are used to compute scores such as the RAND Index (they're named a, b, c, d on that site).

The bottleneck probably are the loops (depth is 3) and the logic and look-up stuff inside of pairwise_indication. I already tried to improve the logical function in it by pulling forward the common cases. I was hoping it can be improved some further. What makes it interminable (if I get the profiling at the bottom of the post correct) are the loops and maybe the set look-ups using in, though I cannot imagine something faster.

import itertools as it
import numpy as np

def pairwise_indication(part1, part2):
    r"""
    Computes the pairwise indicators.

    Parameters
    ----------
    part1, part2 : array, (n)
        The two partition vectors

    Returns
    -------
    n11, n00, n10, n01 : tuple (int)
        Tuple with the counts of incidences

    Examples
    --------
    >>> from cluvap.calc import pairwise_indication
    >>> import numpy as np
    >>> p1 = np.array([1, 2, 3, 3, 1, 2])
    >>> p2 = np.array([3, 2, 3, 3, 1, 2])
    >>> pairwise_indication(p1, p2)
    (2, 10, 1, 2)
    """
    n = len(part1)
    if n != len(part2):
        raise ValueError('Partition shapes do not match')
    # Create clustering as set
    A = to_clust_set(part1)
    B = to_clust_set(part2)
    observations = np.arange(n)
    n11, n00, n10, n01 = 0, 0, 0, 0
    # Count incidences
    # function:  P'QR'S v P'QRS v PQRS v PQR'S (' == not, v == or)
    # condensed: P('QR('S v S) v QR('S v S))
    # P: obs1 in a
    # Q: obs2 in a
    # R: obs1 in b
    # S: obs2 in b
    for obs1, obs2 in it.combinations(observations, 2):
        for a in A:
            for b in B:
                if obs1 in a:
                    if obs2 not in a:
                        if obs1 in b:
                            if obs2 not in b:
                                n00 += 1
                            else:
                                n01 += 1
                    else:
                        if obs1 in b:
                            if obs2 in b:
                                n11 += 1
                            else:
                                n10 += 1

    return n11, n00, n10, n01


def to_clust_set(part):
    """
    Converts a partition to a set of clusters. Noise is considered a cluster.

    Parameters
    ----------
    part : arrays, (n)
        Partitions. A partition ``p`` assigns the ``n``-th observation
        to the cluster ``p[n]``.

    Returns
    -------
    cset : list of sets

    Examples
    --------
    >>> from cluvap.calc import to_clust_set
    >>> import numpy as np
    >>> p = np.array([1, 2, 3, 3, 1, 2])
    >>> to_clust_set(p)
    [set([0, 4]), set([1, 5]), set([2, 3])]
    """
    clusters = set(part)
    obs = np.arange(len(part))

    return [set(obs[part == C]) for C in clusters]


if __name__ == "__main__":
    from time import time
    p1 = np.array([1, 2, 3, 4, 5] * 20)
    p2 = np.array([2, 3, 4, 5, 1] * 20)

    t0 = time()
    pi = pairwise_indication(p1, p2)
    t = (time() - t0) * 1000
    print 'pairwise_indication(p1, p2) in {:.2f} ms:\n'.format(t), pi

Output

pairwise_indication(p1, p2) in 19.00 ms:
(950, 4000, 0, 0)

line profiling

>>> p1 = np.array([1, 2, 3, 4, 5] * 20)
>>> p2 = np.array([2, 3, 4, 5, 1] * 20)
>>> %lprun -f pairwise_indication pairwise_indication(p1, p2)
Timer unit: 4.10918e-07 s

File: cluvap\calc.py
Function: pairwise_indication at line 713
Total time: 0.603277 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   713                                           def pairwise_indication(part1, part2):
...
   736         1           28     28.0      0.0      n = len(part1)
   737         1           17     17.0      0.0      if n != len(part2):
   738                                                   raise ValueError('Partition shapes do not match')
   739                                               # Create clustering as set
   740         1         1455   1455.0      0.1      A = to_clust_set(part1)
   741         1         1154   1154.0      0.1      B = to_clust_set(part2)
   742         1           33     33.0      0.0      observations = np.arange(n)
   743         1           14     14.0      0.0      n11, n00, n10, n01 = 0, 0, 0, 0
   744                                               # Count incidences
   745                                               # (can this be improved using some kind of indexing / logical compression?)
   746      4951        19884      4.0      1.4      for obs1, obs2 in it.combinations(observations, 2):
   747     29700       115151      3.9      7.8          for a in A:
   748    148500       576491      3.9     39.3              for b in B:
   749    123750       505606      4.1     34.4                  if obs1 in a:
   750     24750       100784      4.1      6.9                      if obs2 not in a:
   751     20000        84631      4.2      5.8                          if obs1 in b:
   752      4000        16489      4.1      1.1                              if obs2 not in b:
   753      4000        17915      4.5      1.2                                  n00 += 1
   754                                                                       else:
   755                                                                           n01 += 1
   756                                                               else:
   757      4750        20272      4.3      1.4                          if obs1 in b:
   758       950         4005      4.2      0.3                              if obs2 in b:
   759       950         4185      4.4      0.3                                  n11 += 1
   760                                                                       else:
   761                                                                           n10 += 1
   762
   763         1            5      5.0      0.0      return n11, n00, n10, n01

Answer 1

If I understood what you want to do, this would be a more direct way to compute the same result. At least for the test cases provided, the result is indeed the same.

def pairwise_indication(part1, part2):
    if len(part1) != len(part2):
        raise ValueError('Partition shapes do not match')
    n11, n00, n10, n01 = 0, 0, 0, 0
    for (a1,a2), (b1,b2) in it.combinations(zip(part1, part2), 2):
        if a1 == b1:
            if a2 == b2:
                n11 += 1
            else:
                n10 += 1
        else:
            if a2 == b2:
                n01 += 1
            else:
                n00 += 1
    return n11, n00, n10, n01

On my computer this is 13 times faster than yours.

---

Here's a completely different approach based on the idea that you can partition the partitions with one another (using zip in Python) to obtain a finer partition. (There must be a word for that?) From the size of each subset you can directly calculate the number of pairs that can be formed. Add up the numbers for each partition and subtract the overlap. This is 70 times faster than yours with the given example, and more importantly, operates in linear time, so it scales well to larger data.

from collections import Counter

def pairs(n):
    '''Calculate number of pairs that can be formed from n items'''
    return n * (n - 1) // 2

def partition_pairs(partition):
    '''Calculate number of pairs in subsets of partition'''
    return sum(pairs(x) for x in Counter(partition).values())

def pairwise_indication(part1, part2):
    n = len(part1)
    if n != len(part2):
        raise ValueError('Partition shapes do not match')

    n11 = partition_pairs(zip(part1, part2))
    n10 = partition_pairs(part1) - n11
    n01 = partition_pairs(part2) - n11
    n00 = pairs(n) - n11 - n10 - n01

    return n11, n00, n10, n01