read "itall.txt":
with(queue):

doBredth:=proc(T,vert,ver,fd,lvl,phatlvl)
	local Q,vertex,currentdist,currentlist,val,j,numneighbors;
	fprintf(fd,"\\page\n");
	Q:=new();
        for vertex in ver do
                T[vertex,11]:=false:
        end do:

        currentdist:=-1;
        currentlist:=[];

	enqueue(Q,[vert,0]);

	fprintf(fd,"%d phat %d dirichlet min maxes function %a \\par\n",phatlvl,lvl,vert);

        while not empty(Q) do
                val:=dequeue(Q);
                if not (val[2]=currentdist) then
                        if(currentdist>-1) then
                                fprintf(fd,"dist %d=[%20.20g,%20.20g] \\par\n",currentdist,min(op(currentlist)),max(op(currentlist)));

                	end if;
                        currentlist:=[];
                        currentdist:=val[2];
                end if;
                currentlist:=[op(currentlist),T[val[1],13]];



                if nops(val[1])<2 then
                        numneighbors:=2;
                else
                        numneighbors:=4;
                end if;
                for j from 1 to numneighbors do
                        if not T[T[val[1],j],11] then
                                enqueue(Q,[T[val[1],j],val[2]+1]);
                                T[T[val[1],j],11]:=true;
                        end if;
                end do;
        end do;
end;

cmp_list:=proc(a,b)
	if a[1]<b[1] then return true;
	else return false;
	end if;
end;

lowerlvl:=3;
upperlvl:=4;
lvl:=upperlvl;
phatlvl:=8;
BMTOL:=.0000000000000001;
numrandfuncs:=50;
lowerpath:=cat("../data/level",lowerlvl,"/");
upperpath:=cat("../data/level",upperlvl,"/");
path:=upperpath;

T:=gasket(phatlvl):

ver:=ourVertices(phatlvl):
make_neighbors(T,phatlvl,ver);

lowerevals:=getEigenVals(lowerpath,lowerlvl);
upperevals:=getEigenVals(upperpath,upperlvl);

lowernumberevals:=[];
uppernumberevals:=[];

for i from 1 to nops(lowerevals) do
	lowernumberevals:=[op(lowernumberevals),[lowerevals[i],i]];
end do:

for i from 1 to nops(upperevals) do
	uppernumberevals:=[op(uppernumberevals),[upperevals[i],i]];
end do:

lowernumberevals:=sort(lowernumberevals,cmp_list):
uppernumberevals:=sort(uppernumberevals,cmp_list):

loweruniquelist:=[[lowernumberevals[1][1],[lowernumberevals[1][2]]]]:
upperuniquelist:=[[uppernumberevals[1][1],[uppernumberevals[1][2]]]]:

for i from 2 to nops(lowernumberevals) do
	if(abs(lowernumberevals[i][1]-loweruniquelist[nops(loweruniquelist)][1])<.00000000000000001) then
		loweruniquelist[nops(loweruniquelist)][2]:=[op(loweruniquelist[nops(loweruniquelist)][2]),lowernumberevals[i][2]];
	else
		loweruniquelist:=[op(loweruniquelist),[lowernumberevals[i][1],[lowernumberevals[i][2]]]];
	end if;
end do:

for i from 2 to nops(uppernumberevals) do
	if(abs(uppernumberevals[i][1]-upperuniquelist[nops(upperuniquelist)][1])<.00000000000000001) then
		upperuniquelist[nops(upperuniquelist)][2]:=[op(upperuniquelist[nops(upperuniquelist)][2]),uppernumberevals[i][2]];
	else
		upperuniquelist:=[op(upperuniquelist),[uppernumberevals[i][1],[uppernumberevals[i][2]]]];
	end if;
end do:


for loop from nops(loweruniquelist)+1 to nops(upperuniquelist) do

template:=fopen("../data/template.rtf",READ);
fd:=fopen(cat(path,phatlvl,"phat",lvl,"dirdata",loop-nops(loweruniquelist),".rtf"),WRITE);
intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals",loop-nops(loweruniquelist),".txt"),WRITE);
fprintf(intdata,"l1, l2, and integral data for the %d phat %d extra dirichelt interpolation %d\n",phatlvl,lvl,loop-nops(loweruniquelist));
fclose(intdata);


ch:=fscanf(template,"%c")[1]:
while ch<>"#" do
	fprintf(fd,"%c",ch):
	ch:=fscanf(template,"%c")[1]:
end do:
fclose(template);

#fprintf(fd,"hello \\page 1: 12321=1231 \\par 12321=12321 \\page");

vl:=ourVertices(lvl):
vp:=ourVertices(phatlvl):
for i from 1 to nops(vl) do
	vp:=subsop(1=NULL,vp):
end do:
vl:=subsop(1=NULL,2=NULL,3=NULL,vl):

avel1:=0;
avel2:=0;
aveinteg:=0;
minl1:=1000000000;
maxl1:=0;

for i from 1 to numrandfuncs do
	clearFunc(T,phatlvl,13);
	funcnum:=(rand() mod nops(vp))+1;
		for ulambda from 1 to loop do
			for fnumber in upperuniquelist[ulambda][2] do
				print("fnumber=");
				print(fnumber);
				print("ulambda=");
				print(upperuniquelist[ulambda][1]);
				readFunc(T,12,cat(path,phatlvl,"phat",lvl,"eigen.",fnumber));
				co:=T[vp[funcnum],12]*(3^(lvl+1))/(2*calcbm(upperuniquelist[ulambda][1],[],BMTOL));
				if (abs(co) > 0) then
					funcSumWithConstant(T,phatlvl,co,13,12,13);
				end if;
			end do;
		end do;

	funcProduct(T,phatlvl,13,13,12);
        l2:=simp_rule(T,phatlvl,12);
        funcAbs(T,phatlvl,13,12);
        l1:=simp_rule(T,phatlvl,12);
        integ:=simp_rule(T,phatlvl,13);
	intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals",loop-nops(loweruniquelist),".txt"),APPEND);
	fprintf(intdata,"%a:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",vp[funcnum],l1,l2,integ);
	fclose(intdata);
	avel1:=avel1+l1;
	avel2:=avel2+l2;
	aveinteg:=aveinteg+integ;
	minl1:=min(minl1,l1);
	maxl1:=max(maxl1,l1);

	doBredth(T,vp[funcnum],ver,fd,lvl,phatlvl);

	vp:=subsop(funcnum=NULL,vp);
end do;

avel1:=avel1/numrandfuncs;
avel2:=avel2/numrandfuncs;
aveinteg:=aveinteg/numrandfuncs;

intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals",loop-nops(loweruniquelist),".txt"),APPEND);
fprintf(intdata,"average:L1=%20.20g\tL2=%20.20g\tintegral=%20.20g\n",avel1,avel2,aveinteg);
fprintf(intdata,"minl1=%20.20g\tmaxl1:=%20.20g\n",minl1,maxl1);
fclose(intdata);

fprintf(fd,"}}");
fclose(fd);


end do;