Christian
Christian

Reputation: 577

How to compute fast the containment order of all intersections of a collection of sets

This is a follow-up to How to get all intersections of sets in python fast:

I have a finite collection A = {A1,...Ak} of finite sets Ai of integers and I want to compute in Python the following:

  1. All intersections of subsets of A: F = { intersection of B : B is a subset of A }. This was the above question with a decently fast solution.

  2. a. All pairs (X,Y) for X,Y being sets in F such that X is subset of Y.

    b. All pairs (X,Y) for X,Y being sets in F such that X is subset of Y and there is no set Z in F such that X subset of Z subset of Y. In other words, so that no set Z fits in between X and Y in containment order. Such a pair (X,Y) is called cover.

Why do I want to do that? -- I want to compute face lattices of 10^7 polytopes. In the scenario in mind, the collection A above contains 600 sets. It is indeed the famous 600-cell, the computation currently takes about 6 secs, and I would like that to go down by a factor of 10, if possible.

The 6 secs to get 2.a. is simply done by doing

# this is John Coleman's function from above question's answer
def allIntersections(frozenSets):
    universalSet = frozenset.union(*frozenSets)
    intersections = set([universalSet])
    for s in frozenSets:
        moreIntersections = set(s & t for t in intersections)
        intersections.update(moreIntersections)
    return intersections

def all_intersections(lists):
    sets = allIntersections([frozenset(s) for s in lists])
    return [list(s) for s in sets]


A = [[19, 40, 41, 48], [19, 44, 45, 49], [23, 42, 43, 50], [23, 46, 47, 51], [19, 40, 41, 52], [19, 44, 45, 53], [23, 42, 43, 54], [23, 46, 47, 55], [2, 25, 36, 56], [0, 24, 32, 56], [24, 25, 56, 57], [24, 32, 56, 57], [16, 32, 56, 57], [1, 24, 32, 57], [25, 36, 56, 57], [16, 36, 56, 57], [3, 25, 36, 57], [8, 28, 34, 58], [10, 29, 38, 58], [28, 29, 58, 59], [28, 34, 58, 59], [20, 34, 58, 59], [29, 38, 58, 59], [20, 38, 58, 59], [9, 28, 34, 59], [11, 29, 38, 59], [6, 27, 37, 60], [4, 26, 33, 60], [5, 26, 33, 61], [26, 27, 60, 61], [26, 33, 60, 61], [16, 33, 60, 61], [27, 37, 60, 61], [7, 27, 37, 61], [16, 37, 60, 61], [12, 30, 35, 62], [14, 31, 39, 62], [30, 35, 62, 63], [20, 39, 62, 63], [20, 35, 62, 63], [30, 31, 62, 63], [31, 39, 62, 63], [15, 31, 39, 63], [13, 30, 35, 63], [0, 24, 32, 64], [1, 24, 32, 64], [8, 28, 34, 65], [9, 28, 34, 65], [3, 25, 36, 66], [2, 25, 36, 66], [11, 29, 38, 67], [10, 29, 38, 67], [4, 26, 33, 68], [5, 26, 33, 68], [12, 30, 35, 69], [13, 30, 35, 69], [6, 27, 37, 70], [7, 27, 37, 70], [15, 31, 39, 71], [14, 31, 39, 71], [4, 33, 68, 72], [0, 32, 64, 72], [18, 64, 72, 73], [32, 64, 72, 73], [32, 33, 72, 73], [1, 32, 64, 73], [18, 68, 72, 73], [5, 33, 68, 73], [33, 68, 72, 73], [2, 36, 66, 74], [6, 37, 70, 74], [3, 36, 66, 75], [7, 37, 70, 75], [36, 66, 74, 75], [37, 70, 74, 75], [36, 37, 74, 75], [22, 66, 74, 75], [22, 70, 74, 75], [12, 35, 69, 76], [8, 34, 65, 76], [18, 65, 76, 77], [34, 65, 76, 77], [34, 35, 76, 77], [18, 69, 76, 77], [35, 69, 76, 77], [13, 35, 69, 77], [9, 34, 65, 77], [10, 38, 67, 78], [14, 39, 71, 78], [38, 67, 78, 79], [22, 71, 78, 79], [22, 67, 78, 79], [38, 39, 78, 79], [39, 71, 78, 79], [15, 39, 71, 79], [11, 38, 67, 79], [0, 40, 48, 80], [19, 40, 48, 80], [19, 48, 49, 80], [8, 44, 49, 80], [19, 44, 49, 80], [2, 40, 52, 81], [10, 44, 53, 81], [19, 52, 53, 81], [19, 40, 52, 81], [19, 44, 53, 81], [19, 40, 80, 81], [19, 44, 80, 81], [23, 42, 50, 82], [23, 50, 51, 82], [1, 42, 50, 82], [23, 46, 51, 82], [9, 46, 51, 82], [23, 54, 55, 83], [3, 42, 54, 83], [23, 42, 54, 83], [23, 42, 82, 83], [11, 46, 55, 83], [23, 46, 55, 83], [23, 46, 82, 83], [19, 45, 49, 84], [12, 45, 49, 84], [4, 41, 48, 84], [19, 41, 48, 84], [19, 48, 49, 84], [19, 45, 84, 85], [19, 41, 84, 85], [6, 41, 52, 85], [19, 41, 52, 85], [14, 45, 53, 85], [19, 45, 53, 85], [19, 52, 53, 85], [23, 43, 50, 86], [5, 43, 50, 86], [23, 50, 51, 86], [23, 47, 51, 86], [13, 47, 51, 86], [7, 43, 54, 87], [23, 43, 54, 87], [23, 43, 86, 87], [23, 54, 55, 87], [23, 47, 86, 87], [15, 47, 55, 87], [23, 47, 55, 87], [8, 28, 65, 88], [0, 24, 64, 88], [9, 28, 65, 89], [28, 65, 88, 89], [17, 28, 88, 89], [17, 24, 88, 89], [1, 24, 64, 89], [24, 64, 88, 89], [64, 65, 88, 89], [4, 26, 68, 90], [12, 30, 69, 90], [5, 26, 68, 91], [13, 30, 69, 91], [26, 68, 90, 91], [21, 26, 90, 91], [68, 69, 90, 91], [30, 69, 90, 91], [21, 30, 90, 91], [10, 29, 67, 92], [2, 25, 66, 92], [29, 67, 92, 93], [66, 67, 92, 93], [11, 29, 67, 93], [17, 29, 92, 93], [25, 66, 92, 93], [17, 25, 92, 93], [3, 25, 66, 93], [14, 31, 71, 94], [6, 27, 70, 94], [21, 31, 94, 95], [21, 27, 94, 95], [15, 31, 71, 95], [31, 71, 94, 95], [70, 71, 94, 95], [27, 70, 94, 95], [7, 27, 70, 95], [2, 25, 56, 96], [0, 80, 88, 96], [0, 40, 56, 96], [2, 40, 81, 96], [2, 40, 56, 96], [0, 40, 80, 96], [40, 80, 81, 96], [2, 81, 92, 96], [17, 25, 92, 96], [2, 25, 92, 96], [0, 24, 88, 96], [0, 24, 56, 96], [24, 25, 56, 96], [17, 24, 88, 96], [17, 24, 25, 96], [28, 29, 58, 97], [80, 88, 96, 97], [80, 81, 96, 97], [44, 80, 81, 97], [8, 28, 88, 97], [8, 28, 58, 97], [8, 44, 58, 97], [8, 80, 88, 97], [8, 44, 80, 97], [81, 92, 96, 97], [17, 29, 92, 97], [17, 92, 96, 97], [17, 28, 29, 97], [17, 28, 88, 97], [17, 88, 96, 97], [10, 29, 92, 97], [10, 29, 58, 97], [10, 44, 58, 97], [10, 44, 81, 97], [10, 81, 92, 97], [6, 41, 85, 98], [6, 41, 60, 98], [4, 41, 60, 98], [6, 85, 94, 98], [4, 41, 84, 98], [4, 84, 90, 98], [41, 84, 85, 98], [6, 27, 94, 98], [6, 27, 60, 98], [26, 27, 60, 98], [4, 26, 90, 98], [4, 26, 60, 98], [21, 27, 94, 98], [21, 26, 90, 98], [21, 26, 27, 98], [14, 45, 85, 99], [21, 30, 31, 99], [14, 31, 62, 99], [30, 31, 62, 99], [14, 45, 62, 99], [21, 90, 98, 99], [21, 30, 90, 99], [84, 90, 98, 99], [45, 84, 85, 99], [84, 85, 98, 99], [12, 30, 62, 99], [12, 45, 62, 99], [12, 45, 84, 99], [12, 30, 90, 99], [12, 84, 90, 99], [85, 94, 98, 99], [21, 94, 98, 99], [14, 85, 94, 99], [14, 31, 94, 99], [21, 31, 94, 99], [3, 83, 93, 100], [1, 42, 82, 100], [3, 42, 57, 100], [1, 42, 57, 100], [42, 82, 83, 100], [3, 42, 83, 100], [1, 82, 89, 100], [1, 24, 89, 100], [17, 24, 89, 100], [1, 24, 57, 100], [17, 25, 93, 100], [3, 25, 57, 100], [3, 25, 93, 100], [17, 24, 25, 100], [24, 25, 57, 100], [17, 93, 100, 101], [82, 83, 100, 101], [11, 83, 93, 101], [83, 93, 100, 101], [11, 29, 59, 101], [11, 29, 93, 101], [17, 29, 93, 101], [9, 82, 89, 101], [82, 89, 100, 101], [17, 89, 100, 101], [11, 46, 83, 101], [11, 46, 59, 101], [9, 46, 59, 101], [9, 46, 82, 101], [46, 82, 83, 101], [9, 28, 59, 101], [17, 28, 29, 101], [28, 29, 59, 101], [17, 28, 89, 101], [9, 28, 89, 101], [5, 43, 86, 102], [5, 86, 91, 102], [7, 43, 61, 102], [5, 43, 61, 102], [21, 27, 95, 102], [7, 27, 95, 102], [7, 27, 61, 102], [5, 26, 61, 102], [26, 27, 61, 102], [21, 26, 27, 102], [21, 26, 91, 102], [5, 26, 91, 102], [43, 86, 87, 102], [7, 43, 87, 102], [7, 87, 95, 102], [86, 91, 102, 103], [86, 87, 102, 103], [15, 31, 63, 103], [30, 31, 63, 103], [15, 31, 95, 103], [87, 95, 102, 103], [15, 87, 95, 103], [15, 47, 63, 103], [15, 47, 87, 103], [47, 86, 87, 103], [13, 30, 63, 103], [13, 30, 91, 103], [13, 86, 91, 103], [13, 47, 63, 103], [13, 47, 86, 103], [21, 91, 102, 103], [21, 30, 91, 103], [21, 30, 31, 103], [21, 95, 102, 103], [21, 31, 95, 103], [0, 48, 72, 104], [4, 33, 72, 104], [4, 33, 60, 104], [4, 41, 60, 104], [4, 48, 72, 104], [4, 41, 48, 104], [32, 33, 72, 104], [0, 32, 72, 104], [0, 32, 56, 104], [0, 40, 56, 104], [40, 41, 48, 104], [0, 40, 48, 104], [16, 32, 56, 104], [16, 32, 33, 104], [16, 33, 60, 104], [40, 41, 104, 105], [40, 41, 52, 105], [41, 60, 104, 105], [16, 60, 104, 105], [40, 56, 104, 105], [16, 56, 104, 105], [2, 40, 56, 105], [2, 40, 52, 105], [2, 36, 56, 105], [16, 36, 56, 105], [16, 37, 60, 105], [16, 36, 37, 105], [2, 52, 74, 105], [36, 37, 74, 105], [2, 36, 74, 105], [6, 52, 74, 105], [6, 41, 52, 105], [6, 41, 60, 105], [6, 37, 60, 105], [6, 37, 74, 105], [12, 35, 76, 106], [12, 45, 62, 106], [12, 35, 62, 106], [8, 44, 49, 106], [8, 49, 76, 106], [12, 49, 76, 106], [44, 45, 49, 106], [12, 45, 49, 106], [20, 35, 62, 106], [8, 44, 58, 106], [20, 34, 58, 106], [8, 34, 58, 106], [20, 34, 35, 106], [8, 34, 76, 106], [34, 35, 76, 106], [20, 62, 106, 107], [20, 38, 39, 107], [20, 39, 62, 107], [10, 38, 78, 107], [38, 39, 78, 107], [10, 53, 78, 107], [20, 58, 106, 107], [20, 38, 58, 107], [10, 38, 58, 107], [44, 58, 106, 107], [10, 44, 58, 107], [10, 44, 53, 107], [14, 39, 62, 107], [14, 39, 78, 107], [14, 53, 78, 107], [14, 45, 53, 107], [44, 45, 106, 107], [44, 45, 53, 107], [14, 45, 62, 107], [45, 62, 106, 107], [16, 32, 57, 108], [1, 32, 57, 108], [16, 32, 33, 108], [16, 33, 61, 108], [5, 33, 61, 108], [1, 32, 73, 108], [32, 33, 73, 108], [1, 50, 73, 108], [5, 33, 73, 108], [5, 50, 73, 108], [1, 42, 50, 108], [1, 42, 57, 108], [5, 43, 61, 108], [5, 43, 50, 108], [42, 43, 50, 108], [7, 37, 61, 109], [3, 36, 57, 109], [3, 42, 57, 109], [7, 43, 61, 109], [42, 43, 108, 109], [43, 61, 108, 109], [42, 57, 108, 109], [16, 36, 57, 109], [16, 36, 37, 109], [16, 57, 108, 109], [16, 61, 108, 109], [16, 37, 61, 109], [36, 37, 75, 109], [7, 37, 75, 109], [3, 36, 75, 109], [3, 42, 54, 109], [42, 43, 54, 109], [7, 43, 54, 109], [3, 54, 75, 109], [7, 54, 75, 109], [34, 35, 77, 110], [13, 35, 63, 110], [13, 35, 77, 110], [13, 47, 63, 110], [9, 34, 77, 110], [9, 51, 77, 110], [13, 51, 77, 110], [9, 46, 51, 110], [13, 47, 51, 110], [46, 47, 51, 110], [20, 35, 63, 110], [20, 34, 35, 110], [9, 34, 59, 110], [20, 34, 59, 110], [9, 46, 59, 110], [11, 38, 59, 111], [11, 38, 79, 111], [15, 47, 63, 111], [11, 55, 79, 111], [15, 47, 55, 111], [15, 55, 79, 111], [11, 46, 59, 111], [46, 47, 55, 111], [11, 46, 55, 111], [38, 39, 79, 111], [15, 39, 79, 111], [15, 39, 63, 111], [20, 38, 39, 111], [20, 39, 63, 111], [20, 38, 59, 111], [47, 63, 110, 111], [20, 59, 110, 111], [20, 63, 110, 111], [46, 59, 110, 111], [46, 47, 110, 111], [8, 65, 88, 112], [18, 65, 76, 112], [8, 65, 76, 112], [8, 49, 76, 112], [0, 64, 88, 112], [64, 65, 88, 112], [18, 64, 65, 112], [18, 64, 72, 112], [0, 64, 72, 112], [0, 48, 72, 112], [8, 49, 80, 112], [8, 80, 88, 112], [48, 49, 80, 112], [0, 48, 80, 112], [0, 80, 88, 112], [4, 68, 90, 113], [4, 84, 90, 113], [12, 84, 90, 113], [18, 68, 69, 113], [68, 69, 90, 113], [12, 69, 90, 113], [4, 68, 72, 113], [18, 68, 72, 113], [18, 69, 76, 113], [12, 69, 76, 113], [4, 48, 84, 113], [4, 48, 72, 113], [12, 49, 76, 113], [48, 49, 84, 113], [12, 49, 84, 113], [18, 76, 112, 113], [49, 76, 112, 113], [18, 72, 112, 113], [48, 49, 112, 113], [48, 72, 112, 113], [2, 66, 92, 114], [66, 67, 92, 114], [52, 53, 81, 114], [2, 52, 81, 114], [2, 81, 92, 114], [22, 66, 67, 114], [22, 67, 78, 114], [2, 66, 74, 114], [2, 52, 74, 114], [22, 66, 74, 114], [10, 53, 81, 114], [10, 53, 78, 114], [10, 67, 78, 114], [10, 81, 92, 114], [10, 67, 92, 114], [6, 85, 94, 115], [6, 52, 85, 115], [52, 53, 85, 115], [14, 85, 94, 115], [14, 53, 85, 115], [52, 53, 114, 115], [6, 52, 74, 115], [52, 74, 114, 115], [14, 71, 94, 115], [70, 71, 94, 115], [6, 70, 94, 115], [6, 70, 74, 115], [22, 74, 114, 115], [22, 70, 74, 115], [22, 70, 71, 115], [14, 71, 78, 115], [53, 78, 114, 115], [14, 53, 78, 115], [22, 78, 114, 115], [22, 71, 78, 115], [18, 64, 65, 116], [18, 65, 77, 116], [50, 51, 82, 116], [1, 50, 82, 116], [9, 51, 82, 116], [9, 51, 77, 116], [9, 65, 77, 116], [18, 64, 73, 116], [1, 64, 73, 116], [1, 50, 73, 116], [64, 65, 89, 116], [1, 64, 89, 116], [9, 65, 89, 116], [1, 82, 89, 116], [9, 82, 89, 116], [18, 73, 116, 117], [18, 77, 116, 117], [51, 77, 116, 117], [18, 69, 77, 117], [13, 51, 86, 117], [13, 51, 77, 117], [13, 69, 77, 117], [18, 68, 69, 117], [18, 68, 73, 117], [13, 69, 91, 117], [13, 86, 91, 117], [68, 69, 91, 117], [5, 86, 91, 117], [5, 68, 73, 117], [5, 68, 91, 117], [50, 51, 116, 117], [5, 50, 86, 117], [50, 51, 86, 117], [5, 50, 73, 117], [50, 73, 116, 117], [11, 55, 83, 118], [3, 66, 75, 118], [3, 54, 75, 118], [3, 54, 83, 118], [54, 55, 83, 118], [11, 55, 79, 118], [11, 67, 79, 118], [66, 67, 93, 118], [3, 83, 93, 118], [3, 66, 93, 118], [11, 67, 93, 118], [11, 83, 93, 118], [22, 66, 67, 118], [22, 67, 79, 118], [22, 66, 75, 118], [54, 55, 87, 119], [54, 55, 118, 119], [55, 79, 118, 119], [54, 75, 118, 119], [22, 75, 118, 119], [22, 79, 118, 119], [22, 71, 79, 119], [22, 70, 71, 119], [70, 71, 95, 119], [22, 70, 75, 119], [7, 54, 75, 119], [7, 54, 87, 119], [7, 70, 75, 119], [7, 87, 95, 119], [7, 70, 95, 119], [15, 71, 79, 119], [15, 55, 79, 119], [15, 55, 87, 119], [15, 71, 95, 119], [15, 87, 95, 119]]
from itertools import combinations
F = all_intersections(A) # all intersections: function from other question
                         # takes 415 ms
F = sorted(F,lambda x,y: cmp(len(x),len(y)))
pairs = [ (x,y) for x,y in combinations(F,2) if set(y).issuperset(x) ]
                         # takes ~6 sec

An example is the square with vertices labelled with {1,2,3,4}: the set A is then {{1,2},{2,3},{3,4},{4,1}}, the intersections F are {{},{1},{2},{3},{4},{1,2},{2,3},{3,4},[4,1},{1,2,3,4}}, and the pairs in question are

({},{1}),({},{2}),({},{3}),({},{4}),
({1},{1,2}),({1},{4,1}),
({2},{1,2}),({2},{2,3}),
({3},{2,3}),({3},{3,4}),
({4},{3,4}),({4},{4,1}),
({1,2},{1,2,3,4}),({2,3},{1,2,3,4}),({3,4},{1,2,3,4}),({4,1},{1,2,3,4})

Once you are given the set F, I don't think there is anything better than just comparing the elements. But I am more thinking of an algorithm that computes (1) and (2) at the same time using the knowledge about stuff that was just intersected.

Following the solution by David K below, given the why, there are two more assumptions that can be used:

  1. The resulting order is graded with a unique minimal and a unique maximal element. This is, every maximal chain F0 < F1 < ... < Fm of cover relations has the same length and F0 is the empty set and Fm is the union of the input sets A. We call the set Fi to be of rank i, which is well-defined given the graded-ness.

  2. Every rank M set is the intersection of exactly 2 rank M+1 sets.

Many thanks!

Upvotes: 3

Views: 365

Answers (3)

Christian
Christian

Reputation: 577

Here is an alternative algorithm which only uses the graded-ness assumption above. The idea is that

  1. every intersection except the union of A has a cover which is computable rather quickly

  2. this produces a spanning tree of the order starting at the top element give by union of A

  3. one then computes the ranks successively and compare elements only for ranks that differ by 1

    from collections import defaultdict
    
    # this is John Coleman's function from
    # http://stackoverflow.com/questions/37622153
    # compute all intersections including the empty intersection
    # corresponding to the union of all sets in `A`.
    def all_intersections(A):
        # using frozensets for intersections
        A = map(frozenset,A)
        A.sort(key=len)
    
        # the union of A as the ground set
        universalSet = frozenset.union(*A)
    
        # computing the intersections successively
        intersections = set([universalSet])
        for s in A:
            moreIntersections = set(s & t for t in intersections)
            intersections.update(moreIntersections)
        return intersections
    
    def ranked_pieces(intersections):
        # this is to shortcut the length tests below
        lens = { s:len(s) for s in intersections }
        V = sorted(intersections, lambda x,y:cmp(lens[x],lens[y]))
    
        m = len(V)
    
        lower_covs = defaultdict(set)
        # we first compute the spanning tree...
        for i,x in enumerate(V):
            for j in range(i+1,m):
                y = V[j]
                if lens[x] < lens[y] and y.issuperset(x):
                    lower_covs[y].add(x)
                    # since V is sorted according to the size,
                    # we have surely found a cover and stop
                    break
    
        # ... and then the level sets
        level = set([V[-1]])
        level_sets = [level]
    
        while level:
            level = set.union(*[lower_covs[v] for v in level])
            level_sets.append(level)
    
        # the level sets are ordered backwards now
        # i.e., level_sets[0] is the biggest set
        # and level_sets[-1] is the empty set
        return level_sets
    
    def ranked_pieces_to_covers(ranked_pieces):
        # this is because we want tuples and not sets
        # and we want to compute them only once
        back_ref = { f:tuple(sorted(f)) for i in range(len(ranked_pieces)) for f in ranked_pieces[i] }
        covs = []
        for i in range(len(ranked_pieces)-1):
            high = ranked_pieces[i]
            low  = ranked_pieces[i+1]
    
            for x in low:
                for y in high:
                    if y.issuperset(x):
                        covs.append((back_ref[x], back_ref[y]))
        return covs
    
    def generate_covers(A):
        return ranked_pieces_to_covers(ranked_pieces(all_intersections(A)))
    

For the test data of the 600-cell, the algorithm by David K is much quicker:

    sage: %time X = generatePairs(A)
    CPU times: user 363 ms, sys: 28.2 ms, total: 392 ms
    Wall time: 373 ms
    sage: %time Y = generate_covers(A)
    CPU times: user 1.09 s, sys: 25.3 ms, total: 1.12 s
    Wall time: 1.08 s
    sage: set(X) == set(Y)
    True

but the "denser" the sets are, the better the other algorithm becomes:

The simplex on 10 vertices:

    sage: m = 10; A = [ [i for i in range(m) if i != j] for j in range(m) ]
    sage: %time X = generatePairs(A)
    CPU times: user 151 ms, sys: 4.33 ms, total: 156 ms
    Wall time: 144 ms
    sage: %time Y = generate_covers(A)
    CPU times: user 118 ms, sys: 17.3 ms, total: 136 ms
    Wall time: 106 ms
    sage: set(X) == set(Y)
    True

and the simplex on 14 vertices:

    sage: m = 14; A = [ [i for i in range(m) if i != j] for j in range(m) ]
    sage: %time X = generatePairs(A)
    CPU times: user 56.3 s, sys: 136 ms, total: 56.5 s
    Wall time: 56.6 s
    sage: %time Y = generate_covers(A)
    CPU times: user 31 s, sys: 65.4 ms, total: 31 s
    Wall time: 31 s

I mainly post the algorithm here for documentation, but I would certainly highly appreciate suggestions for improvements.

Upvotes: 0

David K
David K

Reputation: 3132

Here's a function that takes advantage of the assumption that the lists in the input are the facets of an abstract polytope. Instead of taking all intersections of collections of facets, this function assumes the input is a complete list of M-faces (polytopes of rank M) within a polytope of rank M + 1. It then performs a loop in which each iteration takes a complete list of M-faces and produces a complete list of (M-1)-faces, while accumulating all the containment pairs for those two lists of faces.

The main loop of the function intersects each pair of M-faces and builds a structure listing each intersection and the M-faces that contained it. These intersections include all the (M-1)-faces, but also include some lower-rank faces. The lower-rank faces can be identified by observing that each of them is a subset of an (M-1)-face, so any intersection which is a subset of another intersection is eliminated.

A rough breakdown of the running time is 40% to intersect pairs of faces, 40% to keep track of which M-faces contained each of the resulting intersections, 10% to eliminate faces of rank less than M - 1, and 10% to write the containment pairs to the output list. My computer seems to be slower than yours (about 8 seconds instead of 6 or 6.5 seconds for the original function), but the end result of the new function is a list of all the containment pairs between each rank and the next rank, about 10x-15x faster than the original function that produces all containment pairs (including the ones that "skip" ranks).

Note that not every list of lists of integers is valid input for the new function, because there are collections of sets of points that are not facets of an abstract polytope. I did not include code to check the input for correctness.

To check correctness of the output, I added some (rather slow) code to the original function to find all pairs (s,t) in the (original) output list such that pairs of the form (s,u) and (u,t) were also in the list, and then return a revised list with all those pairs removed. I also modified both the new and old functions by invoking sorted() on each list of integers they put in the output so that the output lists would compare correctly. I then confirmed that both functions produced identical output.

By the way, I doubt this function is a pythonic as it could have been. Comments suggesting improvements in that matter are welcome.

from collections import defaultdict
import sys

def generatePairs(A):
    # It is assumed that A consists exactly of all the facets of an abstract
    # polytope of rank N; that is, the abstract polytope is a graded poset
    # in which the minimal element is the empty set and has rank -1, the
    # maximal element is the polytope's body, which has rank N, and A
    # contains all facets of the polytope, which have rank N - 1.
    # Then within the graded poset,
    # each element of rank 0 is a point and has cardinality 1;
    # each element of rank 1 is an edge and has cardiality 2;
    # each element of rank M (where M > 1) is a rank-M polytope and has
    # cardinality at least M + 1, but may have greater cardinality.

    # We start with the facets (rank N-1).
    rank_to_intersect = [frozenset(s) for s in A]

    # Construct the body (rank N).
    polytope_body = list(frozenset.union(*rank_to_intersect))
    body_size = len(polytope_body)

    # covering_pairs will be all the pairs of polytopes (s,t) such that
    # rank(s) + 1 == rank(t) and s is a subset of t. Initially we populate
    # it with just the pairs whose ranks are respectively N-1 and N.
    covering_pairs = [(s, polytope_body) for s in A]

    while (len(rank_to_intersect) > 0) and (len(rank_to_intersect[0]) > 2):
        # For some integer M such that M > 1, rank_to_intersect contains all
        # the polytopes of rank M. At the end of each iteration of the loop,
        # rank_to_intersect will contain all the polytopes of rank M - 1.
        # Also, all the pairs (x,y) where rank(x) = M - 1 and rank(y) = M
        # will have been added to covering_pairs.

        container_map = defaultdict(list)
        while rank_to_intersect:
            s = rank_to_intersect.pop()
            for t in rank_to_intersect:
                x = s & t
                if len(x) > 1:
                    container_map[x].extend([s, t])
                    # Note that the list container_map[x]
                    # may contain duplicates

        # The keys of container_map, consisting of all pairwise
        # intersections of polytopes of rank M, include all polytopes
        # of rank M - 1 but also some polytopes of lower ranks.
        # Any polytope of a lower rank, however, is a subset of
        # a polytope of rank M - 1 that is also in the list.

        min_size   = min([len(s) for s in container_map.keys()])
        max_size   = max([len(s) for s in container_map.keys()])
        size_range = range(min_size, max_size + 1)
        candidates = dict([(i, []) for i in size_range])
        for s in container_map.keys():
            candidates[len(s)].append(s)

        # Repopulate rank_to_intersect with the polytopes of rank M - 1.
        for set_size in size_range:
            larger_sizes = range(set_size + 1, max_size + 1)
            for s in candidates[set_size]:
                if not any(any(t >= s for t in candidates[i])
                           for i in larger_sizes):
                    # We now know that s has rank M - 1, not a lower rank.
                    rank_to_intersect.append(s)

        # Add all the (rank-(M - 1), rank-M) pairs to covering_pairs.
        for s in rank_to_intersect:
            # container_map[s] may contain duplicates; avoid them.
            containers = frozenset(container_map[s])
            covering_pairs.extend([(list(s), list(t)) for t in containers])

    # At the end of the loop, rank_to_intersect contains the rank-1
    # polytopes, that is, the edges.
    # Each edge contains each of its two endpoints.
    points_with_duplicates = []
    for e in rank_to_intersect:
        covering_pairs.extend([([p], list(e)) for p in e])
        points_with_duplicates.extend(e)

    # List the containment pairs of the empty set without duplicating points.
    points = frozenset(points_with_duplicates)
    covering_pairs.extend([([], [p]) for p in points])

    return covering_pairs

Upvotes: 2

beoliver
beoliver

Reputation: 5759

could you just let me know the running time for the following changes?

def all_intersections(lists):
    sets = allIntersections([frozenset(s) for s in lists])
    return list(sets)


A = ...
F = all_intersections(A) 
F.sort(key=len)
pairs = [(x,y) for x,y in combinations(F,2) if y.issuperset(x)]

Upvotes: 0

Related Questions