read "itall.txt":

#set these things to one if you want to do them, 0 if you don't. tolerances are for detecting problems.
#always set PSIMAKINGTOLERANCE to 0. It could save time otherwise, but is just a bad idea.

MAKEBASEEIGENS:=0:
MAKEALLEIGENS:=1:
SAVENORMS:=0:
NORMALIZEEIGENS:=0:
TESTSTUFF:=0:
	CHECKORTHONORMAL:=1:
	TOLNORMALITY:=.0000000001:
	TOLORTHOGONALITY:=.0000000001:
	CHECKEIGENVALS:=1:
	TOLEIGENVALS:=.0000000001:
FATTENEIGENS:=0:
MAKEPSIS:=0:
	PSIMAKINGTOLERANCE:=0;
CHECKPHATEIGENVALS:=0:
	
#these paths should probably all be set to the same thing. Originally they were meant for different things, but it didn't work out.

lvl:=4;
phatlvl:=5;
path:=cat("../data/level",lvl,"/");
inpath:=cat("../data/level",lvl,"/");
outpath:=cat("../data/level",lvl,"/");
numEigens:=(3^(lvl+1)-3)/2;

T:=gasket(phatlvl):


#MAKEBASEEIGENS: makes your base eigenfunctions you will need. (note: level 1 eigen functions are special. you will need to
#	copy them from our data (they were made by hand) to the path you are using. 

if(MAKEBASEEIGENS=1) then
	print(outpath);
	makeOrthoEigen5(T,lvl,outpath);
	makeOrthoEigen6(T,lvl,outpath);
	print("All done making basic orthonormal Eigen Set");
end if:

#MAKEALLEIGENS: takes your base eigenfunctions, and spectral decimates them to give you a full basis of (discrete) eigenfunctions
#	at the level you're working at. 	

if(MAKEALLEIGENS=1) then
	genalllevel(T,lvl,path);
	print("All done making full orthogonal Eigen Set");
end if:

#SAVENORMS: saves the norms of your eigenfunctions (before you normalize them) in case you want to know them.

if (SAVENORMS=1) then
	fd:=fopen(cat(path,lvl,"squarednorms.txt"),WRITE);

	for i from 1 to numEigens do
        	readFunc(T,12,cat(path,lvl,"eigen.",i));
        	fprintf(fd,"%20.20g\n",discInnerProduct(T,lvl,12,12));
	end do;

	fclose(fd);
	print("All done Saving Norms");
end if:

#NORMALIZEEIGENS: normalize all your eigenfunctions. (So your basis is orthoNORMAL)

if (NORMALIZEEIGENS=1) then
	for i from 1 to numEigens do
       		readFunc(T,12,cat(path,lvl,"eigen.",i));
	        normalizeFunc(T,lvl,12);
        	saveFunc(T,12,lvl,cat(path,lvl,"eigen.",i));
	end do;
	print("All done Normalizing");
end if:

#TESTSTUFF: Tests various things about the eigenbasis to make sure they're okay. You can figure out the details.

if (TESTSTUFF=1) then
	startEigen:=12;
	tol2:=TOLEIGENVALS;
	tol1:=TOLORTHOGONALITY;
	tol0:=TOLNORMALITY;

	for i from 1 to numEigens do
	        readFunc(T,i+startEigen-1,cat(path,lvl,"eigen.",i));
	end do;

	if (CHECKORTHONORMAL=1) then
		for i from 1 to numEigens do
		        for j from i to numEigens do
        		        m:=discInnerProduct(T,lvl,i+startEigen-1,j+startEigen-1);
                        
		                if (i=j) then
		                        if(abs(m-1) > tol0) then
		                        printf("(%d,%d)=%20.20g\n",i,j,m);
        	        	        end if;
        		        else    if(abs(m) > tol1) then                
                	        	printf("(%d,%d)=%20.20g\n",i,j,m);
                		        end if;
        		        end if;
	
       			 end do;
		end do;
		print("All done Checking Orthonormality");	
	end if;

	if(CHECKEIGENVALS=1) then

		fd:=fopen(cat(path,lvl,"lambdalist.txt"), READ);
		for i from 1 to numEigens do
        		discLaplac(T,lvl,i+startEigen-1,startEigen+numEigens);
       			funcQuotient(T,lvl,startEigen+numEigens,i+startEigen-1,startEigen+numEigens);
        		lambda:=fscanf(fd,"%g\n")[1];
        		lambda:=-1*lambda;
       			m:=compareFuncList(T,lvl,startEigen+numEigens,[0,lambda]);
        		if (abs(m) > tol2) then 
                		printf("max variance of (laplacian(%d)/%d) =%20.20g\n",i,i,m);
        		end if;
		end do;

		fclose(fd);

		print("All done Checking Eigens");

	end if;

end if:


#FATTENEIGENS: fills in your eigenfunctions to the phatlvl. This lets you see what the eigenfunctions (of the actual (not discrete))
#	SG Laplacian look like. It chooses -1's for epsilons (except the first one for lambda=6's which must be +1) to get the lowest
#	eigenvalues. 


if (FATTENEIGENS=1) then

	fd:=fopen(cat(inpath,lvl,"lambdalist.txt"), READ);
	df:=fopen(cat(outpath,phatlvl,"phat",lvl,"lambdalist.txt"), WRITE);
	
	for i from 1 to numEigens do
	        lambda:=fscanf(fd,"%g\n")[1];
	        readFunc(T,12,cat(inpath,lvl,"eigen.",i));
	        if (lambda=6) then 
        	        if (phatlvl > lvl) then
        	                fixedeps:=[1,op(makeEpsilons(phatlvl - lvl -1))];
        	                newAlgoWrapper(T,lvl,12,lambda,fixedeps);
        	        	newlambda:=makeLambda(lambda,fixedeps);
			end if;
       		else    newAlgoWrapper(T,lvl,12,lambda,makeEpsilons(phatlvl - lvl));
			newlambda:=makeLambda(lambda,makeEpsilons(phatlvl - lvl));
        	end if;
	        fprintf(df,"%20.20g\n",newlambda);
		saveFunc(T,12,phatlvl,cat(outpath,phatlvl,"phat",lvl,"eigen.",i));
	end do;
	
	fclose(df);
	fclose(fd);
	print("All done Fattening Eigens");	
	
end if:

#MAKEPSIS: makes the Psi functions from the Eigenfunctions. PSIMAKINGTOLERANCE should be set to 0 to avoid error.

if (MAKEPSIS=1) then
	tol:=PSIMAKINGTOLERANCE;

	vl:=ourVertices(lvl):
	vl:=subsop(1=NULL,2=NULL,3=NULL,vl):
	vp:=ourVertices(phatlvl):

	for i from 1 to nops(vl) do
	        clearFunc(T,phatlvl,13);
	        for j from 1 to nops(vl) do

			readFunc(T,12,cat(inpath,phatlvl,"phat",lvl,"eigen.",j));
	                co:=T[vl[i],12];
	                if (abs(co) > tol) then
        	                funcSumWithConstant(T,phatlvl,co,13,12,13)
        	        end if;
       		end do;
        
		saveFunc(T,13,phatlvl,cat(outpath,phatlvl,"phat",lvl,"psi.",i));
	end do;

	print("All done Making Psis");
end if:

#CHECKPHATEIGENCALS: Another diagnostic for your eigenfunctions. Checks the Laplacian down to the phatlvl to see that it agrees with
#	what the list says it should be. (the list is written out by FATTENEIGENS).

if(CHECKPHATEIGENVALS=1) then

	df:=fopen(cat(path,phatlvl,"phat",lvl,"lambdaerrs.txt"),WRITE);
	fd:=fopen(cat(path,phatlvl,"phat",lvl,"lambdalist.txt"), READ);
	
	for i from 1 to numEigens do
	   
		readFunc(T,startEigen,cat(path,phatlvl,"phat",lvl,"eigen.",i));
	
	       	discLaplac(T,phatlvl,startEigen,startEigen+1);
	        funcQuotient(T,phatlvl,startEigen+1,startEigen,startEigen+2);
	        lambda:=fscanf(fd,"%g\n")[1];
	        lambda:=-1*lambda;
	        m:=compareFuncList(T,phatlvl,startEigen+2,[0,lambda]);
	        fprintf(df,"Laplacian error for function %d =\t%20.20g\n",i,m);
	                
	end do;

	fclose(fd);
	fclose(df);

	print("All done Checking Phat Eigens");

end if;


print("Have a nice Day");