functionsSage.sage

import itertools as itt
import time

load("gRaysConesS00.sage")
# this contains: rays37, lineality37, coneReps, Matroids37_orbits

load("Gr37MSD0.sage")
# this contains MSDs, a list of all the matroid subdivisions. 


rays_index={rays37[i]:i for i in range(len(rays37))}

#rays_index={rays[i]:i for i in range(len(rays))}
# A dictionary for fast lookup of index of a ray from "raysS00.sage".



Bases37=[(0,1,2),(0,1,3),(0,2,3),(1,2,3),(0,1,4),(0,2,4),(1,2,4),(0,3,4),(1,3,4),(2,3,4),(0,1,5),(0,2,5),(1,2,5),(0,3,5),(1,3,5),(2,3,5),(0,4,5),(1,4,5),(2,4,5),(3,4,5),(0,1,6),(0,2,6),(1,2,6),(0,3,6),(1,3,6),(2,3,6),(0,4,6),(1,4,6),(2,4,6),(3,4,6),(0,5,6),(1,5,6),(2,5,6),(3,5,6),(4,5,6)];
# This is a list of all triples of numbers (i,j,k) 0<= i<j<k<=6.

def fileToLists(fileName):
    with open(fileName,'r') as inFile:
        fStr = inFile.readlines()
    return [tuple(map(int, cone.split(" "))) for cone in fStr]
# Input: "fileName" is a string for the name of a file. Each line of this file should contain a space separated list of integers.
#Output: returns a list of tuples, each tuple is the integers in a single line. 


def Pluckers(k,n):
# Input: positive integers k < n, Output: list of all k-tuples of [0,1,...,n-1] ordered in revLex
    r = range(n)[::-1]
    c = itt.combinations(r, k)
    l = [i[::-1] for i in c]
    rl = l[::-1]
    a=[]
    for e in rl:
        a.append(tuple(e))
    return a

def groupFromGenerators(n,gens):
# Input: a positive integer "n", and a set "gens" of tuples, each tuple is a permutation of (0,1,...,n-1). 
# Output: a list of tuples, each tuple is a permutation of (0,1,...,n-1). These are the elements of S_n generated by the generators "gens" 
    def one_iteration(n,gens,to_act,group,exhausted):
        if to_act == set():
            return list(group)
        else:
            new=set.union(*[set([tuple([g[s[i]] for i in range(n)]) for s in to_act]) for g in gens])
            group.update(new)
            exhausted.update(to_act)
            to_act = new - exhausted
            return one_iteration(n,gens,to_act,group,exhausted)

    return one_iteration(n,gens,gens.copy(),gens.copy(),gens.copy())

def g_inv(g):
# Input: a tuple "g" which is a permutation of (0,1,...,n-1), for n=len(g).
# Output: a tuple that is the inverse of the permutation g.
    inv={g[i]:i for i in range(len(g))}
    return tuple([inv[gi] for gi in range(len(g))])


def s_on_ijk(s,triple):
# Input: a permutation "s"  and a tuple "triple" (i,j,k) representing a basis of a rank 3 matroid.
# Output: a new tuple permuted by s, namely (s[i],s[j],s[k]) (reordered in increasing order).
    return tuple(sorted([s[triple[0]],s[triple[1]],s[triple[2]]]))

def s_on_ray(s,ray):
# Input: a ray "ray" and a tuple "s" which is a permutation on (0,1,...,len(ray)-1).
# Output: new ray with coordinates permuted by s.
    g=g_inv(s)
    return tuple([ray[g[i]] for i in range(len(ray))])

def s_to_NBases(n,Bases,s):
# Input: a positive integer n (the size of the base set of the matroid), the list of "Bases" of a rank 3 matroid on [n], and a tuple "s" representing a permutation 
#  of (0,1,...,n-1).
# Output: a tuple which is a permutation of (0,1,...,len(Bases)-1) given by the induced permutation of s on the list Bases. More specifically, "Bases" is originally ordered in revlex, so this corresponds to the permutation (0,1,...,len(Bases)-1).  Then we use "s_on_ijk" to let s act on each basis, and record its original position in revlex. This gives a permutation of (0,1,...,len(Bases)-1).
    Bases_index={Bases[i]:i for i in range(len(Bases))}

    def s_on_NBases(s,i):
    # Input: and a tuple "s" representing a permutation of (0,1,...,len(Bases)-1).
    # Output: the position of (s[i],s[j],s[k]) in "Bases".
        p_ijk = Bases[i]
        p_sijk= s_on_ijk(s,p_ijk)
        return Bases_index[p_sijk]

    return tuple([s_on_NBases(s,i) for i in range(len(Bases))])


def s_to_Nrays(n,Bases,s,RAYS):
# Input: a positive integer n (the size of the base set of the matroid), the list of "Bases" of a rank 3 matroid on [n], a tuple "s" representing a permutation 
#  of (0,1,...,n-1), and a list "RAYS" of rays, each ray is a tuple.
# Output: a tuple which is a permutation of (0,1,...,len(RAYS)-1) given by the induced permutation of s on the RAYS. More specifically, "RAYS" is a given fixed list, and his corresponds to the permutation (0,1,...,len(Bases)-1).  Then we use "s_on_ray" to let s act on each ray, and record its position in RAYS. This gives a permutation of (0,1,...,len(RAYS)-1).

    sBases=s_to_NBases(n,Bases,s)
    rays_index = {RAYS[i]:i for i in range(len(RAYS))}

    def g_on_ray_index(g,i):
    # Input: and a tuple "g" representing a permutation of (0,1,...,len(RAYS)-1).
    # Output: the position of s_on_ray(g,ray) in "RAYS".
        ray=RAYS[i]
        gray=tuple(s_on_ray(g,ray))
        return rays_index[gray]

    return tuple([g_on_ray_index(sBases,i) for i in range(len(rays))])

def G_in_NRays(n,Bases,gensG,RAYS):
# Input: a positive integer n (the size of the base set of the matroid), the list of "Bases" of a rank 3 matroid on [n], a set "gensG" tuples "s" which are permutations of (0,1,...,n-1) (This is a set of generators of a group G), and a list "RAYS" of rays, each ray is a tuple.

# Output: a list of tuples giving a group "GROUP." This is a subgroup of the permutations on (0,1,...,len(RAYS)-1) generated by the elements of "gensG" where each such generator is converted to a permutation on (0,1,...,len(RAYS)-1).

    G=groupFromGenerators(n,gensG)

    Bases_index={Bases[i]:i for i in range(len(Bases))}
    rays_index={RAYS[i]:i for i in range(len(RAYS))}

    def s_on_NBases(s,i):
    # Input: and a tuple "s" representing a permutation of (0,1,...,len(Bases)-1).
    # Output: the position of (s[i],s[j],s[k]) in (0,1,...,len(Bases)-1) with respect to revLex.
        p_ijk = Bases[i]
        p_sijk= s_on_ijk(s,p_ijk)
        return Bases_index[p_sijk]

    NBases=len(Bases)
    G_in_NBases=[tuple([s_on_NBases(s,i) for i in range(NBases)]) for s in G]

    def g_on_ray_index(g,i):
    # Input: and a tuple "g" representing a permutation of (0,1,...,len(RAYS)-1).
    # Output: the position of s_on_ray(g,ray) in "RAYS".
        ray=RAYS[i]
        gray=tuple(s_on_ray(g,ray))
        return rays_index[gray]

    GROUP=[tuple([g_on_ray_index(g,i) for i in range(len(RAYS))]) for g in G_in_NBases]

    return GROUP

def g_on_tuple(g,x):
# Input: a tuple "g" that is a permutation of (0,1,...,N-1), and a tuple "x" consisting of numbers from {0,1,...,N-1}, sorted in increasing order.
# Output: a new tuple derived from "x" given by the permutation "g", sorted in increasing order.
    return tuple(sorted([g[i] for i in x]))

def G_orbits(G,X):
# Input: a group "G" as a list of permutations of (0,1,...,N-1) (giving a group under composition), and "X" a set of tuples of range(N), each tuple is ordered in increasing order
# Output: a dictionary, keys = representatives, and values = the orbit of the representative.
    remaining = X.copy()
    orbits = {}
    while len(remaining)>0:
        x=remaining.pop()
        orbits[x]=set([x])
        for g in G:
            gx=g_on_tuple(g,x)
            orbits[x].add(gx)
            remaining.discard(gx)
    return orbits

def act_by_G(G,Xreps):
# Input: a group "G" as a list of permutations of (0,1,...,N-1) (giving a group under composition), and "Xreps" a set of tuples of range(N), each tuple is ordered in increasing order. 
# Output: a dictionary, keys = x in Xreps, and values = the orbit of x by the action of G.
    orbits={}
    for x in Xreps:
        orbits[x]=set([g_on_tuple(g,x) for g in G])
    return orbits

def star(cones, face):
# Inupt: a set "cones" where each cone in cones is a tuple of numbers in increasing number, where each number represents a ray in a list of rays, and a "face" which is a cone represented in a similar way (want this to be lower dimensional cone than those in "cones").
# Output: a subset of "cones" that contain "face" as a face.
# Typical use: *cones* will be the list of *all* maximal cones, so this function returns the list of all maximal cones containing "face". 
    star_face=set()
    for cone in list(cones):
        if set(face)<set(cone):
            star_face.add(cone)
    return list(star_face)

def star_dict(cones, faceReps):
# Inupt: a set "cones" where each cone in cones is a tuple of numbers in increasing number, where each number represents a ray in a list of rays, and a set/list "faceReps" of faces which are cones represented in a similar way (want these to be G-orbit representatives of faces of a given dimension).
# Output: a dictionary, keys = face in faceReps, values = all cones in "cones" containing face.

# Typical use: "cones" will be the set of *all* maximal cones, and faceReps will be representatives (up to symmetry) of all non-maximal cones of a given dimension. 
    stars={}
    for face in faceReps:
        stars[face] = star(cones, face)
    return stars


def test_pair(starDict, rays, lineality, dim, cone, pair):
# Input: "starDict" is a dictionary where the keys are the dimensions of the cones (mod lineality), and the value a dimension is another dictionary similar to the output of "star_dict" when  "cones" is the set of *all* maximal cones, and faceReps will be representatives (up to symmetry) of all non-maximal cones of the given dimension.
#"rays" is a list of the rays of the fan, this should be a list, where each list is a tuple representing a ray, "lineality" is a a list consisting of a basis of the lineality space (similar format to "rays"), "dim" is the dimension of "cone" minus the dimension of the lineality space, "cone" is a tuple of numbers in increasing order, where the number corresponds to the position of ray in "rays", and "pair" is a 2-uple of nonnegative integers from 0,...,len(starDict[dim][cone])
    
# Output: Computes the maximal cones "starDict[dim][cone][pair[i]]" (i=0,1) as Polyhedron, and compute the dimension of the intersection of one with (-1) times the other. Here, the maximal cones are viewed in the the space N_R/span(cone). We account for this by taking the lineality space to be the lineality space of the fan + the span of "cone". 
    star_cone = starDict[dim][cone]
    max_cone_0 = star_cone[pair[0]]
    max_cone_1 = star_cone[pair[1]]
    lineality_space = [vector(l) for l in lineality] + [vector(rays[i]) for i in cone]
    Max_Cone_0 = Polyhedron(rays = [vector(rays[i]) for i in max_cone_0], lines = lineality_space, base_ring=QQ)
    Max_Cone_1 = Polyhedron(rays = [vector(rays[i]) for i in max_cone_1], lines = lineality_space, base_ring=QQ)
    a = Max_Cone_0.intersection(-Max_Cone_1)
    return a.dim()




S7 = groupFromGenerators(7,set([(1,0,2,3,4,5,6), (1,2,3,4,5,6,0)]))
# this is a set of tuples, it consists of all permutations of 0,...,6. 

S7OnRayIndices = G_in_NRays(7,Bases37,set([(1,0,2,3,4,5,6), (1,2,3,4,5,6,0)]),rays37)
# This is the action of S7 on (0,...,len(rays37)-1) i.e. the action of S7 on the ray indices.


S7_in_NBases=[s_to_NBases(7,Bases37,s) for s in S7]
# This is a copy of S7 in S35 given by the action of S7 on the triples ijk for i,j,k in [7].

matroid37_orbit_reps_list = [tuple([i for i in range(35) if Matroids37_orbits[j][i]==1]) for j in range(108)]
# this is a list of the rank 3 matriods on [7], up to the action of S7. The matroids are listed in the order in which they appear in "Matroids37_orbits" in the file "extended-Trop-37-input.sage". However, here each matroid is represented by a tuple of numbers among 0, 1, ..., 34. the number i is in the tuple iff the i-th triple ijk is a basis of the matroid (the triples are ordered in revlex).  


matroid37_orbit_reps = set(matroid37_orbit_reps_list) 
# This is just "matroid37_orbit_reps_list", but as a set. 

matroid37_all = set.union(*act_by_G(S7_in_NBases,matroid37_orbit_reps).values())
# this is the set of all rank 3 matroids on [7]. each matroid is represented as a tuple as in matroid37_orbit_reps.

Matroids_all_bases_list=list(matroid37_all)
# this is just matroid37_all, but as a list.

Matroids_all_bases_list_index = {}
# this is a dictionary, the keys are the matroids as in "Matroids_all_bases_list" and the value at a matroid is its index in "Matroids_all_bases_list".
for i in range(len(Matroids_all_bases_list)):
    Matroids_all_bases_list_index[Matroids_all_bases_list[i]] = i

   
def s_on_matroid_index(n,s,Mi):
# Input: a positive integer n (the size of the base set of the matroid "Mi"), and a tuple "s" representing a permutation of (0,1,...,n-1), and Mi is an integer in range(len(Matroids_all_bases_list)) representing the position of a matroid in Matroids_all_bases_list.
# Output: the index of a matroid in Matroids_all_bases_list given by the action of s on the bases of the matroid corresponding to Mi.    
    M = Matroids_all_bases_list[Mi]
    sB = s_to_NBases(n,Bases37,s)
    sM = tuple(sorted([sB[j] for j in M]))
    return Matroids_all_bases_list_index[sM]


def s_on_MSD_index(n,s,MSDi):
# Input: a positive integer n (the size of the base set of the matroid "Mi"), a tuple "s" representing a permutation of (0,1,...,n-1), a tuple of strictly increasing numbers in range(len(Matroids_all_bases_list)), each number is the index of a matroid in Matroids_all_bases_list. Typical use: the matroids corresponding to Mi are the matroids of a matroid subdivision.
# Output: This applies  s_on_matroid_index(n,s,Mi) to each Mi in MSDi, returns a tuple in increasing order.
    return tuple(sorted([s_on_matroid_index(n,s,Mj) for Mj in MSDi]))


def s_on_ray_index(n,s,i):
# Input: a positive integer n (=len(s)),  a tuple "s" representing a permutation of (0,1,...,n-1), an integer i in range(len(rays37)).
# Output: This permutes the coordinates of rays37[i] via s, and returns the index of the corresponding ray in rays37. 
    sB = s_to_NBases(n,Bases37,s)
    r = rays37[i]
    sr = s_on_ray(sB,r)
    return rays_index[sr]
    

def s_on_cone(n,s,c):
    # Input: a positive integer n (=len(s)),  a tuple "s" representing a permutation of (0,1,...,n-1), tuple "cone" of integers among range(len(rays37)) in increasing order representing a cone. 
    # Output: This permutes the coordinates of rays37[i] via s, and returns the index of the corresponding ray in rays37. 
    return tuple(sorted([s_on_ray_index(n,s,i) for i in c]))


def coneSpan(rays, lineality, cone):
    return span(QQ, [rays[i] for i in cone]+lineality)



def intersectLinearSpansRelDim(starDict, rays, lineality, sigma):
    sigmaSpan = coneSpan(rays,lineality,sigma)
    sigmaSpanDim = sigmaSpan.dimension()
    maxCones = starDict[sigma]
    intersectConesSpan = coneSpan(rays,lineality,maxCones[0])
    i=0
    while intersectConesSpan.dimension() > sigmaSpanDim and i<len(maxCones):
        newConeSpan = coneSpan(rays,lineality,maxCones[i])
        intersectConesSpan = intersectConesSpan.intersection(newConeSpan)
        i+=1
    return intersectConesSpan.dimension() - sigmaSpanDim




MSDs_index = [ tuple(sorted([ Matroids_all_bases_list_index[tuple([Bases37.index(B) for B in M])] for M in msd]))  for msd in MSDs]
# this is a list of all the matroid subdivisions from MSDs, each matroid subdivision is a list of their indices in the list Matroids_all_bases_list


coneRepsAll = []
for i in range(1,7):
    coneRepsAll += coneReps[i]
# this the list of all cone representatives as tuples of rays37 indices. Note, this list is in the same order as MSDs, i.e., MSDs_index[i] is the matroid subdivision corresponding to coneRepsAll[i].