Contact Us!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

Sage worksheet for paper "Hook formulas for skew shapes" by Alejandro Morales, Igor Pak and Greta Panova.

Project: excited_skew
Views: 35
Tweet
Share
Share
# Supplementary Sage worksheet for the article "Hook formulas for skew shapes", Alejandro H. Morales, Igor Pak, Greta Panova
# Version 1.0 2015-12-28
# Tested on Version 6.9.beta7 of SageMathCloud
version()
'SageMath Version 6.9.beta7, Release Date: 2015-09-17' 
# Code Block 1: Excited diagrams
class Board(object):

    def __init__(self, r, c, lam, entries, dist_entries = Set([]), other_dist_entries = Set([])):
        self.r = r
        self.c = c
        self.lam = lam
        self.lamc = Partition(lam).conjugate()
        self.entries = entries
        if  not dist_entries.issubset(self.entries):
            raise ValueError("%s must be a subset of %s" % (dist_entries, entries))
        else:
            self.dist_entries = dist_entries
        if  not other_dist_entries.issubset(self.entries):
            raise ValueError("%s must be a subset of %s" % (other_dist_entries, entries))
        else:
            self.other_dist_entries = other_dist_entries


    def __repr__(self):
        return str(self.print2Dtext())

    def __str__(self):
        return str(self.print2Dtext())

    def __eq__(self, other):
        if self.r == other.r and self.c == other.c and self.entries == other.entries and self.dist_entries == other.dist_entries:
            return True
        else:
            return False

    def __ne__(self, other):
        return (not self.__eq__(other))

    def __hash__(self):
        return hash(self.__repr__())

    def getRows(self):
        return self.r

    def getCols(self):
        return self.c

    def getEntries(self):
        return self.entries

    def getNonDistEntries(self):
        return self.entries.symmetric_difference(self.dist_entries)

    def getDistEntries(self):
        return self.dist_entries

    def getComp(self):
        allent = getAllBoard(self.r, self.c)
        return Board(self.r,self.c, allent.symmetric_difference(self.entries))

    def intersect(self, b):
        if not (b.r == self.r and b.c == self.c):
            raise ValueError("board %s does not have the same dimensions as self" % (b,))
        else:
            return Board(self.r,self.c, self.entries.intersection(b.entries))

    def is_in_board(self, t):
        return Set([t]).issubset(self.entries)

    def is_dist_in_board(self, t):
        if not self.is_in_board(t):
            raise ValueError("%s is not even in the board" % (t,))
        else:
            return Set([t]).issubset(self.dist_entries)

    def is_other_dist_in_board(self,t):
        if not self.is_in_board(t):
            raise ValueError("%s is not even in the board" % (t,))
        else:
            return Set([t]).issubset(self.other_dist_entries)

    def add_dist_entry(self, t):
        if not self.is_in_board(t):
            raise ValueError("%s is not even in the board" % (t,))
        else:
            self.dist_entries = self.dist_entries.union(Set([t]))

    def add_dist_entries(self, ents):
        for ent in ents:
            self.add_dist_entry(ent)

    def sizeBoard(self):
        return self.entries.cardinality()

    def print2Dtext(self):
        mat = self.diagram_matrix()
        string1 = "+"
        s = ""
        for i in range(self.c):
            string1 += "---+"
        s += string1
        s += "\n"
        for i in range(self.r):
            rowstring="|"
            for j in range(self.c):
                if mat[i][j]==0:
                    rowstring+= "   |"
                elif mat[i][j]==1:
                        rowstring+= "||||"
                elif mat[i][j]==2:
                        rowstring+= "|X||"
                elif mat[i][j]==3:
                        rowstring+= "|P||"
            s += rowstring
            s += "\n"
            s += string1
            if i < self.r - 1:
                s += "\n"
        return s

    #AsciiArt(self.print2Dtext().splitlines())
    # sage 5.11beta0

    def get_hooks(self):
        HKS = []
        for u in self.dist_entries:
            i,j = u[0], u[1]
            HKS.append(self.lam[i]+self.lamc[j]-j-i-1)
        return HKS

    def get_hooks_complement(self):
        HKS = []
        for u in self.entries.symmetric_difference(self.dist_entries):
            i,j = u[0], u[1]
            HKS.append(self.lam[i]+self.lamc[j]-j-i-1)
        return HKS

    def get_RPPhooks(self):
        HKS = []
        for u in self.other_dist_entries:
            i,j = u[0], u[1]
            HKS.append(self.lam[i]+self.lamc[j]-j-i-1)
        return HKS


    def diagram_matrix(self):
        def mat_ent(i,j):
            entry = 0
            if self.is_in_board((i,j)):
                entry += 1
                if self.is_dist_in_board((i,j)):
                    entry += 1
                if self.is_other_dist_in_board((i,j)):
                    entry += 2
            return entry
        return matrix(QQ, self.r, self.c, lambda i, j: mat_ent(i,j))

    def diagram_matrix2(self, lam):
        def mat_ent(i,j):
            entry = 2
            if self.is_in_board((i,j)):
                entry -= 1
                if self.is_dist_in_board((i,j)):
                    entry -= 1
            return entry
        return Tableau([[mat_ent(i,j) for j in range(lam[i])] for i in range(len(lam))])



    def excited_move(self, u):
        new_entries = self.dist_entries
        new_entries = new_entries.symmetric_difference(Set({u}))
        new_entries = new_entries.union(Set({(u[0]+1,u[1]+1)}))
        return Board(self.r, self.c, self.lam, self.entries, dist_entries = new_entries)


    def excited_peak_move(self, u):
        new_entries = self.dist_entries
        peaks = self.other_dist_entries
        new_entries = new_entries.symmetric_difference(Set({u}))
        new_entries = new_entries.union(Set({(u[0]+1,u[1]+1)}))
        peaks = peaks.symmetric_difference(peaks.intersection(Set({(u[0],u[1]+1),(u[0]+1,u[1])})))
        peaks = peaks.union(Set({u}))
        return Board(self.r, self.c, self.lam, self.entries, dist_entries = new_entries, other_dist_entries= peaks)

    def get_active_squares(self):
        active_squares = []
        for u in self.dist_entries:
            if self.is_in_board((u[0]+1,u[1])) and self.is_in_board((u[0], u[1]+1)) and self.is_in_board((u[0]+1, u[1]+1)):
                a = (u[0]+1,u[1])
                b = (u[0], u[1]+1)
                c = (u[0]+1,u[1]+1)
                if not(self.is_dist_in_board(a)) and not(self.is_dist_in_board(b)) and not(self.is_dist_in_board(c)):
                    active_squares.append(u)
        return active_squares

import Queue
def excited_diagrams(lam ,mu):
    b = skew_Young_diagram(lam, mu)
    ED = []
    seenED = [b]
    Q = Queue.Queue()
    Q.put(b)
    while not Q.empty():
        curr = Q.get()
        active_sq = curr.get_active_squares()
        for u in active_sq:
            new_b = curr.excited_move(u)
            if not(new_b in Set(seenED)):
                seenED.append(new_b)
                Q.put(new_b)
        if not (curr in Set(ED)):
            ED.append(curr)
    return ED

def excited_peak_diagrams(lam, mu):
    b = skew_Young_diagram(lam ,mu)
    ED = []
    seenED = [b]
    Q = Queue.Queue()
    Q.put(b)
    while not Q.empty():
        curr = Q.get()
        active_sq = curr.get_active_squares()
        for u in active_sq:
            new_b = curr.excited_peak_move(u)
            if not(new_b in Set(seenED)):
                seenED.append(new_b)
                Q.put(new_b)
        if not (curr in Set(ED)):
            ED.append(curr)
    return ED

def pleasant_diagrams(ptnlam,ptnmu, defn = 2):
    if defn == 1:
        return pleasant_diagrams_def1(ptnlam, ptnmu)
    if defn == 2:
        return pleasant_diagrams_def2(ptnlam,ptnmu)

def pleasant_diagrams_def1(ptnlam, ptnmu):
    r = len(ptnlam)
    c = ptnlam[0]
    ptnlament = Young_diagram(r,c, ptnlam).entries
    AD = []
    for ss in ptnlament.subsets():
        b = Board(r,c, ptnlam, ptnlament, dist_entries=ss)
        ab = b.diagram_matrix2(ptnlam)
        pp = IHG(ab,[0,]*len(ptnlam))
        if isSkewSupp(pp,ptnmu):
            AD.append(b)
    return AD

def pleasant_diagrams_def2(ptnlam, ptnmu):
    r = len(ptnlam)
    c = ptnlam[0]
    ptnlament = Young_diagram(r,c, ptnlam).entries
    ED = excited_diagrams(ptnlam, ptnmu)
    AD = []
    for d in ED:
        for ss in d.getNonDistEntries().subsets():
            AD.append(ss)
    SAD = set(AD)
    return [Board(r,c, ptnlam, ptnlament, dist_entries=ss) for ss in SAD]



def num_pleasant_diagrams(ptnlam, ptnmu, v=2):
    ED = excited_peak_diagrams(ptnlam, ptnmu)
    tot = 0
    for d in ED:
        stat = add(ptnlam)-add(ptnmu) - (d.other_dist_entries).cardinality()
        tot += v^stat
    return tot




def isSkewSupp(pp,ptnmu):
    for i in range(len(ptnmu)):
        for j in range(ptnmu[i]):
            if pp[i][j] <> 0:
                return False
    return True



def getAllBoard(r, c):
    return Set([(i,j) for i in range(r) for j in range(c)])


def Young_diagram(r, c, ptn, complement=False):
    entries = []
    for j in xrange(len(ptn)):
        entries.extend([(j,i) for i in xrange(ptn[j])])
    board = Board(r, c,ptn, Set(entries))
    if complement:
        board = board.getComp()
    return board



def skew_Young_diagram(ptnlam, ptnmu):
    entries = []
    r = len(ptnlam)
    c = ptnlam[0]
    ptnlament = Young_diagram(r,c, ptnlam).entries
    ptnmuent = Young_diagram(r,c, ptnmu).entries
    if ptnmuent.issubset(ptnlament):
        board = Board(r,c, ptnlam, ptnlament, dist_entries=ptnmuent)
    return board

def NHLF(lam, mu, verbosity = 0):
    ans = 0
    str = ""
    n = sum(lam)-sum(mu)
    ED = excited_diagrams(lam, mu)
    for ed in ED:
        ans += 1/mul([h for h in ed.get_hooks_complement()])
        if verbosity == 1:
            print ed, "hooks: ", ed.get_hooks_complement()
            print
    return factorial(n)*ans

def a(D, lam):
    conj_lam = Partition(lam).conjugate()
    stat = 0
    for ent in D.entries.symmetric_difference(D.dist_entries):
        i, j = ent[0], ent[1]
        stat += conj_lam[j] - i -1
    return stat

def qExcitedSSYT(lam,mu):
    ans = 0
    n = sum(lam)
    m = sum(mu)
    ED = excited_diagrams(lam, mu)
    for ed in ED:
        ans += x^a(ed, lam)/mul([qnum(i) for i in ed.get_hooks_complement()])
    return ans

def qPleasantRPP(lam,mu):
    PD = pleasant_diagrams(lam, mu)
    tot = 0
    for S in PD:
        hks = S.get_hooks()
        tot+= x^add(hks)/prod(1-x^h for h in hks)
    return tot


def qExcitedRPP(lam,mu):
    ans = 0
    n = sum(lam)
    m = sum(mu)
    EDs = excited_peak_diagrams(lam, mu)
    for ed in EDs:
        ans += x^add(ed.get_RPPhooks())/mul([qnum(i) for i in ed.get_hooks_complement()])
    return ans


def qnum(n):
    return 1-x^n

def qfact(n):
    return mul([qnum(i) for i in range(1,n+1)])
Error in lines 155-155 Traceback (most recent call last): File "/cocalc/lib/python3.11/site-packages/smc_sagews/sage_server.py", line 1250, in execute exec( File "", line 1, in <module> ModuleNotFoundError: No module named 'Queue' 
# Code Block 2: generating series of skew SSYT and RPP using P-partitions
def mydes(t):
    descents = []
    #whatpart gives the number for which t is a partition
    whatpart = sum(i for i in t.shape()[0])-sum(i for i in t.shape()[1])
    #now find the descents
    for i in range(1, whatpart):
        #find out what row i and i+1 are in (we're using the
        #standardness of self here)
        for j in range(len(t)):
            if t[j].count(i+1) > 0:
                break
            if t[j].count(i) > 0:
                descents.append(i)
                break
    return descents

def mymaj(t):
    return sum(mydes(t))


def gs_st(lam):
    Ts = StandardTableaux(lam)
    return sum([x^t.standard_major_index() for t in Ts])

def gs(lam,mu):
    Ts = StandardSkewTableaux([lam,mu])
    ans = 0
    for t in Ts:
        #print t, myword(t), mydes(t)
        ans += x^(mymaj(t))
    return ans

# column strict, i.e. SSYT
def gs_cs_rp(lam,mu):
    return gs(lam,mu)/qfact(add(lam)-add(mu))

# skew RPP
def gs_rp(lam,mu):
    P = SkewPartition([lam,mu]).cell_poset()
    n = P.cardinality()
    P2 = P.relabel({c:n-i for (i,c) in enumerate(P)})
    P3 = Poset((P2.list(),[(r[1],r[0]) for r in P2.relations()]))
    return getFpoly(P3)

def getFpoly(P):
    Ls = P.linear_extensions()
    n = P.cardinality()
   # print [L for L in Ls]
    return add([x^Permutation(le).major_index() for le in Ls])/mul([(1-x^i) for i in range(1,n+1)])



# Code Block 3: Hillman Grassl corresondence
def HG(T,mu = None):
    if not(mu):
        mu = [0,]*len(T)
    NT = [pt[:] for pt in T]
    HD = [[0,]*len(pt) for pt in T]
    ddone = False
    while not(ddone):
        if findHookPath(NT,mu):
            NT, HD = HG_step(NT,HD,mu)
        else:
            ddone = True
    return Tableau(HD)


def ppHG(T,smptn = None):
    SkewTableau(T).pp()
    print
    HG(T,mu=smptn).pp()

def HG_step(T,HD,mu):
    NT = [pt[:] for pt in T]
    NHD = [pt[:] for pt in HD]
    p = findHookPath(T,mu)
   # print "path",p
    if p:
        for e in p:
            NT[e[0]][e[1]] -= 1
        p_st = p[0]
        p_end = p[-1]
       # print "add hook", p_end[0],p_st[1]
        NHD[p_end[0]][p_st[1]] += 1
        return NT, NHD
    else:
        return None

def findHookPath(T,mu):
    st = findStartPath(T,mu)
    P = [st]
    curr = st
    if curr:
        while can_go_above(T,curr) or len(T[curr[0]])-1 > curr[1]:
           # print "LOOP", can_go_above(T,curr), len(T[curr[0]])-1, curr[1]
           # print "CURR", curr
            if can_go_above(T,curr):
                curr = [curr[0]-1,curr[1]]
                P.extend([curr])
            else:
                curr = [curr[0],curr[1]+1]
                P.extend([curr])
        return P
    else:
        return None

def can_go_above(T,ent):
   # print T, ent
    if ent[0] == 0:
        return False
    else:
        val = T[ent[0]][ent[1]]
        val_above = T[ent[0]-1][ent[1]]
        if val_above:
            if val == val_above:
                return True
            else:
                return False
        else:
            return False

def findStartPath(T, mu):
    # find lowest nonzero row
    curr_r = len(T)-1
    curr_c = mu[curr_r]
    ddone = false
    while not(ddone):
      #  print curr_r,curr_c, "length",len(T[curr_r])-1, curr_c,len(T[curr_r])-1 > curr_c
        if T[curr_r][curr_c] == 0:
            if len(T[curr_r])-1 > curr_c:
                curr_c += 1
            elif curr_r > 0:
                curr_r -= 1
                curr_c = mu[curr_r]
            else:
                return None
        else:
            ddone = true
    if ddone:
        return [curr_r,curr_c]

def IHG(HD,mu=None):
    if not(mu):
        mu = [0,]*len(HD)
    orderhks = getOrderHooks(HD)
    lam = [len(D) for D in HD]
    T = [[None,]*mu[i]+[0,]*(lam[i]-mu[i]) for i in range(len(lam))]
    for hk in orderhks:
        T = Rev_HG_step(T,hk,lam,mu)
    return SkewTableau(T)

def ppIHG(HD,smptn):
    lam = [len(D) for D in HD]
    SkewTableau(HD).pp()
    T = IHG(HD,mu=smptn)
    print
    for r in range(1,len(lam)):
        print r, min([T[r][j]-T[r-1][j] for j in range(smptn[r-1],lam[r])])
    T.pp()

def ppIHGstopcol(HD,mu,stpc=[]):
    lam = [len(D) for D in HD]
    SkewTableau(HD).pp()
    HDs = []
    for c in stpc:
        cHD = cHD = [[0,]*c+[HD[i][j] for j in range(c,len(HD[i]))] for i in range(len(HD))]
        print "column ", c
        ppIHG(cHD,mu)

def ppIHGstopcol(HD,mu,stpc=[],jc=[]):
    lam = [len(D) for D in HD]
    SkewTableau(HD).pp()
    HDs = []
    for c in stpc:
        cHD = [[0,]*c+[HD[i][j] for j in range(c,len(HD[i]))] for i in range(len(HD))]
        print "column ", c
        SkewTableau(cHD).pp()
        T = IHG(cHD,mu)
        print
        if len(jc) >0:
            for r in range(1,len(lam)-1):
                print r, min([T[r][j]-T[r-1][j] for j in range(jc[r],lam[r])])
        T.pp()

def getOrderHooks(HD):
    L = []
    lam = [len(D) for D in HD]
    lamc = Partition(lam).conjugate()
    for i in range(len(lamc)):
        for j in range(lamc[len(lamc)-i-1]):
            for k in range(HD[j][len(lamc)-i-1]):
                L.append([j,len(lamc)-i-1])
    return L


def Rev_HG_step(T,hk,lam,mu):
    NT = [pt[:] for pt in T]
    p = findRevHookPath(T,hk,lam,mu)
   # print "path",p
    for e in p:
        NT[e[0]][e[1]] += 1
    return NT

def findRevHookPath(T,hk,lam,mu):
    st = [hk[0],lam[hk[0]]-1]
    P = [st]
    curr = st
    lamc = Partition(lam).conjugate()
    while can_go_down(T,curr,lamc) or curr[1] > max(mu[curr[0]],hk[1]):
        if can_go_down(T,curr,lamc):
            curr = [curr[0]+1,curr[1]]
            P.extend([curr])
        else:
            curr = [curr[0],curr[1]-1]
            P.extend([curr])
    return P


def can_go_down(T,ent,lamc):
    if ent[0] == lamc[ent[1]]-1:
        return False
    else:
        val = T[ent[0]][ent[1]]
        val_down = T[ent[0]+1][ent[1]]
        if val == val_down:
            return True
        else:
            return False
Error in lines 84-91 Traceback (most recent call last): File "/cocalc/lib/python3.11/site-packages/smc_sagews/sage_server.py", line 1251, in execute compile(block + '\n', File "<string>", line 7 print r, min([T[r][j]-T[r-Integer(1)][j] for j in range(smptn[r-Integer(1)],lam[r])]) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ SyntaxError: Missing parentheses in call to 'print'. Did you mean print(...)? 


# SECTION 3
# Example 3.1
ED = excited_diagrams([2,2,2,1], [1,1,0,0])
for e in ED:
    e
    print

NHLF([5,5,5,5,5],[2,2,0,0,0])
69283500 
NHLF([2,2,2,1],[1,1,0,0])
9 
# verify NHLF answer
StandardSkewTableaux([[2,2,2,1],[1,1,0,0]]).cardinality()
9 
NHLF([2,2,2,1],[1,1,0,0],verbosity=1)
+---+---+ ||X|||||| +---+---+ ||X|||||| +---+---+ ||||||||| +---+---+ ||||| | +---+---+ hooks: [3, 1, 2, 1, 3] +---+---+ ||X|||||| +---+---+ ||||||||| +---+---+ ||||||X|| +---+---+ ||||| | +---+---+ hooks: [3, 1, 4, 2, 3] +---+---+ ||||||||| +---+---+ ||||||X|| +---+---+ ||||||X|| +---+---+ ||||| | +---+---+ hooks: [3, 1, 4, 5, 3] 9 




# Example 3.2
qExcitedSSYT([2,2,2,1],[1,1,0,0]).show()
−x8(x5−1)(x4−1)(x3−1)2(x−1)−x6(x4−1)(x3−1)2(x2−1)(x−1)−x4(x3−1)2(x2−1)(x−1)2\displaystyle -\frac{x^{8}}{{\left(x^{5} - 1\right)} {\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)}^{2} {\left(x - 1\right)}} - \frac{x^{6}}{{\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)} {\left(x - 1\right)}} - \frac{x^{4}}{{\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)} {\left(x - 1\right)}^{2}}
# verify we get generating series of SSYT
qExcitedSSYT([2,2,2,1],[1,1,0,0]).series(x,15)
1*x^4 + 2*x^5 + 5*x^6 + 9*x^7 + 16*x^8 + 25*x^9 + 38*x^10 + 55*x^11 + 78*x^12 + 106*x^13 + 142*x^14 + Order(x^15) 
gs_cs_rp([2,2,2,1],[1,1,0,0]).series(x,15)
1*x^4 + 2*x^5 + 5*x^6 + 9*x^7 + 16*x^8 + 25*x^9 + 38*x^10 + 55*x^11 + 78*x^12 + 106*x^13 + 142*x^14 + Order(x^15) 
# Section 3.2: Flagged tableaux

def getFlaggedTableaux(mu,flag):
    def is_max_ent(rpp,lam, maxe):
        isgood = True
        i = 0
        while isgood and i<len(lam):
            if rpp[i][lam[i]-1] > maxe:
                isgood = False
                break
            i+=1
        return isgood
    Ts = SemistandardTableaux(mu,max_entry=max(flag))
    gT = []
    for T in Ts:
        isgood = True
        i = 0
        while isgood and i<len(mu):
            if T[i][mu[i]-1] > flag[i]:
                isgood = False
                break
            i+=1
        if isgood:
            gT.append(T)
    return gT

def numFlaggedTableaux(mu,flag):
    k = len(mu)
    M = matrix(k,k,lambda i, j: binomial(flag[i]+mu[i]-i+j-1,mu[i]-i+j))
    print M
    return M.determinant();



# Figure 1
FT = getFlaggedTableaux([1,1],[3,3])
for T in FT:
    T.pp()
    print
1 2 1 3 2 3 
numFlaggedTableaux([2,2],[5,6])
[15 35] [ 6 21] 105 
# Example Figure 2 and Corollary 3.7
skew_Young_diagram([14,12,12,10,10,10,10,8,8,4],[7,4,4,4,2])
+---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||X|||X|||X|||X|||X|||X|||X|||||||||||||||||||||||||||||| +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||X|||X|||X|||X|||||||||||||||||||||||||||||||||| | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||X|||X|||X|||X|||||||||||||||||||||||||||||||||| | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||X|||X|||X|||X|||||||||||||||||||||||||| | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||X|||X|||||||||||||||||||||||||||||||||| | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||||||||||||||||||||||||||||||||||||||||| | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||||||||||||||||||||||||||||||||||||||||| | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||||||||||||||||||||||||||||||||| | | | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||||||||||||||||||||||||||||||||| | | | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ ||||||||||||||||| | | | | | | | | | | +---+---+---+---+---+---+---+---+---+---+---+---+---+---+ 


numFlaggedTableaux([7,4,4,4,2],[4,8,8,8,9])
[ 120 165 220 286 364] [ 120 330 792 1716 3432] [ 36 120 330 792 1716] [ 8 36 120 330 792] [ 0 0 1 9 45] 7932960 




# SECTION 5: The Hillman-Grassl and the RSK correspondence

# first figure
HG([[0,0],[0,1],[0,2],[1]]).pp()
0 0 0 1 0 1 1 
HG([[0,0],[0,1],[1,2],[2]]).pp()
0 0 0 1 1 0 1 
HG([[0,0],[0,1],[2,2],[3]]).pp()
0 0 1 0 1 0 1 
IHG([[0,0],[0,1],[0,1],[1]]).pp()
0 0 0 1 0 2 1 
IHG([[0,0],[0,1],[1,0],[1]]).pp()
0 0 0 1 1 2 2 
IHG([[0,0],[1,0],[1,0],[1]]).pp()
0 0 0 1 2 2 3 
# Example 5.9
pi = [[0,1,3,4],[1,3,5,6],[3,6,7],[3]]
Tableau(pi).pp()
A = HG(pi)
A.pp()
0 1 3 4 1 3 5 6 3 6 7 3 0 2 1 1 1 1 1 2 2 1 1 0 
A1vf = [[1,2,0],[1,1,1]]
Tableau(A1vf).pp()
print
RSK(A1vf)[0].pp()
1 2 0 1 1 1 1 1 2 2 3 2 
A0vf = [[1,2,0],[1,1,1],[1,1,2]]
Tableau(A0vf).pp()
print
RSK(A0vf)[0].pp()
1 2 0 1 1 1 1 1 2 1 1 1 2 2 3 3 2 2 3 
# SECTION 6: Hillman-Grassl map on Skew RPP

# Section 6.1: Pleasant diagrams

# Example 6.3
PD = pleasant_diagrams([2,2],[1,0], defn = 1)
len(PD)
for S in PD:
    S
    print
12 +---+---+ ||X|||||| +---+---+ ||||||||| +---+---+ +---+---+ ||||||||| +---+---+ ||||||X|| +---+---+ +---+---+ ||X|||X|| +---+---+ ||||||||| +---+---+ +---+---+ ||||||X|| +---+---+ ||||||X|| +---+---+ +---+---+ ||X|||||| +---+---+ ||X|||||| +---+---+ +---+---+ ||||||||| +---+---+ ||X|||X|| +---+---+ +---+---+ ||X|||||| +---+---+ ||||||X|| +---+---+ +---+---+ ||X|||X|| +---+---+ ||X|||||| +---+---+ +---+---+ ||||||X|| +---+---+ ||X|||X|| +---+---+ +---+---+ ||X|||X|| +---+---+ ||||||X|| +---+---+ +---+---+ ||X|||||| +---+---+ ||X|||X|| +---+---+ +---+---+ ||X|||X|| +---+---+ ||X|||X|| +---+---+ 
# Example of Theorem 6.4
qPleasantRPP([2,2],[1,0]).show()
−x7(x3−1)(x2−1)2+2 x5(x3−1)(x2−1)−x3x3−1+x4(x2−1)2−x5(x2−1)2(x−1)−2 x2x2−1+2 x3(x2−1)(x−1)−xx−1+1\displaystyle -\frac{x^{7}}{{\left(x^{3} - 1\right)} {\left(x^{2} - 1\right)}^{2}} + \frac{2 \, x^{5}}{{\left(x^{3} - 1\right)} {\left(x^{2} - 1\right)}} - \frac{x^{3}}{x^{3} - 1} + \frac{x^{4}}{{\left(x^{2} - 1\right)}^{2}} - \frac{x^{5}}{{\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}} - \frac{2 \, x^{2}}{x^{2} - 1} + \frac{2 \, x^{3}}{{\left(x^{2} - 1\right)} {\left(x - 1\right)}} - \frac{x}{x - 1} + 1
# verify we get generating series of RPP

qPleasantRPP([2,2],[1,0]).series(x,15)
1 + 1*x + 3*x^2 + 4*x^3 + 6*x^4 + 8*x^5 + 11*x^6 + 13*x^7 + 17*x^8 + 20*x^9 + 24*x^10 + 28*x^11 + 33*x^12 + 37*x^13 + 43*x^14 + Order(x^15) 
gs_rp([2,2],[1,0]).series(x,15)
1 + 1*x + 3*x^2 + 4*x^3 + 6*x^4 + 8*x^5 + 11*x^6 + 13*x^7 + 17*x^8 + 20*x^9 + 24*x^10 + 28*x^11 + 33*x^12 + 37*x^13 + 43*x^14 + Order(x^15) 


# Example of Theorem 6.11, characterization of pleasant diagrams

# defn = 1 using inverse HG
len(pleasant_diagrams([2,2],[1,0], defn = 1))
12 
# defn = 2 using subsets of complements of excited diagrams
len(pleasant_diagrams([2,2],[1,0], defn = 2))
12 
# Example of Theorem 6.15: number of pleasant diagrams [4,4,4]/[2]
num_pleasant_diagrams([4,4,4],[2,0,0])
2816 
len(pleasant_diagrams([4,4,4],[2,0,0]))
2816 
# slow
len(pleasant_diagrams([4,4,4],[2,0,0], defn = 1))
2816 
# Example of Corollary 6.10: generating series of skew RPP in terms of excited diagrams
# Figure 7
qExcitedRPP([4,4,4], [2,0,0]).show()
x10(x6−1)(x5−1)2(x4−1)3(x3−1)3(x2−1)+x9(x6−1)(x5−1)2(x4−1)2(x3−1)3(x2−1)2+x8(x5−1)2(x4−1)3(x3−1)3(x2−1)2+x6(x6−1)(x5−1)2(x4−1)2(x3−1)2(x2−1)2(x−1)+x5(x5−1)2(x4−1)3(x3−1)2(x2−1)2(x−1)+1(x5−1)(x4−1)3(x3−1)3(x2−1)2(x−1)\displaystyle \frac{x^{10}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{3} {\left(x^{3} - 1\right)}^{3} {\left(x^{2} - 1\right)}} + \frac{x^{9}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{2} {\left(x^{3} - 1\right)}^{3} {\left(x^{2} - 1\right)}^{2}} + \frac{x^{8}}{{\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{3} {\left(x^{3} - 1\right)}^{3} {\left(x^{2} - 1\right)}^{2}} + \frac{x^{6}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{2} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}} + \frac{x^{5}}{{\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{3} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}} + \frac{1}{{\left(x^{5} - 1\right)} {\left(x^{4} - 1\right)}^{3} {\left(x^{3} - 1\right)}^{3} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}}
# verify we get generating series of RPP
qExcitedRPP([4,4,4], [2,0,0]).series(x,15)
1 + 1*x + 3*x^2 + 6*x^3 + 12*x^4 + 20*x^5 + 37*x^6 + 59*x^7 + 99*x^8 + 153*x^9 + 239*x^10 + 356*x^11 + 534*x^12 + 769*x^13 + 1109*x^14 + Order(x^15) 
gs_rp([4,4,4],[2,0,0]).series(x,15)
1 + 1*x + 3*x^2 + 6*x^3 + 12*x^4 + 20*x^5 + 37*x^6 + 59*x^7 + 99*x^8 + 153*x^9 + 239*x^10 + 356*x^11 + 534*x^12 + 769*x^13 + 1109*x^14 + Order(x^15) 


# SECTION 7: Hillman-Grassl map on skew SSYT

# Example 7.9
qExcitedSSYT([4,4,4,2], [3,1,0,0]).show()
x13(x7−1)(x6−1)2(x5−1)(x4−1)(x3−1)(x2−1)2(x−1)2+x12(x6−1)2(x5−1)2(x4−1)(x3−1)(x2−1)2(x−1)2+x11(x6−1)(x5−1)2(x4−1)2(x3−1)(x2−1)2(x−1)2+x10(x6−1)(x5−1)2(x4−1)(x3−1)2(x2−1)2(x−1)2+x9(x6−1)(x5−1)2(x3−1)2(x2−1)3(x−1)2+x9(x5−1)2(x4−1)2(x3−1)2(x2−1)2(x−1)2+x8(x5−1)2(x4−1)(x3−1)2(x2−1)3(x−1)2\displaystyle \frac{x^{13}}{{\left(x^{7} - 1\right)} {\left(x^{6} - 1\right)}^{2} {\left(x^{5} - 1\right)} {\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}^{2}} + \frac{x^{12}}{{\left(x^{6} - 1\right)}^{2} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}^{2}} + \frac{x^{11}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{2} {\left(x^{3} - 1\right)} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}^{2}} + \frac{x^{10}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}^{2}} + \frac{x^{9}}{{\left(x^{6} - 1\right)} {\left(x^{5} - 1\right)}^{2} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{3} {\left(x - 1\right)}^{2}} + \frac{x^{9}}{{\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)}^{2} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{2} {\left(x - 1\right)}^{2}} + \frac{x^{8}}{{\left(x^{5} - 1\right)}^{2} {\left(x^{4} - 1\right)} {\left(x^{3} - 1\right)}^{2} {\left(x^{2} - 1\right)}^{3} {\left(x - 1\right)}^{2}}
# verify we get generating series of SSYT
qExcitedSSYT([4,4,4,2], [3,1,0,0]).series(x,20)
1*x^8 + 4*x^9 + 11*x^10 + 26*x^11 + 55*x^12 + 108*x^13 + 198*x^14 + 347*x^15 + 582*x^16 + 946*x^17 + 1490*x^18 + 2291*x^19 + Order(x^20) 
gs_cs_rp([4,4,4,2], [3,1,0,0]).series(x,20)
1*x^8 + 4*x^9 + 11*x^10 + 26*x^11 + 55*x^12 + 108*x^13 + 198*x^14 + 347*x^15 + 582*x^16 + 946*x^17 + 1490*x^18 + 2291*x^19 + Order(x^20) 
# SECTION 8: Excited diagrams and SSYT of thick strips delta_n+2k/delta_n

# Figure 12
ED = excited_diagrams([4,3,2,1], [2,1,0,0])
for D in ED:
    D
    print
+---+---+---+---+ ||X|||X|||||||||| +---+---+---+---+ ||X|||||||||| | +---+---+---+---+ ||||||||| | | +---+---+---+---+ ||||| | | | +---+---+---+---+ +---+---+---+---+ ||X|||||||||||||| +---+---+---+---+ ||X|||||||X|| | +---+---+---+---+ ||||||||| | | +---+---+---+---+ ||||| | | | +---+---+---+---+ +---+---+---+---+ ||X|||X|||||||||| +---+---+---+---+ ||||||||||||| | +---+---+---+---+ ||||||X|| | | +---+---+---+---+ ||||| | | | +---+---+---+---+ +---+---+---+---+ ||X|||||||||||||| +---+---+---+---+ ||||||||||X|| | +---+---+---+---+ ||||||X|| | | +---+---+---+---+ ||||| | | | +---+---+---+---+ +---+---+---+---+ ||||||||||||||||| +---+---+---+---+ ||||||X|||X|| | +---+---+---+---+ ||||||X|| | | +---+---+---+---+ ||||| | | | +---+---+---+---+ 
FT = getFlaggedTableaux([2,1],[2,3])
for T in FT:
    T.pp()
    print
1 1 2 1 1 3 1 2 2 1 2 3 2 2 3 
# Example 8.16
NHLF([3,2,1],[1,0,0],verbosity=1)
+---+---+---+ ||X|||||||||| +---+---+---+ ||||||||| | +---+---+---+ ||||| | | +---+---+---+ hooks: [3, 1, 3, 1, 1] +---+---+---+ ||||||||||||| +---+---+---+ ||||||X|| | +---+---+---+ ||||| | | +---+---+---+ hooks: [3, 1, 3, 5, 1] 16 
# Example Corollary 8.2
[len(excited_diagrams([2,1],[0,0])),len(excited_diagrams([3,2,1],[1,0,0])),len(excited_diagrams([4,3,2,1],[2,1,0,0])),len(excited_diagrams([5,4,3,2,1],[3,2,1,0,0]))]
[1, 2, 5, 14] 

# SECTION 9: Pleasant diagrams and SSYT of thick strips delta_n+2k/delta_n



# Example Theorem 9.1

[num_pleasant_diagrams([2,1],[0,0]),num_pleasant_diagrams([3,2,1],[1,0,0]),num_pleasant_diagrams([4,3,2,1],[2,1,0,0]),num_pleasant_diagrams([5,4,3,2,1],[3,2,1,0,0])]
[8, 48, 352, 2880] 
# Schroder numbers
sch = [1, 1, 3, 11, 45, 197, 903, 4279, 20793, 103049, 518859, 2646723, 13648869, 71039373, 372693519, 1968801519, 10463578353, 55909013009, 300159426963, 1618362158587, 8759309660445, 47574827600981, 259215937709463, 1416461675464871]



[2^3*sch[1], 2^4*sch[2], 2^5*sch[3], 2^6*sch[4]]
[8, 48, 352, 2880] 


# Conjecture 9.3: number of pleasant diagrams of delta_n+2k/delta_n
def Schdetk(n,k):
    M = matrix(QQ,k,k,lambda i, j: plsch(n+i+j))
    print M
    return M.determinant()*2^(k*(k-1)/2);

def plsch(n):
    return 2^(n+2)*sch[n]

num_pleasant_diagrams([5,4,3,2,1],[1,0,0,0,0])
28672 
Schdetk(2,2)
[ 48 352] [ 352 2880] 28672 
num_pleasant_diagrams([7,6,5,4,3,2,1],[1,0,0,0,0,0,0])
251658240 
Schdetk(2,3)
[ 48 352 2880] [ 352 2880 25216] [ 2880 25216 231168] 251658240 
num_pleasant_diagrams([6,5,4,3,2,1],[2,1,0,0,0,0])
1163264 
Schdetk(3,2)
[ 352 2880] [ 2880 25216] 1163264 
num_pleasant_diagrams([8,7,6,5,4,3,2,1],[2,1,0,0,0,0,0,0])
49257906176 
Schdetk(3,3)
[ 352 2880 25216] [ 2880 25216 231168] [ 25216 231168 2190848] 49257906176 
# Conjecture 9.5: determinental identity of generating series of RPP of shape delta_n+2k/delta_n

# normalized qEuler numbers \tilde{E}^*_(2m+1)(q) computed using Stembridge's poset package.

odd_nqE = {1:(x+1)/((1-x)*(1-x^2)*(1-x^3)), 2:(x^7+2*x^6+2*x^5+3*x^4+3*x^3+2*x^2+2*x+1)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)), 3:(x^17+3*x^16+5*x^15+9*x^14+14*x^13+19*x^12+25*x^11+29*x^10+31*x^9+31*x^8+29*x^7+25*x^6+19*x^5+14*x^4+9*x^3+5*x^2+3*x+1)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)*(1-x^6)*(1-x^7)),4:(x^31+4*x^30+9*x^29+19*x^28+36*x^27+61*x^26+98*x^25+147*x^24+207*x^23+277*x^22+354*x^21+432*x^20+505*x^19+568*x^18+613*x^17+637*x^16+637*x^15+613*x^14+568*x^13+505*x^12+432*x^11+354*x^10+277*x^9+207*x^8+147*x^7+98*x^6+61*x^5+36*x^4+19*x^3+9*x^2+4*x+1)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)*(1-x^6)*(1-x^7)*(1-x^8)*(1-x^9)),5:(1+5*x+x^49+14*x^2+74*x^4+34*x^3+459*x^7+266*x^6+145*x^5+1167*x^9+749*x^8+3479*x^12+2507*x^11+1744*x^10+6094*x^14+4675*x^13+11437*x^17+9519*x^16+7720*x^15+15332*x^19+13402*x^18+19991*x^22+18716*x^21+17135*x^20+21346*x^24+20885*x^23+19991*x^27+20885*x^26+21346*x^25+17135*x^29+18716*x^28+13402*x^31+15332*x^30+7720*x^34+9519*x^33+11437*x^32+4675*x^36+6094*x^35+1744*x^39+2507*x^38+3479*x^37+749*x^41+1167*x^40+266*x^43+459*x^42+34*x^46+74*x^45+145*x^44+5*x^48+14*x^47)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)*(1-x^6)*(1-x^7)*(1-x^8)*(1-x^9)*(1-x^10)*(1-x^11)),6:(x^71+6*x^70+20*x^69+55*x^68+134*x^67+294*x^66+599*x^65+1147*x^64+2079*x^63+3597*x^62+5973*x^61+9557*x^60+14787*x^59+22191*x^58+32377*x^57+46022*x^56+63849*x^55+86584*x^54+114916*x^53+149446*x^52+190618*x^51+238670*x^50+293573*x^49+354975*x^48+422183*x^47+494135*x^46+569408*x^45+646256*x^44+722662*x^43+796414*x^42+865213*x^41+926786*x^40+978990*x^39+1019948*x^38+1048146*x^37+1062517*x^36+1062517*x^35+1048146*x^34+1019948*x^33+978990*x^32+926786*x^31+865213*x^30+796414*x^29+722662*x^28+646256*x^27+569408*x^26+494135*x^25+422183*x^24+354975*x^23+293573*x^22+238670*x^21+190618*x^20+149446*x^19+114916*x^18+86584*x^17+63849*x^16+46022*x^15+32377*x^14+22191*x^13+14787*x^12+9557*x^11+5973*x^10+3597*x^9+2079*x^8+1147*x^7+599*x^6+294*x^5+134*x^4+55*x^3+20*x^2+6*x+1)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)*(1-x^6)*(1-x^7)*(1-x^8)*(1-x^9)*(1-x^10)*(1-x^11)*(1-x^12)*(1-x^13))  ,7:(1+7*x+73603123*x^49+27*x^2+223*x^4+83*x^3+2500*x^7+1198*x^6+538*x^5+9279*x^9+4933*x^8+48864*x^12+29084*x^11+16738*x^10+126194*x^14+79630*x^13+433694*x^17+293978*x^16+194912*x^15+888247*x^19+626718*x^18+2274766*x^22+1690672*x^21+1236097*x^20+3932187*x^24+3013188*x^23+8038646*x^27+6419141*x^26+5058656*x^25+12138816*x^29+9939284*x^28+17474843*x^31+14649146*x^30+27752777*x^34+24045984*x^33+20611779*x^32+35829478*x^36+31696298*x^35+48741542*x^39+44423789*x^38+40094498*x^37+57010351*x^41+52965692*x^40+64215233*x^43+60788559*x^42+71626980*x^46+69701332*x^45+67210232*x^44+73603123*x^48+72938719*x^47+31696298*x^62+35829478*x^61+40094498*x^60+24045984*x^64+27752777*x^63+17474843*x^66+20611779*x^65+9939284*x^69+12138816*x^68+14649146*x^67+6419141*x^71+8038646*x^70+3013188*x^74+3932187*x^73+5058656*x^72+1690672*x^76+2274766*x^75+888247*x^78+1236097*x^77+293978*x^81+433694*x^80+626718*x^79+126194*x^83+194912*x^82+29084*x^86+48864*x^85+79630*x^84+9279*x^88+16738*x^87+2500*x^90+4933*x^89+223*x^93+538*x^92+1198*x^91+27*x^95+83*x^94+7*x^96+x^97+72938719*x^50+69701332*x^52+71626980*x^51+64215233*x^54+67210232*x^53+52965692*x^57+57010351*x^56+60788559*x^55+44423789*x^59+48741542*x^58)/((1-x)*(1-x^2)*(1-x^3)*(1-x^4)*(1-x^5)*(1-x^6)*(1-x^7)*(1-x^8)*(1-x^9)*(1-x^10)*(1-x^11)*(1-x^12)*(1-x^13)*(1-x^14)*(1-x^15))}

def qEuler2detk(n,k):
    def N(n,k):
        return k*(k-1)*(6*n+8*k-1)/6
    print matrix(k,k,lambda i, j: n+2*k-2-i-j)
    M = matrix(k,k,lambda i, j: odd_nqE[n+2*k-2-i-j])
    return M.determinant()*x^(-N(n,k));



# n=2, k=1
LHS = qExcitedRPP([3,2,1],[1,0,0])

RHS = qEuler2detk(2,1)
[2] 


(LHS/RHS).series(x,10)
1 + Order(x^10) 


# n=2, k=2
LHS = qExcitedRPP([5,4,3,2,1],[1,0,0,0,0])

RHS = qEuler2detk(2,2)
[4 3] [3 2] 
(LHS/RHS).series(x,10)
1 + Order(x^10) 
# n=3, k=2
LHS = qExcitedRPP([6,5,4,3,2,1],[2,1,0,0,0,0])

RHS = qEuler2detk(3,2)
[5 4] [4 3] 
(LHS/RHS).series(x,10)
1 + Order(x^10) 
# n=3, k=3
RHS = qEuler2detk(3,3)
[7 6 5] [6 5 4] [5 4 3] 
LHS = qExcitedRPP([8,7,6,5,4,3,2,1],[2,1,0,0,0,0,0,0])

# slow
(LHS/RHS).series(x,3)
Error in lines 2-2 Traceback (most recent call last): File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 942, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "sage/symbolic/expression.pyx", line 4000, in sage.symbolic.expression.Expression.series (/projects/sage/sage-6.10/src/build/cythonized/sage/symbolic/expression.cpp:24300) sig_on() File "sage/ext/interrupt/interrupt.pyx", line 88, in sage.ext.interrupt.interrupt.sig_raise_exception (/projects/sage/sage-6.10/src/build/cythonized/sage/ext/interrupt/interrupt.c:924) raise KeyboardInterrupt KeyboardInterrupt 

︠2fd4b727-bd5c-4188-aed9-dc0cb59c6c52︠

︠26a208de-d2e9-4608-99e3-0f4cf5ab44bf︠

︠3a31367a-2365-4a87-8a33-98dda4447991︠

︠a59a678b-3119-4a0e-a36e-fd2c1cf08d54︠

︠a12de0f8-de4a-472c-a84c-af50110c5ba3︠

︠134fc128-d2e9-4740-bc36-e3c746398b84︠