%auto
##############################################################
# INITIALIZE THIS CELL
##############################################################

from sage.graphs.independent_sets import IndependentSets
from sage.combinat.cyclic_sieving_phenomenon import *

##############################################################
# code for maximal orthogonal pairs
##############################################################
def inc(a,b):
    if a[0].issubset(b[0]) and a[1].issubset(b[1]):
        return True
    else:
        return False
def inc1(a,b):
    if a[0].issubset(b[0]):
        return True
    else:
        return False
def mops(g):
    sets=[]
    edges=[(i[0],i[1]) for i in g.edges()]
    for s in Subsets(g):
        for t in Subsets(Set(g).difference(s)):
            no=False
            for i in s:
                for j in t:
                    if (i,j) in edges:
                        no=True
                        break
                if no==True:
                    break
            if no==False:
                sets.append((s,t))
    P=Poset([sets,inc])
    Pmax=P.maximal_elements()
    return Pmax
def poset_of_mops(g):
#naively, and slowly, computes the poset of maximal orthogonal pairs from the directed graph g
    return Poset([mops(g),inc1])
def complete_mop(g,s):
#naively, and slowly, computes the ways to complete s to a maximal orthogonal pair
    sets=[]
    edges=[(i[0],i[1]) for i in g.edges()]
    for t in Subsets(Set(g).difference(s)):
        no=False
        for i in s:
            for j in t:
                if (i,j) in edges:
                    no=True
                    break
            if no==True:
                break
        if no==False:
            sets.append((s,t))
    P=Poset([sets,inc])
    return P.maximal_elements()
##############################################################
# code for left modular lattices
##############################################################
def is_left_modular(L,H=False,verbose=False):
#checks if the subset of elements H of the lattice L is left modular; if verbose=True, then gives a list of failures
    if H==False:
        H=L
    for x in H:
        for z in L:
            for y in L.principal_lower_set(z):
                if (L.join(y, L.meet(x, z)) != L.meet(L.join(y, x), z)):
                    if verbose==False:
                        return False
                    else:
                        out+=[(y,x,z)]
    if verbose==False:
        return True
    else:
        return out
def join_irr_label(L,c,j):
#given a left modular lattice L, a left modular chain c, and a join-irreducible element, returns the corresponding label
    for i in range(len(c)):
        if L.le(j,c[i]):
            return i
def meet_irr_label(L,c,m):
#given a left modular lattice L, a left modular chain c, and a meet-irreducible element, returns the corresponding label
    for i in range(len(c)-1,-1,-1):
        #print (m,c[i])
        if L.ge(m,c[i]):
            return i+1
def downward_labels(L,c,z):
#given a trim lattice L, a left modular chain c, and an element z, returns the downward labels for z
    J=L.join_irreducibles()
    JL=dict([(j,join_irr_label(L,c,j)) for j in J])
    cov=L.lower_covers(z)
    return [min([JL[j] for j in J if L.join([y,j])==z]) for y in cov]
def upward_labels(L,c,z):
#given a trim lattice L, a left modular chain c, and an element z, returns the upward labels for z
    M=L.meet_irreducibles()
    ML=dict([(m,meet_irr_label(L,c,m)) for m in M])
    cov=L.upper_covers(z)
    return [max([ML[m] for m in M if L.meet([y,m])==z]) for y in cov]
def cover_label(L,c,edge):
#given a trim lattice L, a left modular chain c, and a covering relation, returns the label for that relation
    return min([JL[j] for j in J if L.join([edge[0],j])==edge[1]])
##############################################################
# trim lattices
##############################################################
def is_trim(L):
#checks if the lattice L is trim; computes all maximal length chains of L, which can take a while.  If it is true, it returns a (True,list of maximal chains)
    J=L.join_irreducibles()
    M=L.meet_irreducibles()
    C=L.maximal_chains()
    n=max(map(len,C))
    #print n,len(J),len(M)
    mC=[c for c in C if len(c)==n]
    #print is_left_modular(mC[0])
    if len(J)==len(M) and len(J)==n-1 and is_left_modular(L,mC[0]):
        return (True,mC)
    else:
        return (False,False)
#def trim_rowmotion(L,c):
#    Ple=P.linear_extension()
#    def cyc_act(d):
#        return tuple(togud_seq(d,Ple))
#    map(len,orbit_decomposition(L,cyc_act))
#def is_sub(a,b):
#    return a.issubset(b)
#def noncrossing_lattice(L):
#returns the noncrossing partition lattice for the lattice L represented as mops; in general, we only expect this to have nice properties when L is a twisted Cambrian lattice
#    N=[]
#    for i in L:
#        N+=[i[0].intersection(trim_rowmotion(L,i)[1])]
#    return Poset([N,is_sub])
##############################################################
# tight repellent independent pairs
##############################################################
def complete_trip(G,u):
#given an independent set u in the digraph G, return the trip of the form (J,u)
    m=Set([])
    H=DiGraph(G)
    H.reverse_edges(H.edges())
    P=Poset(G)
    Q=Poset(H)
    m=Set([])
    for k in P.linear_extension():
        if Set([i[1] for i in H.edges_incident(k)]).intersection(m)==Set([]) and Set([i[1] for i in G.edges_incident(k)]).intersection(u)==Set([]) and k not in u:
            m=m.union(Set([k]))
    return (m,u)
def trip_to_mop(G,du):
#computes the mop corresponding to a trip
    pass
def maximal_trip(G):
#given a directed graph G, return the maximal trip of the form (J,[])
    return complete_trip(G,Set([]))
def minimal_trip(G):
#given a directed graph G, return the minimal trip of the form ([],M)
    H=DiGraph(G)
    H.reverse_edges(H.edges())
    P=Poset(G)
    Q=Poset(H)
    m=Set([])
    for k in Q.linear_extension():
        if Set([i[1] for i in G.edges_incident(k)]).intersection(m)==Set([]):
            m=m.union(Set([k]))
    return (Set([]),m)
def flip_up(G,du,j):
#given a directed graph G, a trip (D,U) and an element j in u, return flip_j(D,U)
    a=du[0]
    b=du[1]
    if j not in b:
        return du
    elif j in b:
        return flip(G,du,j)
def flip(G,du,j):
#given a directed graph G, a trip (D,U) and an element j of G, return flip_j(D,U)
    H=DiGraph(G)
    H.reverse_edges(H.edges())
    P=Poset(G)
    Q=Poset(H)
    a=du[0]
    b=du[1]
    if j not in b and j not in a:
        return du
    elif j in b:
        a2=Set([k for k in a if not(Q.le(k,j))]+[j])
        b2=Set([k for k in b if not(Q.ge(k,j))])
    elif j in a:
        a2=Set([k for k in a if not(Q.le(k,j))])
        b2=Set([k for k in b if not(Q.ge(k,j))]+[j])
    #print "first step: ",(a2,b2)
    for k in P.linear_extension():
        if not(Q.ge(k,j)) and Set([i[1] for i in G.edges_incident(k)]).intersection(b)==Set([]) and Set([i[1] for i in H.edges_incident(k)]).intersection(a2)==Set([]) and (k not in b2):
            a2=a2.union(Set([k]))
    for k in Q.linear_extension():
        if not(Q.le(k,j)) and Set([i[1] for i in H.edges_incident(k)]).intersection(a)==Set([]) and Set([i[1] for i in G.edges_incident(k)]).intersection(b2)==Set([]) and (k not in a2):
            b2=b2.union(Set([k]))
    return (a2,b2)
def flips(G,du,s):
#given a directed graph G, a trip (D,U) and a sequence s of elements j of G, return (\prod_{j \in s} flip_j)(D,U)
    dv=tuple(du)
    for i in s:
        dv=flip(G,dv,i)
    return dv
def make_flip_poset(G):
#given a directed graph G, return the corresponding independence poset
    edges=[]
    els=Set([])
    last=[]
    new=[minimal_trip(G)]
    while new!=[]:
        last=new
        els=els.union(Set(last))
        new=[]
        for a in last:
            for i in a[1]:
                b=flip(G,a,i)
                edges.append([a,b])
                if b not in new:
                    new=new+[b]
    return Poset((list(els),edges))
def make_flip_tree(G):
#given a directed graph G, return the corresponding flip tree
    edges=[]
    els=Set([])
    last=[]
    new=[(min_ind(),-1)]
    while new!=[]:
        last=new
        els=els.union(Set(last))
        new=[]
        for a in last:
            for i in a[0][1]:
                j=Q.linear_extension().index(i)
                if a[1]==-1 or j>a[1]:
                    b=tog(a[0],i)
                    edges.append([a,(b,j)])
                    if b not in new:
                        new=new+[(b,j)]
    return Poset((list(els),edges))
def rowmotion(L,x):
#compute rowmotion on the trip x
    lower=dict([(l,l[0]) for l in L])
    upper=dict([(l[1],l) for l in L])
    return upper[lower[x]]

def tog(G,x):
    H=G
    for i in H.edges():
        if x in i:
            H.reverse_edge(i)
    return(H)
    #ed=G.edges_incident(x)
    #return G.reverse_edges(ed,inplace=False)

##############################################################
# Galois graphs
##############################################################
def semidistributive_Galois_graph(L):
#returns the Galois graph for a semidistributive lattice L
    lower=dict([(l,Set(L.canonical_joinands(l))) for l in L])
    J=L.join_irreducibles()
    M=L.meet_irreducibles()
    meet_to_join=dict([(m,L.meet([l for l in J if L.le(l,L.upper_covers(m)[0]) and not(L.le(l,m))])) for m in M])
    upper=dict([(Set([meet_to_join[i] for i in L.canonical_meetands(l)]),l) for l in L])
    return DiGraph([ [j for j in J] ,[(j,meet_to_join[m]) for j in J for m in M if j!=meet_to_join[m] and not(L.le(j,m))]])
def trim_Galois_graph(L,c):
#returns the Galois graph for a trim lattice L, using the chain c
        J=L.join_irreducibles()
        M=L.meet_irreducibles()
        JL=dict([(j,join_irr_label(L,c,j)) for j in J])
        lower=dict([(l,Set(downward_labels(L,c,l))) for l in L])
        return DiGraph([ [join_irr_label(L,c,j) for j in J] ,[(join_irr_label(L,c,j),meet_irr_label(L,c,m)) for j in J for m in M if join_irr_label(L,c,j)!=meet_irr_label(L,c,m) and not(L.le(j,m))]])
def Weyl_word_Galois_graph(W,w):
#returns the Galois graph for the word w in the Weyl group W; when w is the c-sorting word for w_o, this is the Galois graph for the Cambrian lattice
    edges=[]
    for i in range(len(w)):
        r=[t.reflection_to_root() for t in W.inversion_sequence(w[i:])]
        for j in range(1,len(r)):
            if (r[j]-r[0]).is_positive_root():# or (r[j]-r[0])==0:
                edges.append([j+i,i])
    return DiGraph(edges)
def Cambrian_Galois_graph(W,c):
    return Weyl_word_Galois_graph(W,W.w0.coxeter_sorting_word(c))
def twisted_Galois_graph(W,c,v):
    w=v.inverse()
#Galois graph for twisted Cambrian; here w is c-sortable with w \geq c in weak order (i.e. ``positive'')
    G0=Weyl_word_Galois_graph(W,W.w0.coxeter_sorting_word(c))
    wr=W.inversion_sequence(W.w0.coxeter_sorting_word(c))
    wc=[wr.index(i) if i in wr else -1 for i in W.inversion_sequence(list(c)+W.w0.coxeter_sorting_word(c))]
    #print wc
    for i,j in enumerate(wc):
          if wc.index(j)!=i:
                wc[i]=False
    wc=wc[len(c):]
    #print wc
    s=Set([wr.index(i) for i in w.inversions()])
    #print s
    x=complete_mop(G0,s)[0]
    #print x
    ed=[(i[0],i[1]) for i in G0.edges()]
    #if Set(range(len(c))).issubset(x[0]):
    G1=G0.subgraph(x[0])
    G2=G0.subgraph(x[1])
    G=G1.union(G2)
    for g2 in G2:
        for g1 in G1:
            if wc[g1]!=False and ((g2,wc[g1]) in ed or g2==wc[g1]):
                G.add_edge(g1,g2)
    return G
    #else:
    #    return False

##############################################################
# Cambrian lattices
##############################################################
def cambrian_intervals(W,c):
#returns all Cambrian intervals in the Cambrian lattice for W,c
    WW=W.weak_lattice(facade=True)
    sort=[w for w in W if w.is_coxeter_sortable(c)]
    intervals=dict([])
    for w in W:
        t=WW.subposet([s for s in sort if s.weak_le(w)]).maximal_elements()[0]
        if t in intervals:
            intervals[t]=intervals[t]+(w,)
        else:
            intervals[t]=(w,)
    return intervals
def Cambrian_lattice(W,c):
    return LatticePoset(make_flip_poset(Cambrian_Galois_graph(W,c)))
def twisted_Cambrian_lattice(W,c,w):
#returns the twisted Cambrian; here w is c-sortable with w \geq c in weak order (i.e. ``positive'')
    return make_flip_poset(twisted_Galois_graph(W,c,w))
def uniq(input):
    output = []
    for x in input:
        if x not in output:
            output.append(x)
    return output

# Builds an example of the independence poset for the directed graph on [a]x[b]
# This can take a long time if you increase a or b
a=2
b=3
G=DiGraph([[(i,j) for i in range(a) for j in range(b)],[[(i,j),(i+1,j)] for i in range(a-1) for j in range(b)]+[[(i,j),(i,j+1)] for i in range(a) for j in range(b-1)]])
#view(G)
IndPoset=make_flip_poset(G)
IndPoset.show() # show the independence poset (note that arrows point towards *larger* elements with Sage conventions)

# Shows an example of a flip on the maximal element of the independence poset IndPoset
# (which must contain the soures of the digraph G)
p_max=IndPoset.maximal_elements()[0]
print "p_max =",p_max
print "flip_g(p_max) at a source g =",flip(G,p_max,G.sources()[0]),"\n"

# Shows an example of a flip on the minimal element of the independence poset IndPoset
# (which must contain the sinks of the digraph G)
p_min=IndPoset.minimal_elements()[0]
print "p_min =",p_min
print "flip_g(p_min) at a sink g =",flip(G,p_min,G.sinks()[0])

# Examples: Building posets of maximal orthogonal pairs (mops) and tight orthogonal pairs (tops)
# These posets are isomorphic iff the top poset is a trim lattice iff the top poset is a lattice

G=DiGraph([[0,1],[1,2]])
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "Posets of mops and tops are isomorphic?",M.is_isomorphic(F)
print "Poset of tops is a lattice?",F.is_lattice()
print "The maximal length of a chain in the poset of mops is",max(map(len,M.maximal_chains()))
print "The maximal length of at chain in the poset of tops is",max(map(len,F.maximal_chains()))
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))
print "\n"

G=DiGraph([[0,1],[1,2],[2,3]])
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "Posets of mops and tops are isomorphic?",M.is_isomorphic(F)
print "Poset of tops is a lattice?",F.is_lattice()
print "The maximal length of a chain in the poset of mops is",max(map(len,M.maximal_chains()))
print "The maximal length of at chain in the poset of tops is",max(map(len,F.maximal_chains()))
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))
print "\n"

G=DiGraph([[0,1],[1,2],[2,3],[3,4]])
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "Posets of mops and tops are isomorphic?",M.is_isomorphic(F)
print "Poset of tops is a lattice?",F.is_lattice()
print "The maximal length of a chain in the poset of mops is",max(map(len,M.maximal_chains()))
print "The maximal length of at chain in the poset of tops is",max(map(len,F.maximal_chains()))
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))
print "\n"

G=DiGraph([[0,1],[1,2],[2,3],[3,4],[4,5]])
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "Posets of mops and tops are isomorphic?",M.is_isomorphic(F)
print "Poset of tops is a lattice?",F.is_lattice()
print "The maximal length of a chain in the poset of mops is",max(map(len,M.maximal_chains()))
print "The maximal length of at chain in the poset of tops is",max(map(len,F.maximal_chains()))
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))
print "\n"

# Example: distributive lattices
P=Poset([[0,1],[[0,1],[0,2]]])
L=P.order_ideals_lattice()
trim=is_trim(L)
print "The distributive lattice L is trim: ",trim[0]
G=trim_Galois_graph(L,trim[1][0]) # compute a Galois graph from a trim lattice
G2=DiGraph([[1,0],[2,0]]) # constructing a new Galois graph G2
print "The Galois graph of L is isomorphic to a certain directed path: ",G.is_isomorphic(G2)
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "L is isomorphic to the lattice of maximal orthogonal sets: ",M.is_isomorphic(L)
print "L is isomorphic to the poset of tight orthogonal sets: ",F.is_isomorphic(L)
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))

# Example: Cambrian lattices
W=WeylGroup(['A',2])
c=W.from_reduced_word([1,2])
L=W.m_cambrian_lattice(c,1).relabel()
trim=is_trim(L)
print "The Cambrian lattice L of type A_2 is trim: ",trim[0]
G=trim_Galois_graph(L,trim[1][0]) # compute a Galois graph from a trim lattice
G2=DiGraph([[0,1],[1,2]]) # constructing a new Galois graph G2
print "The Galois graph of L is a directed path: ",G.is_isomorphic(G2)
M=poset_of_mops(G) # make the lattice of maximal orthogonal pairs corresponding to G
F=make_flip_poset(G) # make the poset of tight orthogonal pairs coresponding to G
print "L is isomorphic to the lattice of maximal orthogonal sets: ",M.is_isomorphic(L)
print "L is isomorphic to the poset of tight orthogonal sets: ",F.is_isomorphic(L)
print "The rowmotion orbits on the poset of tops have lengths",map(len,orbit_decomposition(F,lambda x: rowmotion(F,x)))