## Last edited 11/22/2018 by Joseph Doolittle and Bennet Goeckner.
##
## The following is our counterexample to Conjecture 2.4 from Stanley's "A combinatorial decomposition of acyclic simplicial complexes."
## We will first show that the complex is 2-fold acyclic (i.e., link(F) is acyclic for all F such that |F|<2).
## Then we will show, by use of a linear program, that its face poset cannot be fully decomposed into disjoint rank 2 boolean intervals.
##
## Below we display our example's f-vector. Our example has 20 vertices and 99 facets, and it is pure and three-dimensional.
## More details appear in our preprint.
##
Example=SimplicialComplex({ 
        #Gamma
        ('A','B','C','E'),('B','C','E','F'),('B','C','D','F'),('A','B','C','G'),('B','C','G','H'),('B','C','D','H'),('A','B','E','G'),('B','E','F','G'),('B','F','H','G'),
        #Delta1
        ('A','B','C','J1'),('A','B','I1','J1'),('B','C','I1','J1'),('B','C','D','I1'),('C','D','I1','J1'),
        ('A','B','E','K1'),('A','B','I1','K1'),('B','E','I1','K1'),('B','E','F','I1'),('E','F','I1','K1'),
        ('A','B','G','L1'),('A','B','I1','L1'),('B','G','I1','L1'),('B','G','H','I1'),('G','H','I1','L1'),
        ('C','D','F','K1'),('C','D','J1','K1'),('C','F','J1','K1'),('C','E','F','J1'),('E','F','J1','K1'),
        ('C','D','H','L1'),('C','D','J1','L1'),('C','H','J1','L1'),('C','G','H','J1'),('G','H','J1','L1'),
        ('E','F','G','L1'),('E','F','K1','L1'),('F','G','K1','L1'),('F','G','H','K1'),('G','H','K1','L1'),
        #Delta2
        ('A','B','C','J2'),('A','B','I2','J2'),('B','C','I2','J2'),('B','C','D','I2'),('C','D','I2','J2'),
        ('A','B','E','K2'),('A','B','I2','K2'),('B','E','I2','K2'),('B','E','F','I2'),('E','F','I2','K2'),
        ('A','B','G','L2'),('A','B','I2','L2'),('B','G','I2','L2'),('B','G','H','I2'),('G','H','I2','L2'),
        ('C','D','F','K2'),('C','D','J2','K2'),('C','F','J2','K2'),('C','E','F','J2'),('E','F','J2','K2'),
        ('C','D','H','L2'),('C','D','J2','L2'),('C','H','J2','L2'),('C','G','H','J2'),('G','H','J2','L2'),
        ('E','F','G','L2'),('E','F','K2','L2'),('F','G','K2','L2'),('F','G','H','K2'),('G','H','K2','L2'),
        #Delta3
        ('A','B','C','J3'),('A','B','I3','J3'),('B','C','I3','J3'),('B','C','D','I3'),('C','D','I3','J3'),
        ('A','B','E','K3'),('A','B','I3','K3'),('B','E','I3','K3'),('B','E','F','I3'),('E','F','I3','K3'),
        ('A','B','G','L3'),('A','B','I3','L3'),('B','G','I3','L3'),('B','G','H','I3'),('G','H','I3','L3'),
        ('C','D','F','K3'),('C','D','J3','K3'),('C','F','J3','K3'),('C','E','F','J3'),('E','F','J3','K3'),
        ('C','D','H','L3'),('C','D','J3','L3'),('C','H','J3','L3'),('C','G','H','J3'),('G','H','J3','L3'),
        ('E','F','G','L3'),('E','F','K3','L3'),('F','G','K3','L3'),('F','G','H','K3'),('G','H','K3','L3'),
        })
print "The counterexample's f-vector is" , Example.f_vector()

#Example=SimplicialComplex({(0, 1, 8, 9), (2, 6, 7, 13), (5, 7, 10, 11), (1, 6, 7, 12), (4, 5, 13, 14), (2, 3, 8, 9), (2, 6, 9, 11), (1, 2, 3, 12), (0, 1, 7, 19), (0, 1, 4, 14), (0, 1, 7, 11), (2, 3, 16, 17), (4, 5, 7, 15), (0, 1, 12, 15), (0, 1, 4, 18), (0, 1, 2, 9), (6, 7, 13, 15), (0, 1, 16, 19), (1, 7, 12, 15), (1, 2, 8, 9), (1, 2, 3, 16), (1, 2, 6, 7), (5, 6, 7, 10), (1, 2, 3, 6), (2, 5, 9, 10), (2, 3, 9, 10), (2, 3, 6, 19), (1, 2, 16, 17), (1, 4, 5, 7), (4, 5, 18, 19), (2, 4, 5, 17), (0, 1, 2, 13), (6, 7, 18, 19), (0, 1, 2, 17), (2, 3, 12, 13), (1, 4, 8, 10), (0, 1, 2, 7), (6, 7, 16, 19), (2, 3, 13, 15), (1, 6, 7, 16), (2, 6, 7, 9), (0, 1, 4, 10), (4, 5, 7, 11), (2, 6, 13, 15), (0, 1, 7, 15), (1, 5, 6, 7), (1, 4, 5, 16), (2, 4, 5, 9), (4, 5, 8, 10), (1, 2, 3, 8), (1, 4, 5, 8), (1, 2, 12, 13), (0, 1, 2, 4), (2, 5, 13, 14), (2, 3, 13, 14), (0, 1, 8, 10), (1, 4, 5, 12), (2, 6, 7, 17), (2, 3, 6, 15), (2, 3, 5, 18), (5, 7, 18, 19), (0, 1, 8, 11), (1, 7, 8, 11), (5, 6, 7, 18), (2, 3, 9, 11), (1, 2, 3, 5), (2, 3, 17, 19), (6, 7, 17, 19), (2, 5, 17, 18), (2, 3, 17, 18), (6, 7, 9, 11), (1, 2, 4, 5), (4, 5, 16, 18), (4, 5, 10, 11), (0, 1, 12, 13), (1, 7, 16, 19), (2, 4, 5, 13), (6, 7, 10, 11), (4, 5, 9, 10), (2, 3, 5, 14), (1, 4, 12, 14), (4, 5, 7, 19), (4, 5, 14, 15), (4, 5, 17, 18), (1, 4, 16, 18), (1, 6, 7, 8), (6, 7, 14, 15), (0, 1, 16, 18), (5, 6, 7, 14), (6, 7, 12, 15), (5, 7, 14, 15), (6, 7, 8, 11), (0, 1, 4, 7), (0, 1, 16, 17), (2, 6, 17, 19), (2, 3, 5, 10), (0, 1, 12, 14), (2, 3, 6, 11), (4, 5, 12, 14)})
#print "The counterexample's f-vector is" , Example.f_vector()

## The complex itself and the links of all vertices are acyclic; thus the complex is 2-fold acyclic:
##
print "Homology of complex :", Example.homology()
for v in Example.vertices():
    print "Homology of vertex" , v , ":" , Example.link([v]).homology()

## The following linear program tries to decompose a given 3-dimensional complex's face poset into rank 2 boolean intervals (here called ``diamonds'').
##
def DiamondDecomp(Q): #Counts the maximum number of diamonds that can appear in a disjoint decomposition. Specific to 3 dimensional complexes.
    Diamonds = []
    for f in Q.faces()[3]: #Adds a diamond for every possible interval which is topped by a 3-face.
        for e in Subsets(f.set()):
            if len(e)==2:
                D=[]
                for vs in Subsets(f.set().difference(set(e))):
                    D.append(set(e).union(vs))
                Diamonds.append(list(D))
    for t in Q.faces()[2]: #Adds a diamond for every possible interval which is topped by a 2-face.
        for v in t.set():
            D=[]
            for vs in Subsets(t.set().difference(set([v]))):
                D.append(set([v]).union(vs))
            Diamonds.append(list(D))
    for e in Q.faces()[1]: #Adds a diamond for every possible interval which is topped by a 1-face.
        D=[]
        for vs in Subsets(e.set()):
            D.append(set(vs))
        Diamonds.append(list(D))

    p = MixedIntegerLinearProgram(maximization=True)
    v = p.new_variable(binary=True) #Makes it a binary program.
    for sigma in Q.face_iterator(): #For each face of the complex, constrains so that the number of diamonds containing that face is less than or equal to 1.
        p.add_constraint(sum(v[i] * (1 if set(sigma) in Diamonds[i] else 0) for i in range(len(Diamonds))) <= 1)
    p.set_objective(sum(v[i] for i in range(len(Diamonds)))) #Maximize the number of the diamonds.
    result=p.solve()
    return result

## The f-polynomial of our example is f(t) = (1+t)^2*(1+18t+99t^2), so we would need 1+18+99=118 rank 2 boolean intervals to fully decompose the face poset.
## However, we will see below that the largest decomposition the linear program can find contains only 117 rank 2 boolean intervals, therefore the complex is a counterexample to Stanley's conjecture.
##
print "The largest possible decomposition of our example into rank 2 boolean intervals contains only" , DiamondDecomp(Example) , "intervals. It would need 118 intervals to be fully decomposed."