# Compute stats on two class vectors

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


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

• Similarly on mine (15.3 times faster). That's already a big benefit. Could you be more concrete on the second suggestion? I doubt that would be even faster, would it? Mar 31, 2014 at 17:29
• @embert My second suggestion wasn't that good after all, but see my third suggestion. Mar 31, 2014 at 20:16
• I don't yet fully grab why this works, but apparently it does a good job. Apr 3, 2014 at 6:55