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
andto_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