=========
Maximum Likelihood Degree
=========

Below is the computation (in Sage) for ML-degree for linear concentration models. It computes the ML-degree for general subspaces of the space of n x n symmetric matrices for n up to 8.


# Cartan torus fixed points in Gausian Moduli and their compasses via recursive built
# and calculation of ML-degree
# divided into sage cells in case yiuwant to run it in a notebook
#################################################
# CELL 1, initialization
# setting the value of M
M=7
# a path = a H<SL(V) fixed point in space of orbits of C* = sequence of C*xH fixed points of invariant sections of O(1) on LG(m,2m)
# recursive construction of paths and their compasses i.e weights of H on co-normal
# P = [[path],[compass]], m+1  = length of path = number of fixed points/sections; binom{m+1}{2}-1 = number of elements in the compasss 
# nesting: [list for m = [ pair=  [path=[list of m+1 fixed pts[m weights]]], [compass=[list of binom{m+1}{2}-1 [m weights]]] end of pair] end for m]
# checking: sum of entries in a vector in the compass = 0, sum in the i-th entry of the path = -m + 2*i, i = 0..m
# first entries are constructed manually
P=[]
P.append([[], []])                                    # m = 0, empty path, empty compass, recursion does not go that deep, it is because 1-th element is 0
P.append([ [[ [-1],[1] ],  [[]]] ])                    # m = 1, one path, empty compass
# the first non-trivial case is m=2, see example 4.2: GM = P^2 
P.append([ [ [[-1,-1],[1,-1],[1,1]], [[-1,1],[-2,2]]],   [ [[-1,-1],[-1,1],[1,1]], [[1,-1],[2,-2]] ] ,   [ [[-1,-1],[0,0],[1,1]], [[-1,1],[1,-1]]]   ])
##
print "calculating ML degree using Gauss Moduli for m=", M
print "\ninitializing the list L with L[m] containing information for GM_m"
# uncomment this part if you want to see the first result
print "m = 2 done, with 3 elements in the grid (fixed points + compasses)"
#for path in P[2]:
#    print "fixed point", path[0], "with compass", path[1]
# the first and the last coordinate of the fixed point is negligible to go on with induction
#################################################
# CELL 2, construction for m=3, trivial loop left to be used later
# DO NOT run this cell twice as this would destroy recursion, re-run from the beginning if needed
# initialization
LT=[]
for m in [3]:
    b=[]  # beginning point of every part = m copies of -1
    for i in range(m):
        b.append(-1)
#    ero=[] #zero vector of length m-4
    zero=[] # zero vector of length m-3
    zzero=[] # zzero vector of length m-2
    if m > 4:
        for i in range(m-4):
            ero.append(0)
    for i in range(m-3):
        zero.append(0)
    for i in range(m-2):
        zzero.append(0) 
        
######################        
############ constructing paths starting with line segments        
    LT=[] # paths with a line at the beginning
    for i in [0..m-1]: # will insert the weight 1 at i-th place
        for t in P[m-1]: # take a [path, compass] of length m-1
            v=[b] # insert the beginning of the new path v
            for s in t[0]: #extend the path by adding 1 at i-th place at every fixed point    
                w=copy(s)
                w.insert(i,1)
                v.append(w)
                
### now constructing compasses                
            cv=[] # compass of the path v
            
            for cs in t[1]: #extend the compass by adding 0 at i-th place at every fixed point   
                cw=copy(cs)
                cw.insert(i,0)
                cv.append(cw)
               
            for jj in [0..m-2]: #type A vector in the compass, (m-2) copies
                typeA=copy(zzero)
                typeA.insert(jj,1)
                typeA.insert(i,-1)        
                cv.append(typeA)
             
            typeB=copy(zzero)            
            for jjj in [0..m-2]: #constructing type B vector in the compass, 1 copy  
                if t[0][1][jjj]==1:
                    typeB.insert(jjj,2)
                    typeB.insert(i,-2)        
                    cv.append(typeB)  
                     
            typeC=copy(zero) # this works for m=3 only
            typeC=[1,1]
            if t[0][1][1]==0:
                typeC.insert(i,-2)
                cv.append(typeC)                       
            LT.append([v,cv])
   
#### CT paths starting with a "conic", adding two zeroes and then two ones at ith and jth place of a path of length m-2
    CT=[] # 
    for i in [0..m-2]:
        for j in [0..i]:
            for t in P[m-2]: # take a  pair [path, compass] in P[m-2]
                v=[b]   # initialize a new path
                w=copy(t[0]) # take path
                w1=copy(w[0])
                w1.insert(i,0)
                w1.insert(j,0)
                v.append(w1)
                for k in [0..len(w)-1]:
                    wk=copy(w[k])
                    wk.insert(i,1)
                    wk.insert(j,1)
                    v.append(wk)
                ########## path completed
                
                cv=[]
       
                if len(t[1])>1:
                    for cs in t[1]: #extend the compass by adding 0 at i-th and j-th place, binom{m-1}{2}-1 copies
                        cw=copy(cs)
                        cw.insert(i,0)
                        cw.insert(j,0)
                        cv.append(cw)
          # adding -1 at i or j and 1 at some other one, 2(m-2) copies, type A
                for k in range(len(zero)+1):
                    cwwi=copy(zero)
                    cwwj=copy(zero)
                    cwwi.insert(k,1)
                    cwwi.insert(i,-1)
                    cwwi.insert(j,0)
                    cwwj.insert(k,1)
                    cwwj.insert(i,0)
                    cwwj.insert(j,-1)
                    cv.append(cwwi)
                    cv.append(cwwj)
                cwwij=copy(zzero) # adding (1,-1) and (-1,1)
                cwwij.insert(i,-1)
                cwwij.insert(j,1)
                cv.append(cwwij)
                cwwji=copy(zzero)
                cwwji.insert(i,1)
                cwwji.insert(j,-1)
                cv.append(cwwji)
                #################
                typeB=copy(zero)            
                for jjj in [0..m-3]: #constructing type B vector in the compass 
                    if t[0][1][jjj]==1:
                        typeB.insert(jjj,2)
                        typeB.insert(i,-1) 
                        typeB.insert(j,-1)      
                        cv.append(typeB)   
                #type C for conics does not exists for m = 3
                
                CT.append([v,cv])
P.append(LT+CT)
######################
print "m = 3, done, with 12 elements in the grid"
# check, uncomment this part if you want to see an example
#print "for m = 3 the points in GM and their compasses are as follows"
#for path in P[3]:
#    print "fixed point\n", path[0], "\nwith compass\n", path[1]
#####################################
# CELL 3, main loop
#######################
# now we can do the recursive construction for your choice of M, realistically not bigger than 9
# initial variables first
for m in [4..M]:
    b=[]  # beginning point of every part = m copies of -1
    for i in range(m):
        b.append(-1)
    ero=[] #zero vector of length m-4
    zero=[] # zero vector of length m-3
    zzero=[] # zzero vector of length m-2
    for i in range(m-4):
        ero.append(0)
    for i in range(m-3):
        zero.append(0)
    for i in range(m-2):
        zzero.append(0) 
        
######################        
############ constructing paths starting with line segments        
    LT=[] # paths with a line at the beginning
    for i in [0..m-1]: # will insert the weight 1 at i-th place
        for t in P[m-1]: # take a [path, compass] of length m-1
            v=[b] # insert the beginning of the new path v
            for s in t[0]: #extend the path by adding 1 at i-th place at every fixed point    
                w=copy(s)
                w.insert(i,1)
                v.append(w)
                
### now constructing compasses                
            cv=[] # compass of the path v
            
            for cs in t[1]: #extend the compass by adding 0 at i-th place at every fixed point   
                cw=copy(cs)
                cw.insert(i,0)
                cv.append(cw)
               
            for jj in [0..m-2]: #constructing type A vector in the compass, (m-2) copies
                typeA=copy(zzero)
                typeA.insert(jj,1)
                typeA.insert(i,-1)        
                cv.append(typeA)
             
            typeB=copy(zzero)            
            for jjj in [0..m-2]: #constructing type B vector in the compass, 1 copy  
                if t[0][1][jjj]==1:
                    typeB.insert(jjj,2)
                    typeB.insert(i,-2)        
                    cv.append(typeB)   
            
            typeC=copy(zero) # type C
            j0=0
            while j0 <= m-3:
                if t[0][1][j0]!=0:
                    j0+=1
                else:
                    break
            j1=j0+1
            while j1 < m-2:
                if t[0][1][j1]!=0:
                    j1+=1
                else:
                    break
            if j0!=m-2:
                typeC.insert(j0,1)
                typeC.insert(j1,1)
                typeC.insert(i,-2)
                cv.append(typeC)
                    
                   
            LT.append([v,cv])
    
    
#### CT paths starting with a "conic", adding two zeroes and then two ones at ith and jth place of a path of length m-2
    CT=[] # 
    for i in [0..m-2]:
        for j in [0..i]:
            for t in P[m-2]: # take a  pair [path, compass] in P[m-2]
                v=[b]   # initialize a new path
                w=copy(t[0]) # take path
                w1=copy(w[0])
                w1.insert(i,0)
                w1.insert(j,0)
                v.append(w1)
                for k in [0..len(w)-1]:
                    wk=copy(w[k])
                    wk.insert(i,1)
                    wk.insert(j,1)
                    v.append(wk)
                ########## path completed
                
                cv=[]
       
                if len(t[1])>1:
                    for cs in t[1]: #extend the compass by adding 0 at i-th and j-th place, binom{m-1}{2}-1 copies
                        cw=copy(cs)
                        cw.insert(i,0)
                        cw.insert(j,0)
                        cv.append(cw)
          # adding -1 at i or j and 1 at some other one, 2(m-2) copies
                for k in range(len(zero)+1):
                    cwwi=copy(zero)
                    cwwj=copy(zero)
                    cwwi.insert(k,1)
                    cwwi.insert(i,-1)
                    cwwi.insert(j,0)
                    cwwj.insert(k,1)
                    cwwj.insert(i,0)
                    cwwj.insert(j,-1)
                    cv.append(cwwi)
                    cv.append(cwwj)
                cwwij=copy(zzero) # adding (1,-1) and (-1,1)
                cwwij.insert(i,-1)
                cwwij.insert(j,1)
                cv.append(cwwij)
                cwwji=copy(zzero)
                cwwji.insert(i,1)
                cwwji.insert(j,-1)
                cv.append(cwwji)
                #################
                typeB=copy(zero)            
                for jjj in [0..m-3]: #constructing type B vector in the compass 
                    if t[0][1][jjj]==1:
                        typeB.insert(jjj,2)
                        typeB.insert(i,-1) 
                        typeB.insert(j,-1)      
                        cv.append(typeB)   
            
                if m > 3: # this part for m \leq 3 is for nothing
                    typeC=copy(ero) # constructing type C
                    j0=0
                    while j0 <= m-3:
                        if t[0][1][j0]!=0:
                            j0+=1
                        else:
                            break
                    j1=j0+1
                    while j1 < m-2:
                        if t[0][1][j1]!=0:
                            j1+=1
                        else:
                            break
                    if j0!=m-2:
                        typeC.insert(j0,1)
                        typeC.insert(j1,1)
                        typeC.insert(i,-1)
                        typeC.insert(j,-1)
                        cv.append(typeC)
                    
                CT.append([v,cv])

    
#    PPK=LT+CT
#    for path in CT:
#        print path[0]
#        print path[1], "\n ==========="
#        if len(path[1]) < 9:
#            print path[0]
    P.append(LT+CT)
    print "m =", m, "done, with", len(P[m]), "elements in the grid"

print "the data to calculate function phi for m up to", M, "is in the list L"
##########################################################
# CELL 4, flattening the data, downgrading to a sufficiently general C* action 
# choose your m
m=M
# C choice of a C* by projection of the lattice of characters
# you can choose different C but then do the sanisty check below
# a1, a2 pull-backs of O(a1) and O(a2)
dim=(m+1)*m/2 -1 # binom{m+1}{2}-1
C=[]
for i in [0..m-1]:
    C.append(2^i)
A=[1,1]
PPKflat=[]
for path in P[m]:
    Mu1=A[0]*sum(C[i]*path[0][1][i] for i in [0..m-1])
    Mu2=A[1]*sum(C[i]*path[0][m-1][i] for i in [0..m-1])
    Cmps=[]
    for compass in path[1]:
        cflat=sum(C[i]*compass[i] for i in [0..m-1])
        Cmps.append(cflat)
    PPKflat.append([[Mu1,Mu2],Cmps])
print "\ndowngrading action to C* with weights", C, "\nresult in list PPKflat"
#sanity check for another choice of C
#for uv in PPKflat: #checking if the choice is sufficiently general i.e. no zeroes in compass
#    if prod(a for a in uv[1])==0:
#            print "oj"
#################################################################
# CELL 5, calculating phi via equivariant cohomology
def denomhz(a): return prod(u for u in a);
print "\ncalculating function phi(a,m) via equivariant cohomology for m = ", m, "\n===================="
for a in [0..dim]:
     SUM=(-1)^dim*sum((b[0][0])^a*(b[0][1])^(dim-a)/denomhz(b[1]) for b in PPKflat);
     W=expand(factor(SUM));
     print "a = ", a, " phi(a,m) = ", W, "\n================"
###########################################################