A certain polyhedron-related computation: any way to get it quicker?

edit

(This is a white flag sort of question.) I have a numeric conjecture that I was trying to verify. The conjecture is this:

• take a polyhedron P whose vertices are a poset
• extend this to a (non-reflexive) operation on all faces by saying that F is smaller than G if max vertex of F is
  smaller than min vertex of G
• take an algebra R of rational functions on generators x_F for all faces of P, plus extra generator t
• define the following endomorphism f of R:
  • for a face F, let Ch(F) be the set of all chains in F, where for a chain C = F1 <... < Fn its total codimension is codim(F,C) = dimF - sum dim Fi
  • for a chain C = F1 < ... < Fn, define its monomial m(C) by taking the product of either x_Fi, if dim F_i>0, or of x_Fi/(1-x_Fi), if Fi is a vertex, and then times t^{codim(F,C)}
  • define f by sending x_F to the sum of m(C) for all chains C in Ch(C)
• additionally define an endomorhism I by sending x_F to -x_F and t to -t

CONJECTURE: very often f*I will be an involution.

The conjecture is proved for simplices, cubes, and all polygons. It's also disproved in some cases. A key example where I need to know the answer is the associahedron, with Tamari order on vertices.

I've asked several questions here and, thanks to tmonteil, wrote a code that checks the conjecture for the 3D associahedron -- and... my code can't run in any reasonable amount of time! Does the problem above sound like something computationally approachable, given that the number of chains in the 3-dimensional face of associahedron is equal to 22663? How would you write this code?

Attached, my (very bad) code:

---

import itertools

A = Polyhedron(vertices = [[1,2,6,1],[1,6,2,1],[2,1,6,1],[1,2,3,4],[1,4,1,4],[1,6,1,2],[2,1,3,4],[4,1,4,1],[4,3,2,1],[4,3,1,2],[3,1,2,4],[3,2,1,4],[4,1,2,3],[4,2,1,3]])

Fempty = list(A.face_generator())
F = []
for i in Fempty:
    if i.dim() != -1:
        F.append(i)

R = Frac(PolynomialRing(QQ,x,len(F)+1))
R.inject_variables()

face = {x:f for x,f in list(zip(R.gens(),F))}  
x = {f:x for x,f in list(zip(R.gens(),F))}

V = VectorSpace(QQ,4)
W = Cone([[1,-1,0,0],[1,0,-1,0],[1,0,0,-1],[0,1,-1,0],[0,1,0,-1],[0,0,1,-1]])


elms = A.vertices()
rels = []
for i in A.faces(1):
    Li = i.vertices()
    if vector(list(Li[0]))-vector(list(Li[1])) in W:
        rels.append([Li[0],Li[1]])
    else:
        rels.append([Li[1],Li[0]])

P = Poset((elms, rels), cover_relations = True, facade = True)

def flarger(f1,f2):
    L1 = f1.vertices()
    L2 = f2.vertices()
    for v in L1:
        for w in L2:
            if (P.lt(w,v) or v == w) == False:
                return False
    return True

bigelms = A.faces(0)+A.faces(1)+A.faces(2)+A.faces(3)

bigrels = []

for i in bigelms:
    for j in bigelms:
        if (flarger(i,j) and i != j):
            bigrels.append([i,j])

Q = Poset((bigelms, bigrels), cover_relations = False, facade = True)

elll = list(x[i] for i in bigelms)
relll = []
for i in bigrels:
    relll.append([x[i[0]],x[i[1]]])

QQ = Poset((elll, relll), cover_relations = False, facade = True)

def inside(f1,f2):
    L1 = face[f1].vertices()
    L2 = face[f2].vertices()
    for v in L1:
        if v not in L2:
                return False
    return True

def subfaces(feis):
    result = []
    for item in elll:
        if inside(item,feis):
            result.append(item)
    return result

def sfposet(feis):
    Res = Poset(QQ.subposet(subfaces(feis)),facade = True)
    return Res


def chains(myx):
    Huj = sfposet(myx).chains()
    result = []
    for i in Huj:
        if i != []:
            result.append(i)
    return result


def p(feis):
    answer = feis
    if face[feis].dim() == 0:
        answer = feis/(1-feis)
    return answer

def monprod(monlist):
    answer = 1
    for item in monlist:
        answer = answer*p(item)
    return answer

def chdim(ch):
    dim = 0
    for item in ch:
        dim = dim+face[item].dim()
    return dim

def fmap(feis):
    answer = 0
    dimf = face[feis].dim()
    Flist = chains(feis)
    for item in Flist:
        power = dimf - chdim(item)
        answer = answer + monprod(item)*(x45^power)
    return answer


minuslist = []

for feis in F:
    minuslist.append(-x[feis])

minuslist.append(-x45)

fmaplist = []

for feis in F:
    fmaplist.append(fmap(x[feis]))

fmaplist.append(x45)

f = R.hom(fmaplist)
I = R.hom(minuslist)

f*I*f*I