SolverGenerator.txt

read cat(solverFolderName, "/Eqs.txt"): 
read "ExtractMonomialBasis.txt":
read "polysolver.map":
read "ComputeSparsePolyBasis.txt": 
march('open', "convex.mla"):
read "IntegerPoints.txt":
with(convex):
eqs := convert(eqs, list):


selectedhiddenvarind := -1:
if hiddenVarNumber = -1 then
	fr := 1:
	en := numelems(vars):
else
	fr := hiddenVarNumber:
	en := hiddenVarNumber:
end if:

for hiddenvarind from fr to en do

    hiddenVarNumber := hiddenvarind:
    hiddenVar := parse(cat("a", hiddenVarNumber)):
    unHiddenVars := select(v -> v <> hiddenVar, vars):
    printf("Trying out by hidding variable  %s \n.", cat("a", hiddenVarNumber));

    newHiddenVarNumber := max(map(v->parse(substring(v,2..)),vars))+1:
    acthiddenvarnumber := newHiddenVarNumber:
    vars := [cat('a',hiddenVarNumber), op(unHiddenVars),cat('a',newHiddenVarNumber)]:
    eqs := [op(eqs), vars[1] - vars[-1] ]:

    unHiddenVars := vars[..-2]:
    hiddenVarNumber := newHiddenVarNumber:
    hiddenVar := parse(cat("a", hiddenVarNumber)):
    CUDA[Enable](true):
    printf("CUDA enabled : %s \n",CUDA[IsEnabled]());

    if numelems(polycomb) = 0 then polycomb := []: end if:
        randrange := -1000000..1000000:
        randomize():

        Bn := min():
        sparsebasis := []:
        #======================================================================================
        # Iterate through all possible permutations of the variables in question. 
        #======================================================================================

    for varperm in [combinat[permute](vars)[1]] do
            #======================================================================================
            # Convert polynomial equations to matrix form	
            #======================================================================================
    
            M, originalbasis, tempts := monpolymult(eqs, [op(unHiddenVars), hiddenVar], varperm):
    
            numoftempts := numelems(tempts):
            Msymbolic := M:
            #======================================================================================
            # Reduce the matrix to a reduced row echelon form
            #======================================================================================
        Mreduced, reducedeqs, nc := reducepoly(M, originalbasis, noofrowstoreduce, noofdatacoeff, randrange, tempts):
    	numofncs := numelems(nc):

    	#======================================================================================
    	# Random instance for the reduced equations
    	#======================================================================================
    
        for i to noofdatacoeff do  
    		assign(cat('c', i), (i+ithprime(100))/(i-ithprime(10000))):
    	end do:
    	for i to numoftempts do 
    		assign(cat('t', i), tempts[i]): 
    	end do:
    
    	if noofrowstoreduce > 0 then
    		Mpatch:=LinearAlgebra[ReducedRowEchelonForm](convert(Msymbolic[1..noofrowstoreduce,..],Matrix)):
    		M := ArrayTools[Concatenate](1, Mpatch, evalm(convert(Msymbolic[noofrowstoreduce+1..,..], Matrix)) );
    	else
    		M := Msymbolic:
        end:
    	linearMnew := ArrayTools[Reshape](M, numelems(eqs)*numelems(originalbasis), 1):
    	for i to numofncs do 
    		assign(cat('nc', i), linearMnew[nc[i],1]): 
    	end do:
    
    	#======================================================================================
    	# We randomize the hidden variable as well
    	#======================================================================================
        assign(cat('a',hiddenVarNumber), RandomTools[Generate](rational)):
    
        #======================================================================================
        # We find sparse basis for the reduced equations. Here the coefficients are 'nc's.
    	#======================================================================================
        try:
            Toriginal, B, T := extract_mon_basis(reducedeqs, unHiddenVars, sizeofcombs, polycomb, heurisitictemplatesize):
            Tn := convert(map(at -> Size(T[at])[1], [seq(1 .. numelems(T))]), `+`):
            Cr := polytomat(T, B, eval(reducedeqs), unHiddenVars, numelems(B), Tn):           
        	#======================================================================================
        	# We select the newly found basis only if it is smaller than the one previously selected.
            #======================================================================================		
        	if numelems(B) > 0 and numelems(B) < Bn then 
        		sparsebasis := B:
        		Bn := numelems(B):
        		monmultiples := T:
        		selectedhiddenvarind := hiddenvarind:
            end if:
        catch:
            print("error"):
        end try:
        #======================================================================================
        # Unassign everything. As the basis has been found and we do not need sample instance of reduced equations.
        #======================================================================================		
        unassign(cat('a',hiddenVarNumber)):
        for i to numofncs do
    		unassign(cat('nc',i)):
    	end do:
    	for i to numoftempts do
    		unassign(cat('t',i)):
    	end do:	
    	for i to noofdatacoeff do
    		unassign(cat('c',i)):
    	end do:

    end do:
end do:

hiddenVarNumber := acthiddenvarnumber:
hiddenVar := parse(cat("a", hiddenVarNumber)):
unHiddenVars := select(v -> v <> hiddenVar, vars):
vars := [op(unHiddenVars), hiddenVar]:

reducedeqs :=[]:
for ei to numelems(eqs) do 
	reducedeqs := [op(reducedeqs), evalm(Mreduced[ei,..] &* originalbasis)]:
end do:
#======================================================================================
# With the sparse basis obtained we construct a symbolic resultant matrix 
# which is to be solved for eigenvalues.
#======================================================================================		

temp := convert(monmultiples[-1], list, nested = true):
bbs := map(proc (i) options operator, arrow; select(proc (j) options operator, arrow; sparsebasis[j] = [temp[i][1], op(temp[i][2..])] end proc, [seq(1 .. numelems(sparsebasis))])[1] end proc, [seq(1 .. numelems(temp))]):
gbs := convert(convert([seq(1 .. numelems(sparsebasis))], set) minus convert(bbs, set), list):
sparsebasis := sparsebasis[[op(bbs), op(gbs)]]:
Cred := polytomat(convert(monmultiples, list), sparsebasis, reducedeqs, unHiddenVars, Bn, Bn):
degreehidvar := max(map(proc (f) options operator, arrow; degree(f, cat('a',hiddenVarNumber)) end proc, Cred)):
allcs := map(proc (i) options operator, arrow; coeff(Cred, cat('a',hiddenVarNumber)^i) end proc, [seq(1 .. degreehidvar)]):
Cs := convert([coeff(Cred*cat('a',hiddenVarNumber), cat('a',hiddenVarNumber)), op(allcs)], Matrix):

printf("The size of sparse basis is %d.\n", Bn);

#======================================================================================
# Generating the solver and the template to calculate var values from eigen vectors
#======================================================================================
sparsefinalbasis := LinearAlgebra[Transpose](convert(sparsebasis[[seq(1..numelems(bbs))]], Matrix)):
printf("Generating template to be used for extracting values of the un-hidden variables \n");

with(SolveTools): 
unassign('d'):
lambdas := Matrix(Size(sparsefinalbasis)[2], 1, symbol = 'd'):
solForm := Matrix(Size(sparsefinalbasis)[2], Size(vars)[2]-1, symbol = 's'): 
for i to Size(vars)[2]-1 do 
	eVec := Matrix(Size(vars)[2], 1, 0): 
	eVec[i, 1] := 1: 
	polEqs := Matrix([[convert(sparsefinalbasis, Matrix)], [Matrix(1, Size(sparsefinalbasis)[2], 1)]]): 
	temp := LinearAlgebra[LinearSolve](polEqs, eVec, free = 'd'): 
	for j to Size(sparsefinalbasis)[2] do 
		solForm[j, i] := temp[j]:
	end do 
end do: 
for k to ArrayTools[Size](sparsefinalbasis)[2] do 
	assign(lambdas[k, 1], 0): 
end do: 
solForm := evalm(solForm):
################################################################################################
fd := fopen("offline_patch.m", READ,TEXT):
	temp := readbytes(fd, infinity):
fclose(fd):

solverFileName := convert(cat(solverFolderName, "/solver.m"), string):
fd := fopen(solverFileName, WRITE, TEXT):
writeline(fd, "function[PEPsolutions] = solve(data) "):
writeline(fd, cat("hiddenvarnumber = ",selectedhiddenvarind, ";")):
writeline(fd, cat(op(map(i -> cat("c", i, " = data(", i, ");"), [seq(1..noofdatacoeff)])))): 

for i to numoftempts do
	writeline(fd, cat("t",i, " = ",tempts[i] ,";")):
end do:

writeline(fd, convert(cat("M = zeros(",numelems(eqs), ",", numelems(originalbasis),");" ),string)):
for i to numelems(eqs) do
	for j to numelems(originalbasis) do
		if Msymbolic[i,j] <> 0 then
			writeline(fd, cat("M(",i,",",j,") = ", Msymbolic[i,j], ";")):
		end if:
	end do:
end do:

writeline(fd, cat("M = [rref(M(1:", noofrowstoreduce, ",:)); M(", (noofrowstoreduce+1), ":end,:)];")):
writeline(fd, cat(op(map(i -> cat("nc", i, " = M(", nc[i], ");"), [seq(1..numofncs)]))) ): 

writeline(fd, convert(cat("Cs = zeros(",ArrayTools[Size](Cs)[1], ",", ArrayTools[Size](Cs)[2],");" ),string)):
for j to ArrayTools[Size](Cs)[1] do
	for k to ArrayTools[Size](Cs)[2] do
		if Cs[j,k] <> 0 then
		    writeline(fd, convert(cat("Cs(",j, ",", k, ") = ", convert(simplify(Cs[j,k]),string), ";" ),string)):
		end if:
	end do:
end do:
writeline(fd, CodeGeneration[Matlab](convert(solForm,list), resultname="solForm", output=string)):
writeline(fd, cat("allCss = mat2cell(Cs, ", numelems(sparsebasis), ", ones(1,2) * ", numelems(sparsebasis), ");")):
writeline(fd, cat("noOfVars = ",(numelems(vars)-1), "; \nC0 = allCss{1}; \nC1 = allCss{2};")):
writeline(fd, cat("A1 = C0(1:end-", numelems(bbs), ",1:", numelems(bbs), ");")):
writeline(fd, cat("A2 = C0(1:end-", numelems(bbs), ",", (numelems(bbs)+1), ":end);")):
writeline(fd, cat("B1 = C0(end-", numelems(bbs), "+1 : end,1:", numelems(bbs), ");")):
writeline(fd, cat("B2 = C0(end-", numelems(bbs), "+1 : end, ", (numelems(bbs)+1), ":end);")):
writebytes(fd, temp):
writeline(fd, "\nend"):
fclose(fd):
printf("The solver has been generated and saved in %s.\n\n", cat(solverFolderName));
########################################################################################################