# -*- coding: utf-8 -*-
"""RootMultsCode (June 2025, Peter).ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1gBwDu9CSZKgLjYX4BT6aAY-3-j9UmQ-8

# Root Mult Code (June 2025 Peter)
"""

#This is code made to accompany the paper ``Quiver Varieties and Root Multiplicities in Rank 3''
#Code was written by: Patrick Chan and Peter Tingley

#Original file is located at
#    https://colab.research.google.com/drive/1uWtTvfkVePJL6a5zLvhqjl1XBmK1c6k4

from pickle import TRUE
from math import *
class DyckSet():
    '''
    This class of objects encodes a path as in Equation (1.2).
    1^a_1 2^b_1 3^c_1 \cdots  1^a_n 2^b_n 3^c_n is stored as
    [(1,a_1), (2,b_1), (3,c_1), (1,a_2), (2,b_2), (3,c_2),..., (1,a_n), (2,b_n), (3,c_n)]

    Example:
    DyckSet([(1,2),(2,3),(3,1)]) -> [1,1,2,2,2,3] which is equivalent to 1^2 2^3 3^1

    The input to DyckSet comprises of an ordered list of ordered pairs. In explanations,
    an ordered pair (x,y) will be called a uniform section of the DyckSet, x
    the value of a uniform section, and y is the length of the uniform section.

    Remarks:
    - This class recognizes and simplifies inputs, i.e. if two adjacent uniform sections share a common
      value, the class will concatinate the two uniform sections, into one uniform section. e.g.
      DyckSet([(x,a),(x,b)]) -> DyckSet([(x,a+b)]).
      The class does not remove uniform sections of length zero. This will be
      important for indexing.
    - The representation of an object and how it is printed as a string are different.
      When this it is returned, it will be reoresented as a word akin to Example 4.1.
      When printed, it will be represented as a list of ordered pairs.
    '''
    def __init__(self,listoftuples):
        '''
        This simplifies objects in the class by combining adjacent uniform sections that share a common
        value.
        e.g. DyckSet([(1,2),(1,1),(2,1),(2,3),(3,1)]) -> DyckSet([(1,3),(2,4),(3,1)])
        '''
        k=1
        self.values=listoftuples.copy()
        self.values.append((0,0)) #creates buffer for index
        while k+1 <= len(self.values):
            while self.values[k][0] == self.values[k-1][0]:
                self.values[k-1]=(self.values[k-1][0],self.values[k][1]+self.values[k-1][1])
                self.values.pop(k)
            k=k+1
        self.values.pop() #removes buffer

    def __repr__(self):
        '''
        Represents DyckSet([(1,2),(2,3),(3,1)]) as [1,1,2,2,2,3]
        '''
        Output=[]
        for k in range(len(self.values)):
            for l in range(self.values[k][1]):
                Output.append(self.values[k][0])
        return str(Output)

    def __str__(self):
        '''
        When printed, the string of the DyckSet will be returned as a list or ordered pairs.
        Sets the string of DyckSet([(1,2),(2,3),(3,1)]) as [(1,2),(2,3),(3,1)]
        '''
        return str(self.values)

    def add(self,other):
        '''
        Viewing DyckSet as a list of ordered pairs, the add function appendeds the second list to
        the end of the first
        DyckSet([(1,2),(2,3),(3,1])+DyckSet([(3,1),(1,1)]) -> DyckSet([(1,2),(2,3),(3,2),(1,1)])
        '''
        return DyckSet(self.values+other.values)

    def __add__(self,other):
        if isinstance(other,DyckSet):
            return self.add(other)
        else:
            return self.add(DyckSet([other]))

    def __getitem__(self,a):
        return self.values[a]

    def __len__(self):
        return len(self.values)

#################################################################################################
# All functions that check conditions or constuct paths which adhere to a condition in the paper
# must have a way to translate the paper's
# $$\cdots 1^{a_i}2^{b_i}3^{c_i}1^{a_{i+1}}2^{b_{i+1}}3^{c_i} \cdots$$
# to the indexing of these DyckSets in the code.
# All conditions will assume the inputted path satisfies the following:
# (1) The number of ordered pairs is a multiple of three, and these are ordered as
#     (1,a_0), (2,b_0), (3,c_0), (1,a_1), (2,b_1), (3,c_1), \cdots.
# (2) There are no section (1,a_i), (2,b_i), (3,c_i) where b_i and c_i are both zero.
#     This allows use to index in multiples of threes. i.e. for k in range(0,len(path),3),
#     path[k] is the ordered pair (1,a_{k/3}), path[k+1] is (2,b_{k/3}),
#      and path[k+2] is (3,c_{k/3}).
#################################################################################################

def PathCandidates(ones,twos,threes):
    '''
    Constructs all possible paths that follow Condition 1 of Thm 1.1, but with
    the annoying flaw that the final bk and ck may both be 0 - that is, the path
    may end with a 1. Those cases are eliminated in teh ListofDyckPaths
    function.
    '''
    answer=[]
    if twos or threes >=1:
        if ones >= 1:
            for i in range(twos+1):
                for j in range(threes+1):
                    partialperms=PathCandidates(ones-1,twos-i,threes-j)
                    for partial in partialperms:
                        partial=partial+DyckSet([(1,1)]) #Ensures a_i \neq 0
                        if i !=0 or j!=0:
                            partial=partial+DyckSet([(2,i)])+DyckSet([(3,j)])
                            answer.append(partial)
                        else:
                            answer.append(partial)
            return answer
        else:
            return []
    else:
        return [DyckSet([(1,ones)])]

def ListofDyckPaths(ones,twos,threes,s,t):
    '''
    This function returns a list of Dyck Paths given the parameter input.
    This is used in Root_Mult_Estimates().
    '''
    DyckPaths=[]
    Candidates=PathCandidates(ones,twos,threes)
    Candidates=[path for path in Candidates if path[-1][0]!=1]
    for c in Candidates:
        if Condition2(c,s,t,ones,twos,threes):
            if Condition3(c,s,t,ones,twos,threes):
                DyckPaths.append(c)
    return DyckPaths

def List_of_Valid_Paths(ones,twos,threes,s,t):
    '''
    This function is meant for use in small cases, returning a list of all paths
    for the user to comb through. Nothing is contingent on this function.
    '''
    ValidPaths=[]
    Candidates=ListofDyckPaths(ones,twos,threes,s,t)
    for c in Candidates:
        if Condition4(c,s,t,ones,twos,threes):
            if Condition5(c,s,t,ones,twos,threes):
                if Condition6(c,s,t,ones,twos,threes):
                    if RefinedConditions(c,s,t,ones,twos,threes):
                        ValidPaths.append(c)
    return ValidPaths

def TestPath_st(path,s,t):
    A=0
    B=0
    C=0
    for k in range(0,len(path),3):
        A=A+path[k][1]
        B=B+path[k+1][1]
        C=C+path[k+2][1]
    '''
    This function tests a path under the Dyckpath class under all conditions
    (excluding Thm 1.1 Condition 1, which our notations assumes holds),
    printing all conditions broken. This is not used in finding estimates, it
    is here to allow experimentation.
    '''
    passed=True
    if Condition2(path,s,t,A,B,C) == False:
        print('Fails Thm 1.1, Condition 2')
        passed=False
    if Condition3(path,s,t,A,B,C) == False:
        print('Fails Thm 1.1, Condition 3')
        passed=False
    if Condition4(path,s,t,A,B,C) == False:
        print('Fails Thm 1.1, Condition 4')
        passed=False
    if Condition5(path,s,t,A,B,C) == False:
        print('Fails Thm 1.1, Condition 5')
        passed=False
    if Condition6(path,s,t,A,B,C) == False:
        print('Fails Thm 1.1, Condition 6')
        passed=False
    if RefinedConditions(path,s,t,A,B,C) == False:
        print('Fails Thm 5.1')
        passed=False
    if passed==True:
        print('Passes all conditions')

#################################################################################################
# The following condition functions all contain A, B, and C as input variables. This is
# to allow the code to check if subpaths follow the conditions.
#################################################################################################

def Condition2(path,s,t,A,B,C):
    '''
    This function is used to check Thm 1.1: Condition 2: "Rational Dyck Path"/Slope Test.
    Used in List_of_Valid_Paths(), ListofDyckPaths(), TestPath_st()
    '''
    a=0 #a will act as the summation of a_i for i=1,...,k
    b=0 #b will act as the summation of b_i for i=1,...,k
    c=0 #c will act as the summation of c_i for i=1,...,k
    for k in range(0,len(path),3):
        a=a+path[k][1]
        b=b+path[k+1][1]
        c=c+path[k+2][1]
        if (a/((s*b)+(t*c))) < A/((s*B)+(t*C)):
            return False
    return True

def Condition3(path,s,t,A,B,C):
    '''
    This function is used to check Thm 1.1: Condition 3.
    Used in List_of_Valid_Paths(), ListofDyckPaths(), TestPath_st()
    '''
    a=0
    b=0 #b will act as the summation of b_i for i=1,...,k
    c=0 #c will act as the summation of c_i for i=1,...,k
    for k in range(0,len(path),3):
        a=a+path[k][1]
        b=b+path[k+1][1]
        c=c+path[k+2][1]
        if c>0:
            if (a/((s*b)+(t*c)))==(A/((s*B)+(t*C))):
                if b/c < B/C:
                    return False
    return True

def Condition4(path,s,t,A,B,C):
    '''
    This function is used to check Thm 1.1: Condition 4 (both the s and t cases).
    Used in List_of_Valid_Paths(), Root_Mult_Estimates(), TestPath_st()
    '''
    for k in range(0,len(path),3):
          if path[k+1][1] > s*(path[k][1]): #Thm 1.1: Condition 4s: b_i < s*a_i
              return False
          if path[k+2][1] > t*(path[k][1]): #Thm 1.1: Condition 4t: c_i < t*a_i
              return False
    return True

def Condition5(path,s,t,A,B,C):
    '''
    This function is used to check Thm 1.1: Condition 5.
    Used in List_of_Valid_Paths(), Root_Mult_Estimates(), TestPath_st()
    '''
    if len(path)>3:
        for k in range(0,len(path)-4,3):
            nb = min(path[k+1][1],s*(path[k+3][1])-path[k+4][1])
            nc = min(path[k+2][1],t*(path[k+3][1])-path[k+5][1])
            if path[k+3][1] > (s*nb)+(t*nc)-max((1/s)*nb,(1/t)*nc):
                return False
    return True

def Condition6(path,s,t,A,B,C):
    '''
    This function is used to check Thm 1.1: Condition 6.
    This function is used in List_of_Valid_Paths(), Root_Mult_Estimates(), TestPath_st()
    '''
    if len(path)>3:
        for k in range(0,len(path)-4,3):
            nb = min(path[k+1][1],s*(path[k+3][1])-path[k+4][1])
            nc = min(path[k+2][1],t*(path[k+3][1])-path[k+5][1])
            if ((t*nc)+(s*nb)) !=0 :
                if (path[k+3][1]/((t*nc)+(s*nb)) > (1/2+sqrt((((s**2)+(t**2))**2)-4*((s**2)+(t**2)))/(2*((s**2)+(t**2))))):
                    return False
            else:
                print("Thm 1.1, condition 6 results in divsions by 0. Condition 5 should have failed")
                return None
    return True

def RefinedConditions(path,s,t,A,B,C):
    '''
    This function tests both conditions laid out in Thm 5.1.
    Used in List_of_Valid_Paths(), Root_Mult_Estimates(), and TestPath_st().
    '''
    ones=A
    twos=B
    threes=C
    if len(path) > 3:
        for j in range(0,len(path)-3,3):
            Aj=0  #records a_1 + \cdots + a_j, the number of 1s up to j
            Bj=0  #records b_1 + \cdots + b_j, the number of 2s up to j
            Cj=0  #records c_1 + \cdots + c_j, the number of 3s up to j
            Aij=0 #records a_i + \cdots + a_j
            Bij=0 #records b_i + \cdots + b_j
            Cij=0 #records c_i + \cdots + c_j
            #This loop calculates Bj,Bj,Cj for a fixed j.
            for k in range(0,j+1,3):
                Aj=Aj+path[k][1]
                Bj=Bj+path[k+1][1]
                Cj=Cj+path[k+2][1]
            #This loop cycles over all i<j and check the condition for that pair i,j.
            for i in range(j,-1,-3):
                Aij=Aij+path[i][1]
                Bij=Bij+path[i+1][1]
                Cij=Cij+path[i+2][1]

                #These are $n_{A_{ij}}$ from the paper, and similarly for B and C.
                nAij=Aij+path[j+3][1]-path[i][1]
                nBij=min(Bij,s*nAij-(Bij-path[i+1][1]+path[j+4][1]))
                nCij=min(Cij,t*nAij-(Cij-path[i+2][1]+path[j+5][1]))

                #These are $\tilde a_{ij}$, and similarly for b and c.
                ta=Aj-Aij+s*nBij+t*nCij-nAij
                tb=Bj-Bij+nBij
                tc=Cj-Cij+nCij

                if s*tb+t*tc == 0:
                    print("Thm 5.1 is not applicable due to a failure of previous conditions")
                    return None
                if (ta/((s*tb)+(t*tc)))<(ones/((s*twos)+(t*threes))): #Tests if path satisfies Thm 5.1 (1)
                    return False
                elif (ta/((s*tb)+(t*tc)))==(ones/((s*twos)+(t*threes))) and tc > 0: #Tests if path satisfies Thm 5.1 (2)
                    if (tb/tc) < (twos/threes):
                        return False
    return True

#################################################################################################
# This marks the end of the conditions, the following is the main functions of
# this code, creating our estimates.
#################################################################################################

def Root_Mult_Estimates(ones,twos,threes,s,t):
    '''
    This prints the number of Dyck Paths, and the estimate using Thm 1.1, and the estimate
    applying the refined condition.
    '''
    DyckPaths=ListofDyckPaths(ones,twos,threes,s,t)
    Thm1paths=[]
    Thm2paths=[]
    print("The number of Dyck paths is " + str(len(DyckPaths)))
    for p in DyckPaths:
        if Condition4(p,s,t,ones,twos,threes):
            if Condition5(p,s,t,ones,twos,threes):
                if Condition6(p,s,t,ones,twos,threes):
                    Thm1paths.append(p)
    print("The estimate using Theorem 1.1 is " + str(len(Thm1paths)))
    for p in Thm1paths:
        if RefinedConditions(p,s,t,ones,twos,threes):
            Thm2paths.append(p)
    print("The estimate using Theorem 5.1 is " + str(len(Thm2paths)))

#This is Example 4.1
List_of_Valid_Paths(4,3,2,2,1)

#This is the path from Example 4.2.
path=DyckSet([(1,3),(2,3),(3,0),(1,3),(2,2),(3,1),(1,4),(2,5),(3,0)])
TestPath_st(path,2,2)

#This is the path from Example 4.4.
path=DyckSet([(1,4),(2,3),(3,0),(1,3),(2,2),(3,1),(1,4),(2,2),(3,1),(1,4),(2,6),(3,0),(1,6),(2,8),(3,0)])
TestPath_st(path,2,2)

#This is the path from Example 4.4.
path=DyckSet([(1,4),(2,4),(3,0),(1,4),(2,3),(3,1),(1,4),(2,6),(3,0),(1,5),(2,6),(3,0)])
TestPath_st(path,2,2)

#This is the path from Example 6.1.
path=DyckSet([(1,3),(2,2),(3,0),(1,2),(2,1),(3,1),(1,1),(2,2),(3,0)])
TestPath_st(path,2,2)

#This is the one of the paths from Example 6.2.
path=DyckSet([(1,2),(2,2),(3,0),(1,2),(2,1),(3,1),(1,1),(2,1),(3,1),(1,1),(2,1),(3,1)])
TestPath_st(path,2,1)

#This is the one of the minimal case, from Example 6.6, of a path that passes all
#conditions but does not correspond to a stable component.
path=DyckSet([(1,2),(2,1),(3,1),(1,1),(2,1),(3,1),(1,1),(2,1),(3,1)])
TestPath_st(path,2,1)

#This is the path from Example 7.3.
path=DyckSet([(1,30),(2,20),(3,10),(1,15),(2,20),(3,6),(1,33),(2,40),(3,20)])
TestPath_st(path,2,1)

#Finds estiates for the root (8,7,4) when s=2,t=1
Root_Mult_Estimates(8,7,4,2,1)

#this generates a table with the root mutiplicities and estimates for many cases
#at once, with s=t=2. Since the actual multiplicities are calculated elsewhere,
# these must be entered by hand.
# Fr the paper, we made a bigger table, but that code does not run on colab.

def ExampleTable_2_2():
  Cases=[(1,1,1,1),(2,1,1,2),(2,2,1,3),(2,1,2,3),(3,2,1,4),(3,1,2,4),(3,3,1,5),(3,1,3,5),(4,3,1,7),(4,1,3,7),(4,4,1,10),(4,1,4,10),(3,2,2,10),(3,3,2,12),(3,2,3,12),(5,4,1,13),(5,1,4,13),(5,5,1,16),(5,1,5,16),(6,5,1,21),(6,1,5,21),(6,6,1,28),(6,1,6,28),(4,3,2,25),(4,2,3,25),(7,6,1,35),(7,1,6,35),(7,7,1,43),(7,1,7,43),(8,7,1,55),(8,1,7,55),(4,3,3,46),(4,4,3,58),(4,3,4,58),(5,4,2,61),(5,2,4,61),(8,8,1,70),(8,1,8,70)]
  print("Case:     Actual:    Refined Est:     Initial Est:    Total Dyck Paths:")
  for quad in Cases:
    DyckPaths=ListofDyckPaths(quad[0],quad[1],quad[2],2,2)
    Thm1paths=[]
    Thm2paths=[]
    for p in DyckPaths:
        if Condition4(p,2,2,quad[0],quad[1],quad[2]):
            if Condition5(p,2,2,quad[0],quad[1],quad[2]):
                if Condition6(p,2,2,quad[0],quad[1],quad[2]):
                    Thm1paths.append(p)
    for p in Thm1paths:
        if RefinedConditions(p,2,2,quad[0],quad[1],quad[2]):
            Thm2paths.append(p)
    print("("+str(quad[0])+","+str(quad[1])+","+str(quad[2])+")   "+str(quad[3])+"          "+str(len(Thm2paths))+"                "+str(len(Thm1paths))+"              "+str(len(DyckPaths)))

ExampleTable_2_2()

#This generates a table with the root mutiplicities and estimates for many cases
#at once, with s=2, t=1. Again, actual multiplicities must be entered by hand.
#Caution: this code takes some time to run, and crashes on colab
def ExampleTable_2_1():
  Cases=[(2,2,1,2),(3,3,1,3),(4,3,2,5),(4,4,1,5),(5,5,1,7),(5,4,2,11),(6,6,1,11),(5,5,2,15),(7,7,1,15),(6,5,2,22),(8,8,1,22),(6,5,3,30),(9,9,1,30),(7,6,2,42),(10,10,1,42),(7,7,2,56),(11,11,1,56),(7,6,3,77),(8,7,2,77),(12,12,1,77),(7,7,3,101),(9,8,2,135),(8,7,3,176),(9,9,2,176),(8,7,4,231),(8,8,3,231),(10,9,2,231),(9,7,4,297),(9,8,3,385),(11,10,2,385),(9,8,3,385),(11,11,2,490),(9,8,4,627),(9,9,4,792),(10,9,3,792),(10,8,5,1002)]
  print("Case:     Actual:    Refined Est:     Initial Est:    Total Dyck Paths:")
  for quad in Cases:
    DyckPaths=ListofDyckPaths(quad[0],quad[1],quad[2],2,1)
    Thm1paths=[]
    Thm2paths=[]
    for p in DyckPaths:
        if Condition4(p,2,1,quad[0],quad[1],quad[2]):
            if Condition5(p,2,1,quad[0],quad[1],quad[2]):
                if Condition6(p,2,1,quad[0],quad[1],quad[2]):
                    Thm1paths.append(p)
    for p in Thm1paths:
        if RefinedConditions(p,2,1,quad[0],quad[1],quad[2]):
            Thm2paths.append(p)
    print("("+str(quad[0])+","+str(quad[1])+","+str(quad[2])+")   "+str(quad[3])+"          "+str(len(Thm2paths))+"                "+str(len(Thm1paths))+"              "+str(len(DyckPaths)))