##################################################################
# Code relating to data presented in paper
# "Quasisymmetric and Schur expansions of cycle index polynomials"
# by Nick Loehr and Greg Warrington
##################################################################
# First few sections define needed functions. Go to end for usage.
# Code has not been cleaned up much or made user friendly, but
# there are still a few user-friendly functions (see end). You may
# need to change the 'homepath' string below for some functions to
# work. Email me if you have questions!
##################################################################
import sys
from sage.combinat.subset import Subsets
Sym = SymmetricFunctions(QQ)
s = Sym.schur()
m = Sym.monomial()
QSym = QuasiSymmetricFunctions(QQ)
F = QSym.F()
M = QSym.M()
# Insert path to writeable directory here
homepath = '/home/gswarrin/research/polya/'
#######################################################################
# Code for performing basic manipulations
#######################################################################
#######################################################################
def cmp_alpha_to_sub(alpha):
    """ 
    Convert composition as a word (i.e., sequence of partial sums) to subset
    For example, [3,5,9] -> {3,8}
                 [3,5,6] -> {3,8}
                 [1,3,4,5,9] -> {1,4,8,13}
    """ 
    if len(alpha) == 1:
        return "{}"
    ans = [alpha[0]]
    for i in range(1,len(alpha)-1):
        ans.append(ans[i-1]+alpha[i])
    return "{" + ','.join(map(lambda x: repr(x), ans)) + "}"

##########################################################################
def symm_to_F_list(f,n):
    """ convert a symmetric function to its F-expansion as list of coefficients
    compositions are encoded as words (alpha) in order returned by Compositions(n)
    """
    # make sure it's an F-expansion
    f = F(f)       
    fterms = str(f).split('+')
    # convert compositions to list of strings so we can 
    # recognize them from F-expansion of f
    cmps = map(lambda x: str(x), Compositions(n).list())
    ans = [0 for i in cmps]
    # run through the terms in F-expansion, extracting coefficient  and partition
    for tm in fterms:
        if '*' in tm:
            a,b = tm.split('*')
        else:
            a = 1
            b = tm
        tm_coeff = int(a)
        c,d = b.split('F')
        tm_ptn = d
        if tm_ptn[-1] == ' ':
            tm_ptn = tm_ptn[:-1]
        # Find index of composition from F-expansion in list of compositions
        idx = cmps.index(tm_ptn)
        # increment coefficient accordingly
        ans[idx] = ans[idx] + tm_coeff
    return ans

##########################################################################
# Code for finding algebraic expansions
##########################################################################
def get_des_ides_mat(k,N):
    """ Make an array indexed by compositions of size N and size k
        Keep track of how many have a given descent set, inverse descent set pair
    """
    Ncmps = Compositions(N).list()
    kcmps = Compositions(k).list()
    cnt = 0
    arr = [[0 for idx in range(pow(2,N-1))] for jdx in range(pow(2,k-1))]
    for x in SymmetricGroup(N).list():
        px = Permutation(x)
        # figure out actual descent set of permutation (intersected with [k-1])
        cur_desc = filter(lambda y: px[y] > px[y+1], range(k-1))
        cur_desc.append(k-1)
        cmp = [cur_desc[i]-cur_desc[i-1] for i in range(1,len(cur_desc))]
        cmp = [cur_desc[0]+1] + cmp
        # figure out corresponding inverse descent set
        ides = filter(lambda y: px.index(y+1) < px.index(y), range(1,len(px)))
        ides.append(N)
        icmp = [ides[i] - ides[i-1] for i in range(1,len(ides))]
        icmp = [ides[0]] + icmp
        arr[kcmps.index(cmp)][Ncmps.index(icmp)] += 1
    # print matrix(QQ,int(pow(2,k-1)),int(pow(2,N-1)),arr)
    return arr

##########################################################################
def find_alg_expn(f,G,k,N,verbose=False,sset_only=True):
    """ View G <= S_k action on words of length N
        Try to write as sum over FQ indexed by permutations with legal descent set in some subset of [k-1]
    """
    # get F-expansion of cycle index of cyc(G)
    subG = SymmetricGroup(N).subgroup(G)
    cycG = subG.cycle_index()
    Fcoeffs = symm_to_F_list(cycG,N)
    if verbose:
        print Compositions(N)
        print Fcoeffs
    # keeps track of how each subset of [k-1] affects F-expansion
    Texpans = get_des_ides_mat(k,N)
    T = []
    
    def _rec_search_alg_T(curFcoeffs,Texpans,curLC,idx):
        """ curFcoeffs is difference between Fcoeffs and what we've already selected
            Texpans gives matrix encoding how each t in T modifies F-expansion
            curLC is the current linear combination we're using
            idx is the index we're working on
        """
        if sum(curFcoeffs) == 0:
            Gs = subG.gens_small()
            Grep = str(repr(subG)).replace("(Symmetric group of order ","S").replace('! as a permutation group)','')
            if max(curLC) > 1:
                set_type = "MSet"
            else:
                set_type = "SSet"
            print "          +A+ Algebraic proof found ",
            cmps = Compositions(k).list()
            tmpstr = '{'
            for i in range(len(cmps)):
                if curLC[i] > 0:
                    tmpstr += "%s" % (cmp_alpha_to_sub(cmps[i]))
                    print "%d*%s" % (curLC[i],cmp_alpha_to_sub(cmps[i])),
                    if i < len(cmps)-1:
                        print "+",
                        tmpstr += ","
            tmpstr = tmpstr.replace('{}','\\emptyset').\
                        replace('{','\\{').replace('}','\\}')
            print
            f.write("%s\\}\\\\\n" % tmpstr)
            return True
            
        for i in range(idx,len(Texpans)):
            modFcoeffs = [curFcoeffs[idx] for idx in range(len(curFcoeffs))]
            isokay = True
            for j in range(len(modFcoeffs)):
                modFcoeffs[j] -= Texpans[i][j]
                if modFcoeffs[j] < 0:
                    isokay = False
            if isokay == True:
                curLC[i] += 1
                if sset_only == True:
                    if _rec_search_alg_T(modFcoeffs,Texpans,curLC,i+1) == True:
                        return True
                else:
                    if _rec_search_alg_T(modFcoeffs,Texpans,curLC,i) == True:
                        return True
                curLC[i] -= 1
        return False
    return _rec_search_alg_T(Fcoeffs,Texpans,[0 for i in range(len(Texpans))],0)

################################################################################
def find_alg_expansions(k,N,sset_only,pr_to_file=True):
    """ 
    find all algebraic expansions obtained by summing over permutations 
    with permissible descent sets (multisets possibly allowed)
    """
    curSk = SymmetricGroup(k)
    fstr = homepath + ('alg-expn-k%d-N%d' % (k,N))
    if pr_to_file:
        if sset_only == True:
            f = open(fstr + '-ssetonly.txt','w')
        else:
            f = open(fstr + '-mset.txt','w')
    # run through representative subgroup from each conjugacy class
    for G in curSk.conjugacy_classes_subgroups():
        if pr_to_file:
            f.write('\n')
            f.write("G is <%s> <= S%d; |G|=%d\n" % \
            (repr(G.gens_small()),N,G.order()))
        else:
            print 
            print "G is <%s> <= S%d; |G|=%d" % (repr(G.gens_small()),N,G.order())
        # second-to-last argument is simply in regards to how verbose
        # first argument is vestigal - ignore it.
        if not find_alg_expn('',G,k,N,False,sset_only):
            print "          -A- None found"
    if pr_to_file:
        f.close()
    
    
###################################################################################
# Code for working with bijective expansions
###################################################################################
###################################################################################
# code for finding all words that standardize to a given word
###################################################################################
def cmp_to_zwd(cmp):
    """ take word as alpha and convert to a zwd
        eg 131 -> 12223
    """
    ans = []
    for i in range(1,len(cmp)+1):
        ans.extend([i for idx in range(cmp[i-1])])
    return ans
    

def get_legal_zwds(px,cmps,cmpsums,flatwds):
    """ return z with ascents at ides of px
    """
    ides = filter(lambda y: px.index(y+1) < px.index(y), range(1,len(px)))
    # print "ides: ",ides
    ans = []
    for i in range(len(cmpsums)):
        isokay = True
        for j in range(len(ides)):
            if ides[j] not in cmpsums[i]:
                isokay = False
                
        if isokay == True:
            ans.append([flatwds[i][px[idx]-1] for idx in range(len(px))])
    return ans

###################################################################################
# code for computing *all* (even flattened) orbits
###################################################################################   
def find_cmp_orbit(G,perms,ides,c,csum,fwd):
    """ find orbits corresponding to a given composition
        perms are all permutations of length n
        ides are the corresponding inverse descent sets (as subsets)
        c is a composition (as an alpha)
        csum is partial sums of c
        fwd is the flatwd corres to c
    """
    # first find all zwds with given content
    allwds = []
    for i in range(len(perms)):
        isokay = True
        for j in range(len(ides[i])):
            if ides[i][j] not in csum:
                isokay = False
                continue
                
        if isokay == False:
            continue
        allwds.append([fwd[perms[i][idx]-1] for idx in range(len(perms[i]))])
    
    # now split those zwds up into orbits
    orbits = [] 
    while len(allwds) > 0:
        px = allwds.pop(0)
        tmp_orbit = [px]
        # run through group elements
        for g in G:
            pg = list(Permutation(g))
            n = len(pg)
            y = [px[pg.index(i)] for i in range(1,len(px)+1)]
            if y not in tmp_orbit:
                tmp_orbit.append(y)
                allwds.pop(allwds.index(y))
        orbits.append(tmp_orbit)
    return orbits

def get_all_orbits(G,N,perms,ides,cmps,cmpsums,flatwds):
    """
    """
    orbits = []
    stidices = [0]
    stposn = 0
    for i in range(len(cmps)):
        corbit = find_cmp_orbit(G,perms,ides,cmps[i],cmpsums[i],flatwds[i])
        stposn += len(corbit)
        stidices.append(stposn)
        orbits.extend(corbit)
    return orbits,stidices
    

###########################################################################################
# Code for associating LC of orbits to each permutation according to ones it hits
###########################################################################################
def get_orbitLC(perms,G,N,myverbose=False):
    """
    """
    ides = [filter(lambda y: perms[i].index(y+1) < perms[i].index(y), range(1,N))\
                 for i in range(len(perms))]
    cmps = Compositions(N).list()
    cmpsums = [[sum(cmps[i][:j]) for j in range(1,len(cmps[i]))] \
                for i in range(len(cmps))]
    flatwds = [cmp_to_zwd(cmp) for cmp in cmps]
    orbits,stdices = get_all_orbits(G,N,perms,ides,cmps,cmpsums,flatwds)
    if myverbose: 
        for o in orbits:
            tt = []
            for elt in o:
                tmp = [0 for j in range(len(elt))]
                for j in range(len(elt)):
                    if elt[j] > elt[(j+1)%len(elt)]:
                        tmp[j] = 1
                tt.append(tmp)
            print "Orbit: ",o,tt
    LC = []
    for px in perms:
        zwds = get_legal_zwds(px,cmps,cmpsums,flatwds)
        curLC = [0 for i in range(len(orbits))]
        for z in zwds:
            zsort = [len(filter(lambda y: y == i, sorted(z))) for i in range(1,max(z)+1)]
            tmppos = cmps.index(zsort)
            stpos = stdices[tmppos]
            for j in range(stpos,stdices[tmppos+1]):
                if z in orbits[j]:
                    curLC[j] += 1
        LC.append(curLC)
    return LC

########################################################################################################
def reorder_perms_LC(perms,LC):
    """
    """
    seen = []
    nperms = []
    nLC = []
    fnonz = []
    for i in range(len(LC[0])):
        # which LC are non-zero for this orbit
        ilist = filter(lambda y: LC[y][i] > 0 and y not in seen, range(len(LC)))
        nperms.extend([perms[y] for y in ilist])
        nLC.extend([LC[y] for y in ilist])
        fnonz.extend([i for idx in range(len(ilist))])
        seen.extend(ilist)
    return nperms,nLC,fnonz
    

########################################################################################################   
def find_bij_arb_G(f,origG,N):
    """ Look for arbitrary LC of permutations that yields cyc(G)
        These are bijective solutions using arbitrary sets of permutations (not nec defined by descents)
        Works for n=4, but gets bogged down quickly for n=5.
    """
    G = SymmetricGroup(N).subgroup(origG)
    operms = Permutations(N)
    oLC = get_orbitLC(operms,G,N)
    perms,LC,fnonz = reorder_perms_LC(operms,oLC)
    
    nfound = 0
    
    #########################
    def _rec_search_bij(perms,LC,fnonz,curCoeffs,curLC,idx,nleft):
        """ curCoeffs is difference between vector of all 1s and what we've already selected
            pxLC gives matrix encoding which orbits each permutation px hits
            curLC is the current linear combination we're using
            idx is the index we're working on
            npos is the number we've turned on
        """
        if sum(curCoeffs) == 0:
            tans = [perms[i] for i in filter(lambda y: curLC[y] == 1, range(len(curLC)))]
            print "      +S+ bijective FQ expn via IDes and stdn on subset S: ",tans
            tans = str(tans).replace(', ','').replace('[[','\\{').\
                   replace(']]','\\}').replace('[','').replace(']',',')
            f.write("%s & " % (tans))
            sys.stdout.flush()
            return True
        
        if nleft <= 0 or len(LC)-idx < nleft:
            return False
               
        # find position of first nonzero in curCoeffs
        cidx = 0
        while curCoeffs[cidx] == 0:
            cidx += 1
        for i in range(idx,len(LC)):
            if fnonz[i] > cidx:
                return False
            if fnonz[i] < cidx:
                continue
            modCoeffs = [curCoeffs[jdx] for jdx in range(len(curCoeffs))]
            isokay = True
            for j in range(len(modCoeffs)):
                modCoeffs[j] -= LC[i][j]
                if modCoeffs[j] < 0:
                    isokay = False
            if isokay == True:
                curLC[i] += 1
                if _rec_search_bij(perms,LC,fnonz,modCoeffs,curLC,i+1,nleft-1) == True:
                    return True
                curLC[i] -= 1
        return False
    return _rec_search_bij(perms,LC,fnonz,[1 for i in range(len(LC[0]))],\
                           [0 for i in range(len(LC))],0,factorial(N)/G.order())

def get_descent_subsets(N):
    """ return all 2^(2^{N-1}) sets of possible descent sets
    """
    subs = []
    for i in range(N):
        tmp = map(lambda y: list(y), Subsets(range(N-1),i))
        subs.extend(map(lambda x: list(x), list(tmp)))
    
    # Since our set consists of lists, we can't use Subsets for it since not hashable
    # So we'll do recursively
    def _rec_get_powerset(S,curSet,idx,ans):
        """
        """
        if idx == len(S):
            ans.append(curSet)
            return ans
            # print "One sub: ",curSet
            # return
            
        ans = _rec_get_powerset(S,curSet,idx+1,ans)
        nSet = [[y for y in x] for x in curSet]
        nSet.append(S[idx])
        ans = _rec_get_powerset(S,nSet,idx+1,ans)
        return ans
    
    ans = _rec_get_powerset(subs,[],0,[])
    return ans

#####################################################################################################
def defined_by_des(T,perms):
    """ Check if set of permutations is defined according to a particular collection of descent sets
    """
    N = len(perms[0])
    Tdes = []
    for px in T:
        tmp = filter(lambda y: px[y] > px[y+1], range(N-1))
        if tmp not in Tdes:
            Tdes.append(tmp)
    for px in Permutations(N):
        if px not in T and filter(lambda y: px[y] > px[y+1], range(N-1)) in Tdes:
            return False,None
    return True,Tdes
        

####################################################################################################
def get_perms_from_des(k,N,perms,legal_des):
    """ return permutations with descent set cap [k-1] lying in legal_des
    """
    ans = []
    for i in range(len(perms)):
        tmp = filter(lambda y: perms[i][y] > perms[i][y+1], range(k-1)) 
        if tmp in legal_des:
            ans.append(i)
            # print "Appending: ",perms[i]," for ",legal_des
    return ans

####################################################################################################
def find_bij_des_G(f,G,k,N,myverbose=False):
    """ similar to find_LC_that_work
        looks for collections of permutations *defined by descent sets* that yield cyc(G)
    """
    perms = Permutations(N)
    LC = get_orbitLC(perms,G,N,myverbose)
    des_subs = get_descent_subsets(k)
    cnt = 0
    for legal_des in des_subs:
        cnt += 1
        des_perms_indices = get_perms_from_des(k,N,perms,legal_des)
        coeffs = [1 for i in range(len(LC[0]))]
        isokay = True
        for i in des_perms_indices:
            for j in range(len(LC[i])):
                coeffs[j] -= LC[i][j]
                if coeffs[j] < 0:
                    isokay = False
                    continue
            if isokay == False:
                continue
        if isokay == True and max(coeffs) == 0:
            tmpstr = str(map(lambda x: map(lambda y: y+1, x), legal_des))
            tmpstr = tmpstr.replace('[]','\\emptyset').replace('[','\\{').replace(']','\\}')
            print "  +T+ bijective FQ expn via IDes and stdn with descent set T: %s" % (tmpstr)
            f.write("%s & " % (tmpstr))
            return True
    return False
            

####################################################################################################
def vals_to_subs(vals,N):
    """ convert numbers to subsets via binary encoding
    """
    ans = []
    for v in vals:
        binexp = list(Integer(v).binary())
        tmp = []
        for idx in range(len(binexp)):
            if binexp[idx] == '1':
                tmp.append(idx+(N-len(binexp)-1))
        ans.append(tmp)
    return ans

####################################################################################################
def vals_to_csubs(vals,N):
    """ convert numbers to subsets via binary encoding
    """
    ans = []
    for v in vals:
        binexp = list(Integer(v).binary())
        tmp = []
        for idx in range(len(binexp)):
            if binexp[idx] == '1':
                tmp.append(idx+(N-len(binexp)))
        ans.append(tmp)
        # print v,binexp,tmp
    return ans
    

####################################################################################################
def guess_des(G,k,N,vals,myverbose=False):
    """ similar to find_LC_that_work
        looks for collections of permutations *defined by descent sets* that yield cyc(G)
    """
    perms = Permutations(N)
    LC = get_orbitLC(perms,G,N,myverbose)
    des_subs = [vals_to_subs(vals,N)]
    # print "Done with get_orbitLC",G.order()
    cnt = 0
    # des_subs = [[[],[4,5],[3,5],[2,5],[2,3,4,5],[1,3,4,5],[1,2,4,5],[3,4],[2,4],[0,1,2,3,4,5]]]
    # des_subs = [[[],[5,6],[4,6],[3,6],[3,4,5,6],[2,4,5,6],[2,3,5,6],[4,5],[3,5],[1,2,3,4,5,6]]]
    # des_subs = [[[],[0,1,2,3],[1,3],[2,3]]]
    for legal_des in des_subs:
        cnt += 1
        print legal_des
        des_perms_indices = get_perms_from_des(k,N,perms,legal_des)
        # print legal_des,"dpi: ",des_perms_indices
        coeffs = [1 for i in range(len(LC[0]))]
        isokay = True
        for i in des_perms_indices:
            for j in range(len(LC[i])):
                coeffs[j] -= LC[i][j]
                if coeffs[j] < 0:
                    # print "Done Failed on: ",j
                    isokay = False
                    continue
            # print des_perms_indices,"After i:",i,coeffs
            if isokay == False:
                continue
        if isokay == True and max(coeffs) == 0:
            print "  +T+ bijective FQ expn via IDes and stdn with k=%d, descent set T: %s" % (k,map(lambda x: map(lambda y: y+1, x), legal_des))
            # print "Yes!: ",G,map(lambda x: map(lambda y: y+1, x), legal_des)
            return True
        if cnt % 8000 == 0:
            print "Done with: ",cnt," out of",len(des_subs)
    return False

####################################################################################################
# ***
def find_bij_Gcc(f,origG,N):
    """ Call find_bij_arb on all elements of G's conjugacy class (matters for bijective proof)
    """
    curSN = SymmetricGroup(N)
    G = curSN.subgroup(origG)
    coset_reps = [x[0] for x in curSN.cosets(G, side='left')]
    seen = []
    for x in Permutations(N):
        G = curSN.subgroup(origG)
        Gconj = curSN.subgroup(G.conjugate(x))
        if (N >= 5 and seen == []) or (N <= 4 and Gconj not in seen):
            seen.append(Gconj)
            # print G,x,Gconj
            print "\nLooking for proof for G gen'd by %s in S%d; |G|=%d" % (Gconj.gens_small(),N,G.order())
            # for table for paper
            tmpstr = "%s & %d & " % (Gconj.gens_small(),G.order())
            # print "tmpstr is: ",tmpstr
            tmpstr = tmpstr.replace('[','\\langle').replace(']','\\rangle').replace(',','').replace(') (','),(')
            f.write(tmpstr)
            # print "tmpstr is now: ",tmpstr
            # print map(lambda x: Permutation(x), Gconj.list())
            isdone = False
            if find_bij_des_G(f,Gconj,N,N) == False:
                print "  -T- No subset T found" #  for k = %d" % (k)
                f.write(" & ")
            else:
                f.write("\\\\\n")
                isdone = True
            if (N <= 4 and isdone == False):
                if find_bij_arb_G(f,Gconj,N) == False:
                    print "      -S- No subset S found"
                    f.write(" & ")
                else:
                    f.write("\\\\\n")
                    isdone = True
            if isdone == False: 
                if find_alg_expn(f,Gconj,N,N,False,True) == False:
                    print "          -A- No algebraic proof found"
                    f.write("\\\\\n")
                else:
                    isdone = True

####################################################################################################
def get_maj_perm_indices(k,N,perms):
    """ check that permutations with even maj works for (1,n)(2,n-1)...
    """
    ans = []
    for i in range(len(perms)):
        px = Permutation(perms[i])
        maj = 0
        for j in range(k-1):
            if px[j] > px[j+1]:
                maj += (j+1)
        if maj%2 == 0:
            print "Okay: ",px
            ans.append(i)
    return ans

#####################################################################################################
def find_bij_alg(fn,N):
    """ Run through subgroups looking for a bijective proof
    """
    f = open('/home/gswarrin/research/polya/' + fn,'w')
    for G in SymmetricGroup(N).conjugacy_classes_subgroups():
        print "---Begin--- groups in a given conjugacy class ------------------------------------------"
        find_bij_Gcc(f,G,N)
        f.write('\\midrule')
        print "\n----End---- groups in a given conjugacy class ----------------------------------------"
    f.close()
######################################################################
# Code for finding perfect matchings
######################################################################
def comp_bij(n):
    """
    """
    orbs = []
    seen = []
    cmps = Compositions(n)
    for c in cmps:
        if c not in seen:
            for i in range(len(c)):
                seen.append(c[i:] + c[:i])
            orbs.append(c)
    # print "0 []"
    # for k in range(1,n):
    #     print k,
    #     for y in filter(lambda z: len(z) == k, orbs):
    #         print y,
    #     print
    # print n,[1 for idx in range(n)]
    # print orbs
    orbs.append([])
    # orbs.append([1 for idex in range(n)])
    return orbs

def find_perfect_matching(n):
    """ pair up
    """ 
    orbs = comp_bij(n)
    arr = [[0 for idx in range(len(orbs))] for jdx in range(len(orbs))]
    # first pair up those with at most one element greater than 1
    pairs = [] # the perfect matching
    seen = []  # for keeping track of what in orbits we've 
    # pairs.append([[],[n]])
    # seen.append([])
    # seen.append([n])
    # for x in orbs:
    #     print x
    for i in range(len(orbs)):
        x = orbs[i]
        
        # first we get rid of the special cases
        if len(x) == 0:
            arr[i][orbs.index([n])] = 1
            continue
        if len(x) == 1:
            arr[i][orbs.index([])] = 1
            continue    
        if len(x) == n:
            jlist = filter(lambda y: len(orbs[y]) == n-1, range(len(orbs)))
            if len(jlist) != 1:
                print "Oooops: ",jlist
            else:
                arr[i][jlist[0]] = 1
            continue
        if len(filter(lambda y: y > 1, x)) == 1:
            val = n-len(x)+1
            idx = x.index(val)
            if len(x)%2 == 1: # need to up
                if idx == 0:
                    y = x[:-1]
                else:
                    y = x[:idx-1] + [val+1] + x[idx+1:]
            else:
                y = x[:idx] + [1,val-1] + x[idx+1:]
            arr[i][orbs.index(y)] = 1
            continue
         
        # now we take care of the easy cases
        # first we collect some data
        # print "working on",x
        M,m,Midx,midx,al,bl = get_params(x)    
        # go up as far as possible (absorbing 1s)
    
        if 1 == 1:
            print_all_edges(x,M,m,Midx,midx)
        pos = find_flip_idx(x,M,m,Midx,midx)
        if pos in Midx:
            y = x[:pos] + [1,M-1] + x[pos+1:]
        else:
            if pos == 0:
                y = [M] + x[1:-1]
            else:
                y = x[:pos-1] + [M] + x[pos+1:]
        # print x," --?--> ",y
        cnt = 0
        while y not in orbs and cnt <= len(x):
            y = y[1:] + [y[0]]
            cnt += 1
        arr[i][orbs.index(y)] = 1
        
        # print "%s M=%d m=%d Midx=%s midx=%s al=%s bl=%s pos=%d" % (repr(x),M,m,repr(Midx),repr(midx),repr(al),repr(bl),pos)
        # print "---- ",find_flip_idx(x,M,m,Midx,midx)            
    # for i in range(len(orbs)):
    #     for j in range(len(orbs)):
    #         if arr[i][j] == 1:
    #             print orbs[i]," --> ",orbs[j]
            # orbs.pop(orbs.index(x))
    # print "Symmetric?: ",is_symm(arr) # matrix(ZZ,arr).str()
    

def get_params(x):
    """ get M,m,Midx,...
    """
    M = max(x)
    m = M-1 # max(filter(lambda y: y != M, x))
    Midx = filter(lambda y: x[y] == M, range(len(x)))
    midx = filter(lambda y: x[y] == m, range(len(x)))
    # compute the number of 1s before each type
    al = []
    bl = []
    # print M,m,Midx,midx
    for j in Midx:
        cur = 0
        while x[(j-cur-1+len(x))%len(x)] == 1:
         cur += 1
        al.append(cur)
    for j in midx:
        cur = 0
        while x[(j-cur-1+len(x))%len(x)] == 1:
         cur += 1
        bl.append(cur)
    # print " orig: M=%d m=%d Midx=%s midx=%s al=%s bl=%s" % (M,m,repr(Midx),repr(midx),repr(al),repr(bl))
    if len(filter(lambda y: y%2 == 1, al)) >= 1:
        m = M
        midx = [zz for zz in Midx]
        M += 1
        Midx = []
        bl = [zzz for zzz in al] # filter(lambda y: y%2 == 1, al)
        al = []
    Midx = [Midx[j] for j in filter(lambda y: al[y]%2 == 0, range(len(al)))]
    midx = [midx[j] for j in filter(lambda y: bl[y]%2 == 1, range(len(bl)))]
    al = filter(lambda y: y%2 == 0, al)
    bl = filter(lambda y: y%2 == 1, bl)
    # print "  new: M=%d m=%d Midx=%s midx=%s al=%s bl=%s" % (M,m,repr(Midx),repr(midx),repr(al),repr(bl))
    return M,m,Midx,midx,al,bl 
    

def lex_gt(a,b):
    """ returns true if a >=_lex b; assumes same length
    """
    for i in range(len(a)):
        if a[i] < b[i]:
            return False
        if a[i] > b[i]:
            return True
    # they're equal
    return True
    

def find_flip_idx(x,M,m,Midx,midx):
    """ find index where we should flip the index
    """
    # absorb all of the 1s
    w = [x[i] for i in range(len(x))]
    Mjdx = [Midx[i] for i in range(len(Midx))]
    mjdx = [midx[i] for i in range(len(midx))]
    for j in mjdx:
        # print w,"j is: ",j
        if j == 0:
            w = w[:len(w)-2] + [M]
        else:
            w = w[:j-1] + [M] + w[j+1:]
        # print w,"new w"
        for k in range(j,len(x)):
            if k in Mjdx and k >= j:
                Mjdx[Mjdx.index(k)] -= 1
            if k in mjdx and k >= j:
                mjdx[mjdx.index(k)] -= 1
    # print "w: ",w
    # print Mjdx,mjdx
    # now figure out which is lex maximal
    z = [w[i] for i in range(len(w))]
    maxz = [w[i] for i in range(len(w))]
    maxidx = 0
    for j in range(1,len(z)):            
        nz = z[j:] + z[:j]
        if lex_gt(nz,maxz):
            maxz = [nz[k] for k in range(len(nz))]
            maxidx = j
    # figure out which position to flip
    j = maxidx
    while j not in Mjdx and j not in mjdx:
        j = (j+1)%len(w)
    if j in Mjdx:
        return Midx[Mjdx.index(j)]
    if j in mjdx:
        return midx[mjdx.index(j)]
    print "oops"

def lexmax_rep(w):
    """
    """
    z = [w[i] for i in range(len(w))]
    maxz = [w[i] for i in range(len(w))]
    maxidx = 0
    for j in range(1,len(z)):            
        nz = z[j:] + z[:j]
        if lex_gt(nz,maxz):
            maxz = [nz[k] for k in range(len(nz))]
            maxidx = j
    return maxz

def is_symm(arr):
    """ check if matrix is symmetric
    """
    for i in range(len(arr)):
        for j in range(len(arr)):
            if arr[i][j] != arr[j][i]:
                return False
        if sum(arr[i]) != 1:
            return False
        if sum([arr[k][i] for k in range(len(arr))]) != 1:
            return False
        if arr[i][i] > 0:
            return False
    return True

def print_all_edges(x,M,m,Midx,midx):
    """ find index where we should flip the index
    """
    # absorb all of the 1s
    w = [x[i] for i in range(len(x))]
    Mjdx = [Midx[i] for i in range(len(Midx))]
    mjdx = [midx[i] for i in range(len(midx))]
    for j in mjdx:
        w = [x[i] for i in range(len(x))]
        # print w,"j is: ",j
        if j == 0:
            w = w[:len(w)-2] + [M]
        else:
            w = w[:j-1] + [M] + w[j+1:]
        print lexmax_rep(x)," -> ",lexmax_rep(w)
    for j in Mjdx:
        w = [x[i] for i in range(len(x))]
        w = w[:j] + [1,M-1] + w[j+1:]
        print lexmax_rep(x)," -> ",lexmax_rep(w)
        
#############################
# Simple example computations
#############################
# Compute the cycle index of a group in various bases
G = PermutationGroup([[(1,4,2,3)],[(1,4),(2,3)]])
cyc = G.cycle_index()
print m(cyc)
print F(cyc)
print s(cyc)
print
# Find a perfect matching for G_5' 
# (printing vertices as compositions rather than as binary words)
print "Perfect matching for G_5'"
find_perfect_matching(5)
print
# Find a bijective FQ-expansion for a given group
print "Find bijective FQ-expansion for",G
ff = open(homepath + 'test','w')
G = PermutationGroup([[(1,2,3,4)],[(1,2),(3,4)]])
find_bij_arb_G(ff,G,4)
ff.close()
# The following will generate the data for Tables 2, 3 and 4
find_bij_alg('table3',3)
find_bij_alg('table4',4)
find_bij_alg('table5',5)