%auto
typeset_mode(True, display=False)

%md
## Supplementary Sage worksheet for the article <b>On the Okounkov-Olshanski formula for standard tableaux of skew shapes</b>, <a href="https://arxiv.org/abs/2007.05006">arxiv:2007.05006</a> Alejandro H. Morales, Daniel G. Zhu

%md 
### The Okounkov-Olshanski formula states that the number $f^{\lambda/\mu}$ of SYT of skew shape equals
### $f^{\lambda/\mu} = \frac{|\lambda/\mu|!}{\prod_{u \in [\lambda]} h(u)} \sum_{T \in SSYT(\mu,d)} \prod_{u\in [\mu]} (\lambda_{d+1-T(u)} - j+i),$
### where $d=\ell(\lambda)$. 

### Let's look at an example.

def OOF(lam,mu, verbosity=1):
    d = len(lam)
    SSYT = SemistandardTableaux(mu,max_entry=d)
    tot = 0
    for T in SSYT:
        contrib = mul([lam[d+1-T[u[0]][u[1]]-1]-u[1]+u[0] for u in T.cells()])
        if verbosity == 1:
            print "tableaux", T, "conbribution", contrib
        tot += contrib
    return tot*factorial(add(lam)-add(mu))/mul(Partition(lam).hooks())

OOF([2,2,2,1],[1,1])

StandardTableaux([[2,2,2,1],[1,1]]).cardinality()

%md 
It is not so hard to show that this formula is nonnegative. However some terms have zero contribution.

OOF([4,3,2,1],[2,1])

%md
## Number of terms ##
### Our first main result, Theorem 1.2, has two determinant formulas for the number of nonzero terms of the Okounkov-Olshanski formula that we denote by $OOT(\lambda/\mu)$

def OOT1(lam,mu):
    d = len(lam)
    mup = mu +[0,]*(d-len(mu))
    def mat_ent(i,j,lam,mu):
        if lam[i]-mu[j]+(j+1)-1 >= (i+1)-1: 
            return binomial(lam[i] - mu[j]+(j+1)-1,(i+1)-1)
        else: 
            return 0
    M = matrix(d,d,lambda i, j: mat_ent(i,j,lam,mup))
    print M
    return M.determinant();

def OOT2(lam,mu):
    lamc = Partition(lam).conjugate()
    muc = Partition(mu).conjugate()
    mucp = muc + [0,]*(len(lamc)-len(muc))
    k = mu[0]
    def mat_ent(i,j,lamc,muc):
        if lamc[i] >= muc[j] +i-j:
            return binomial(lamc[i],muc[j]+i-j)
        else:
            return 0
    M = matrix(k,k, lambda i,j: mat_ent(i,j,lamc,mucp))
    print M
    return M.determinant()
    
    

% md
### Let's do some examples

OOT1([2,2,2,1],[1,1])

OOT2([2,2,2,1],[1,1])

OOT1([3,2,1],[1])

OOT2([3,2,1],[1])

OOT1([4,3,2,1],[2,1])

OOT2([4,3,2,1],[2,1])

%md
### Our second main result is a formula for the number of terms for thick zigzags $\delta_{n+2k}/\delta_n$ as a determinant of normalized Genocchi numbers $G_{2n}/(2n)!$, In particular $OOT(\delta_{n+2}/\delta_n) = G_{2n+2}$

G = [1, 1, 3, 17, 155, 2073, 38227, 929569, 28820619, 1109652905, 51943281731, 2905151042481, 191329672483963, 14655626154768697, 1291885088448017715, 129848163681107301953, 14761446733784164001387, 1884515541728818675112649, 268463531464165471482681379]

def nG(n):
    return G[n-1]/factorial(2*n)

def staircase(n):
    return range(1,n)[::-1]

def OOT_thick_zigzag(n,k):
    M = matrix(k, k, lambda i, j: nG(n+(i+1)+(j+1)-1))
    print "shape", staircase(n+2*k),"/", staircase(n)
   # print M
    return M.determinant()*mul([factorial(2*(n+i+k-1))/factorial(2*i-1) for i in range(1,k+1)]);
    

% md
Interestingly this number of terms is proportional to the number of SYT of this shape.

n=5
k=3
print "size shape", add(staircase(n+2*k))-add(staircase(n))
factor(StandardTableaux([staircase(n+2*k),staircase(n)]).cardinality()/OOT_thick_zigzag(n,k))

%md 
## Reformulations
### We give several reformulations of the Okounkov-Olshanski formula. One is in terms of <b>reverse excited diagrams</b>. Excited diagrams appear in the Naruse hook length formula for skew shapes. The reverse excited diagram reformulation is 

### $f^{\lambda/\mu} \,=\, \frac{|\lambda/\mu|!}{\prod_{u\in [\lambda]} h(u)} \sum_{D\in \mathcal{RE}(\lambda)}  \prod_{(i,j) \in B(D)} (\lambda_i+d-j)$,

### where $d=\ell(\lambda)$, $B(D)$ are certain cells of $[\lambda]$ (viewed as a shifted skew shape) associated to $D$ and $\lambda_i+d-j$ equals the arm-length of the cell $(i,j)$.

# read text file with code
load("excited_code.sage")

OOF_RED([2,2,2,1],[1,1,0,0],verbosity=1)

%md
### Another reformulation is in terms of Knutson-Tao puzzles: $f^{\lambda/\mu} = \frac{|\lambda/\mu|!}{\prod_{u\in[\lambda]} h(u)} \sum_{P \in \Delta^{\lambda\mu}_\lambda} \prod_{p \in \lozenge(P)}\operatorname{ht}(p)$.

latex.extra_preamble(r'''\usepackage{tikz}''')
vars =  ['y%d' % i for i in range(1,7)]
R = PolynomialRing(QQ, 6, vars)
ps = KnutsonTaoPuzzleSolver("HT")
lam = [2,2,2,1]
mu = [1,1]
solns = ps('101000', '0010011')
tot = 0
for p in solns:
    if p.south_labels() == ('1','0','1','0','0','0'):
        curr = R(p.contribution()).substitute(y1=1,y2=2,y3=3,y4=4,y5=5,y6=6)
        p, curr
        tot += curr
print "expecting 9", tot*factorial(add(lam)-add(mu))/mul(Partition(lam).hooks())

%md 
# $q$-analogues for skew reverse plane partitions
## Stanley-Chen have a $q$-analogue of the Okounkov-Olshanski formula for the generating function of skew semistandard tableaux.
## Our last main results are $q$-analogues of the Okounkov-Olshanski formula for the generating function of skew reverse plane partitions. Let $rpp_{\lambda/\mu}(q) := \sum_{T \in RPP(\lambda/\mu)} q^{|T|}$.
### <b>Theorem (1st $q$-analogue):</b> $\frac{rpp_{\lambda/\mu}(q)}{rpp_{\lambda}(q)} = \sum_{T \in SSYT(\mu, d)} q^{p(T)} \prod_{u \in [\mu]} (1 - q^{w(u, T(u))})$,
### where $p(T)=\sum_{u \in [\mu], m_T(u) \leq k < T(u)} w(u,k)$ and $m_T(u)$  is the minimum $k\leq T(u)$ such that replacing $T(u)$ by $k$ still gives a semistandard tableaux.

%md 
### Example: For $\lambda/\mu = 2221/11$. From the theory of $P$-partitions or computing it directly the answer for $rpp_{\lambda/\mu}(x)/rpp_{\lambda}(x) = 1-x^5-x^6-x^7+x^8+x^{10}$ 

f1 = x^3*(1-x^2)*(1-x^3) + x^4*(1-x^2)*(1-x^3) + x*(1-x^2)*(1-x^3)+x^6*(1-x)*(1-x^3)+x^3*(1-x)*(1-x^3)+x^0*(1-x^1)*(1-x^3)
f1

f1.series(x,12)

%md 
### <b>Theorem (2nd $q$-analogue):</b> $\frac{rpp_{\lambda/\mu}(q)}{rpp_{\lambda}(q)} = \sum_{T \in SSYT(\mu, d)} q^{p^*(T)} \prod_{u \in [\mu]} (1 - q^{w(u, T(u))})$,
### where $p^*(T) = \sum_{u \in [\mu], T(u) < k \leq \min(d, M_T(u))} w(u, k)$, and $M_T(u)$ is the maximum $k\geq T(u)$ such that replacing $T(u)$ by $k$ still gives a semistandard tableaux.

f2 = x^0*(1-x^2)*(1-x^3) + x^2*(1-x^2)*(1-x^3) + x^3*(1-x^2)*(1-x^3)+x^4*(1-x)*(1-x^3)+x^5*(1-x)*(1-x^3)+x^6*(1-x^1)*(1-x^3)
f2

f2.series(x,12)