ExtractMonomialBasis.txt

extract_mon_basis:= proc(eqs, vars, sizeofcombs, polycomb, heurisiticbasislength) 
local tempT, zcols, Cr, Ctemp, rows2rem, rows2kp, rowindices, colindices, cstried, c, Btemp, Ttemp, B1, T1, T2, B2, temp, bbs, gbs, cols2kp, cols2rem, simplex, selectedB, selectedT, counter, newpredeigsz, oldpredeigsz, modifiedCred, found, ti, crowstobedelet, rowColCnts, tempCred, Cred, tn, C, solverrank, testcoeffmat, Tn, ii, baks, k, T, Bset, minmonnorm, isbnpresent, Bn, B, interiorintpoints, hp1s, hlfsps, translatedconvhull, delta, deltatoskip, deltasiterator, combind, i, allBs, minkowskis, ind, PForMink, comb, vci, variouscombinations, A, v, P, EE, ll, indicestotest, indicestoskip:

EE := eval(expand(eqs)):
for i to ArrayTools[Size](EE)[2] do 
	P[i] := newtonpolytope(EE[i], vars):
	v[i] := vertices(P[i]):
	A[i] := lattpts(EE[i], vars):     
end do: 
simplex := map(proc (j) options operator, arrow; map(proc (i) options operator, arrow; `if`(i = j, rand(1..4)(), 0) end proc, [seq(1 .. numelems(vars))]) end proc, [seq(1 .. numelems(vars)+1)]):
P[i] := convhull(op(simplex)):
v[i] := vertices(P[i]):
A[i] := vertices(P[i]):
v := convert(v, list):

if numelems(polycomb) = 0 then 
    variouscombinations := (map(proc (l) options operator, arrow; seq(combinat[choose]([seq(1 .. numelems(P))], l)) end proc, sizeofcombs)):
else
    variouscombinations:= [polycomb]: 
end if:

unassign('minkowskis'); 
for vci to numelems(variouscombinations) do 
	comb := variouscombinations[vci]:
	PForMink := []; 
	for ind in comb do 
		PForMink := [op(PForMink), P[ind]]: 
	end do:
	minkowskis[vci] := minkowskisum(op(PForMink)): 
end do:

selectedB := -1:
selectedT := -1:
allBs := []:
oldpredeigsz := infinity:
for combind in combinat[randperm]([seq(1..numelems(variouscombinations))]) do
	randomize():
	printf("Considering the polytopes : %s \n", convert(variouscombinations[combind],string));
    deltasiterator := combinat[cartprod](map(proc (i1) options operator, arrow; map(proc (j) options operator, arrow; j end proc, [ -0.01, 0.01, 0]) end proc, [seq(1 .. numelems(vars))])):	
    while not deltasiterator[finished] do
       delta := deltasiterator[nextvalue]():
        
		translatedconvhull := minkowskisum(minkowskis[combind], convhull(convert(delta, rational))):
		if dim(translatedconvhull) < numelems(vars) then 
    			next: 
		end if:     
		hp1s := map(rel -> convert(`~`[`*`](lhs(rel), ` $`, vars), `+`) <= rhs(rel), convert(hspaces(translatedconvhull), list)): 

		try
			B := convert(IntegerPoints1(hp1s, vars), list); 
		catch:
			printf("Failed for delta %s \n", convert(delta,string));
			next:
		end try:
		Bn := numelems(B):
        if  Bn = 0 then next: end if:
		
        minmonnorm := norm(convert(ListTools[FindMinimalElement](B, proc (b1, b2) options operator, arrow; norm(convert(b1-b2,Vector)) = 0 end proc),Vector)):
		# if minmonnorm <> 0 then next: end if:

		Bset := convert(B, set); 
		for k to numelems(EE) do 
			T[k] := []:
			baks := map(proc (belement) options operator, arrow; op(select(proc (diffba) options operator, arrow; foldl(`and`, true, op(`~`[`>=`](diffba, ` $`, map(proc (ii) options operator, arrow; 0 end proc, [seq(1 .. numelems(vars))])))) end proc, map(proc (aelement) options operator, arrow; `~`[`-`](belement, aelement) end proc, A[k]))) end proc, B):
			for ii to numelems(A[k])-1 do 
				baks := ListTools[FindRepetitions](baks):
			end do: 
			T[k] := convert(convert(convert(baks, set), list), Matrix): 
		end do:
		T := convert(T,list):
		Tn := convert(map(t->ArrayTools[Size](t,1), T), `+`):
        
		if Tn >= Bn and min(map(t->ArrayTools[Size](t,1), T)) > 0 then 
		        temp := convert(T[-1], list, nested = true): 
        		bbs := map(proc (i1) options operator, arrow; select(proc (j1) options operator, arrow; B[j1] = [temp[i1][1], op(temp[i1][2..])] end proc, [seq(1 .. Bn)])[1] end proc, [seq(1 .. numelems(temp))]): 
		        gbs := convert(convert([seq(1 .. numelems(B))], set) minus convert(bbs, set), list):   
                testcoeffmat := polytomat(T, B, EE, vars, Bn, Tn):			
        	    testcoeffmat := LinearAlgebra[Modular][Mod](ithprime(1000), testcoeffmat, integer[8]):
                solverrank := LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat):

	            if Bn = solverrank then
                    newpredeigsz := Bn - LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat[1..(Tn - numelems(temp)), ..]):

    		        if newpredeigsz <= oldpredeigsz and convert(map(b->numelems(b)>Bn, allBs), `and`) then
                        oldpredeigsz := newpredeigsz:

                        if  LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat[1..Tn - numelems(temp), gbs]) = numelems(gbs) then                                                                                    
                            selectedB := B:
                            selectedT := T:
                            allBs := [op(allBs), B]:
                            printf(" .................... \n"):
                            printf("Coefficient matrix: %d X %d \n", Tn, Bn); 
                            printf("Predicted eigenvalue problem: %d \n", newpredeigsz):                            
            	        	printf("delta =  %{t}2f \n", LinearAlgebra[Transpose](convert(delta,Matrix))):
                            printf(" .................... \n"):
                            
                        end if:
                    end if:
                end if:
    		end if:
        end do:
    end do:

if selectedT = -1 then 
	printf("NO SOLVER FOUND .... \n"):
	return convert([],list), convert([],list), [[]]:
end if:

T := selectedT:
Tn := convert(map(proc (t) options operator, arrow; ArrayTools[Size](t,1) end proc, T), `+`):
B := selectedB:
Bn := numelems(B):
C := polytomat(convert(T, list), B, eval(eqs), vars, Bn, Tn); 
testcoeffmat := LinearAlgebra[Modular][Mod](ithprime(1000), C, integer[8]):
solverrank := LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat):
printf("Monomial basis extracted. \n"):
temp := convert(T[-1], list, nested = true):
oldpredeigsz := Bn - LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat[1..(Tn - numelems(temp)), ..]):

# Removal of rows and columns =================================================================================================================
B1 := B: T1 := convert(T,list):
T2 := T1: B2 := B1:
found := numelems(B):
    while numelems(B2) > heurisiticbasislength do
    B1 := B: T1 := convert(T,list):
    Tn := convert(map(proc (t) options operator, arrow; ArrayTools[Size](t,1) end proc, T1), `+`):    
    cstried :=[]:
    printf("Attempting row-column removal... \n"):    
    ind := 1;    
    
    while numelems(B1) > heurisiticbasislength do
    	randomize():
        Tn := convert(map(proc (t) options operator, arrow; ArrayTools[Size](t,1) end proc, T1), `+`):
    	Ctemp := polytomat(T1, B1, eval(eqs), vars, numelems(B1), convert(map(proc (t) options operator, arrow; ArrayTools[Size](t,1) end proc, T1), `+`)):
        Ctemp := LinearAlgebra[Modular][Mod](ithprime(1000), Ctemp, integer[8]):
    	Btemp := B1:
    	Ttemp := T1:
        
        temp := convert(Ttemp[-1], list, nested = true):
    	bbs := map(i1 -> select(k -> Btemp[k] = [temp[i1][1], op(temp[i1][2..])], [seq(1 .. numelems(Btemp))])[1], [seq(1 .. numelems(temp))]):
    	gbs := convert(convert([seq(1 .. numelems(Btemp))], set) minus convert(bbs, set), list):    
        
        if numelems(cstried) = numelems(B) then break: end if:
        rowindices := [seq(1..Tn)]:
        colindices := [seq(1..numelems(Btemp))]:
					
        c := [op(combinat[randperm]([op(gbs), op(bbs)])[1..1])]:
        cstried	:= [op(cstried), op(c)]:
        

    	# Cols to remove
        rows2rem := convert(convert(map(j -> op(select(i -> Ctemp[i,j] <> 0, [seq(1..ArrayTools[Size](Ctemp)[1])])), c),set),list):	
        rows2kp := select(proc (j) options operator, arrow; not(j in rows2rem) end proc, [seq(1 .. ArrayTools[Size](Ctemp,1))]): 		
                
        # Modified C
        Ctemp := Ctemp[rows2kp, () .. ()]:                
        if ArrayTools[Size](Ctemp,1) < ArrayTools[Size](Ctemp,2) or ArrayTools[Size](Ctemp,1) = 0 then next: end if:
        
        zcols := convert(ArrayTools[AddAlongDimension](map(i -> `if`(i = 0, 0, 1), Ctemp), 1), list):
        cols2kp := select(j -> zcols[j] <> 0, [seq(1 .. numelems(zcols))]):
        cols2rem := select(j -> zcols[j] = 0, [seq(1 .. numelems(zcols))]):     
        Ctemp := Ctemp[() .. (), cols2kp]: 
        Btemp := Btemp[cols2kp]:                
        if numelems(cols2rem) = 0  then next: end if:        

        rowindices := rowindices[rows2kp]:
        Ttemp:=map(k -> Ttemp[k][select(j -> j+convert(map(t->ArrayTools[Size](Ttemp[t],1), [seq(1 .. numelems(Ttemp[..k-1]))]), `+`) in rowindices, [seq(1..ArrayTools[Size](Ttemp[k],1))]),..], [seq(1..numelems(Ttemp))]):	
    	if convert(map(t->ArrayTools[Size](t,1), Ttemp), `*`) = 0  then next: end if:               

        temp := convert(Ttemp[-1], list, nested = true):
        bbs := map(i1 -> select(k -> Btemp[k] = [temp[i1][1], op(temp[i1][2..])], [seq(1 .. numelems(Btemp))])[1], [seq(1 .. numelems(temp))]):
        gbs := convert(convert([seq(1 .. numelems(Btemp))], set) minus convert(bbs, set), list):

        if (numelems(gbs)-LinearAlgebra[Modular][Rank](ithprime(1000), Ctemp[[seq(1..convert(map(t -> ArrayTools[Size](t,1), Ttemp[..-2]), `+`))], gbs]) <> 0)  then next: end if:

        newpredeigsz := numelems(Btemp) - LinearAlgebra[Modular][Rank](ithprime(1000), Ctemp[1..convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, Ttemp[..-2]), `+`), ..]):        
        if newpredeigsz > oldpredeigsz then next: end if:
        B1 := Btemp:
        T1 := Ttemp:    
        cstried := []:            
        
    end do:

    if found > numelems(B1) then 
        T2 := T1:
        B2 := B1:
        found := numelems(B1):        
    end if:
end do:

T:=T2: B:=B2:
Tn := convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, T), `+`):
C := polytomat(convert(T, list), B, eval(eqs), vars, numelems(B), Tn): 
printf("Reduced coefficient matrix: %d X %d \n", Tn, numelems(B)):
# =============================================================================================================
# =============================================================================================================
# For tracking an indice, we add '1' to the end of each element in the list of monomial multiples
tn := 0:

for k to  ArrayTools[Size](T)[2] do
	T[k] := convert((map(proc (t) options operator, arrow; [op(convert(T[k][t], list)), t+tn] end proc, [seq(1 .. ArrayTools[Size](T[k])[1])])), Matrix):
	tn := tn + ArrayTools[Size](T[k])[1]:
end do:
# =============================================================================================================

printf("Removing the excess rows to get a reduced matrix with smaller basis.\n"):
# =============================================================================================================
# =============================================================================================================

for counter from 1 to 500 do
    randomize():
    Cred := LinearAlgebra[Copy](LinearAlgebra[Modular][Mod](ithprime(1000), C, integer[8])):
    T1 := T:
    tempCred := Cred:
    rowColCnts := ArrayTools[Size](Cred):
    crowstobedelet := []:

    ll := max(map(proc (ti) options operator, arrow; ArrayTools[Size](ti)[1] end proc, T1)):
    indicestotest := [seq(1..Tn)]: 
    temp := ListTools[Reverse](ListTools[Reverse]([seq(1..convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, T1), `+`))])[1..ArrayTools[Size](T1[-1],1)]):
        
    indicestoskip := []:
    #indicestoskip := [op(indicestoskip),op(select(i -> convert(convert(T1[-1][i  - (tn-ArrayTools[Size](T1[-1],1)) ,..-2], list), `+`) <= 1, temp))]:

    temp := convert(convert(temp, set) minus convert(indicestoskip,set), list):
    indicestotest := convert(convert(indicestotest,set) minus convert(temp,set), list):
    indicestotest := [op(ListTools[Reverse](temp)), op(indicestotest[1..0]), op(combinat[randperm](indicestotest[1..]))]:
    
    for i in indicestotest do 
	for ti to ArrayTools[Size](T1)[2] do 
		temp := T1[ti]; 
		found := select(proc (ti1) options operator, arrow; T1[ti][ti1][-1] = i end proc, [seq(1 .. ArrayTools[Size](T1[ti])[1])]); 
		if numelems(found) > 0 then 
            temp := LinearAlgebra[DeleteRow](T1[ti], found); 
            break: 
		end if: 
	end do; 
	modifiedCred := LinearAlgebra[DeleteRow](tempCred, [op(crowstobedelet), i]):
	testcoeffmat := LinearAlgebra[Modular][Mod](ithprime(1000), modifiedCred, integer[8]):
	solverrank := LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat):
    
	if ArrayTools[Size](temp)[1] > 0 then 
		if solverrank = rowColCnts[2] then
			
			tempT := T1:
			tempT[ti] := LinearAlgebra[DeleteRow](tempT[ti], found):
			temp := convert(tempT[-1], list, nested = true): 
			bbs := map(proc (i1) options operator, arrow; select(proc (j1) options operator, arrow; B[j1] = [temp[i1][1], op(temp[i1][2..-2])] end proc, [seq(1 .. numelems(B))])[1] end proc, [seq(1 .. numelems(temp))]): 
			gbs := convert(convert([seq(1 .. numelems(B))], set) minus convert(bbs, set), list):
            			
			if numelems(gbs) = LinearAlgebra[Modular][Rank](ithprime(1000), testcoeffmat[1 .. convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, tempT[() .. -2]), `+`), gbs]) then
				Cred := modifiedCred:
				T1[ti] := LinearAlgebra[DeleteRow](T1[ti], found):
				crowstobedelet := [op(crowstobedelet), i]:
				rowColCnts := ArrayTools[Size](Cred):
			end if:
		end if:		
	end if:
end do:

    for k to ArrayTools[Size](T)[2] do
	    T1[k] := convert((map(proc (t) options operator, arrow; [op(convert(T1[k][t][..-2], list))] end proc, [seq(1 .. ArrayTools[Size](T1[k])[1])])), Matrix):
    end do:

    temp := convert(T1[-1], list, nested = true):
    bbs := map(i1 -> select(j1 ->  B[j1] = [temp[i1][1], op(temp[i1][2..])], [seq(1 .. numelems(B))])[1], [seq(1 .. numelems(temp))]): 
    gbs := convert(convert([seq(1 .. numelems(B))], set) minus convert(bbs, set), list):
    Cr := polytomat(convert(T1, list), B, eval(eqs), vars, numelems(B), convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, T1), `+`)):
    temp := Cr[1 .. convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, T1[() .. -2]), `+`), gbs]:
    T2 := T1:
    if ArrayTools[Size](Cred,1) = ArrayTools[Size](Cred,2) then     
        break: 
    end if:
end do:
T := T2:

#print(T, numelems(B), "........................................");
Tn := convert(map(proc (t) options operator, arrow; ArrayTools[Size](t,1) end proc, T), `+`):

# =============================================================================================================
Cr := polytomat(convert(T, list), B, eval(eqs), vars, numelems(B), convert(map(proc (t) options operator, arrow; ArrayTools[Size](t)[1] end proc, T), `+`)):

return selectedT, B, T:
end proc: