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;

lvl:=4;
phatlvl:=8;
BMTOL:=.0000000000000001;
numrandfuncs:=50;
path:=cat("../data2/level",lvl,"/");

T:=gasket(phatlvl):

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

template:=fopen("../data/template.rtf",READ);
fd:=fopen(cat(path,phatlvl,"phat",lvl,"dirdata.rtf"),WRITE);
intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals.txt"),WRITE);
fprintf(intdata,"l1, l2, and integral data for the %d phat %d extra dirichelt\n",phatlvl,lvl);
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:=1000000000000000000000;
maxl1:=0;

for i from 1 to numrandfuncs do
	clearFunc(T,phatlvl,13);
	fd2:=fopen(cat(path,lvl,"lambdalist.txt"), READ);
	funcnum:=(rand() mod nops(vp))+1;
	for j from 1 to nops(vl) do
		ulambda:=fscanf(fd2,"%g\n")[1];
		readFunc(T,12,cat(path,phatlvl,"phat",lvl,"eigen.",j));
		co:=T[vp[funcnum],12]*(3^(lvl+1))/(2*calcbm(ulambda,[],BMTOL));
		if (abs(co) > 0) then
			funcSumWithConstant(T,phatlvl,co,13,12,13);
		end if;
	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);
	minl1:=min(minl1,l1);
	maxl1:=max(maxl1,l1);
        integ:=simp_rule(T,phatlvl,13);
	intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals.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;

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

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

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

intdata:=fopen(cat(path,phatlvl,"phat",lvl,"extradirintegrals.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);