findFuncMaxDifferenceWithoutZero:=proc(T,lvl,funcnum,epsilon)
        local ver,vertex,ma,mi,didit;
        #ma:=T[[0],funcnum];
        #mi:=T[[0],funcnum];
        didit:=0;
        ver:=vertices(lvl);
        for vertex in ver do
                if abs(T[vertex,funcnum])>epsilon then
                        if didit=1 then
                                if T[vertex,funcnum]<mi then mi:=T[vertex,funcnum];
                                end if;
                                if T[vertex,funcnum]>ma then ma:=T[vertex,funcnum];
                                end if;
                        else
                                ma:=T[vertex,funcnum];
                                mi:=T[vertex,funcnum];
                                didit:=1;
                        end if;
                end if;
        end do;
        if didit=1 then
                RETURN(ma-mi);
        else
                RETURN(0);
        end if;
end; 

LocFindFuncMaxDifferenceWithoutZero:=proc(T,lvl,ver,funcnum,epsilon)
        local vertex,ma,mi,didit;
        #ma:=T[[0],funcnum];
        #mi:=T[[0],funcnum];
        didit:=0;
        for vertex in ver do
                if abs(T[vertex,funcnum])>epsilon then
                        if didit=1 then
                                if T[vertex,funcnum]<mi then mi:=T[vertex,funcnum];
                                end if;
                                if T[vertex,funcnum]>ma then ma:=T[vertex,funcnum];
                                end if;
                        else
                                ma:=T[vertex,funcnum];
                                mi:=T[vertex,funcnum];
                                didit:=1;
                        end if;
                end if;
        end do;
        if didit=1 then
                RETURN(ma-mi);
        else
                RETURN(0);
        end if;
end;

rp:=proc(L,m)
 option remember;
 local L1, i,L2;
   if m = 0 then RETURN([]);
   else
   if m = 1 then
	L1:=[];
	for i in L do
		L1:=[op(L1),[i]];
	end do; 
	RETURN(L1);
   else
     L1 := rp(L,m-1);
     L2:=[];
     #[op(map(add_one,L1,0)),op(map(add_one,L1,1)),op(map(add_one,L1,2))]
     for i in L do
	L2:=[op(L2),op(map(add_one,L1,i))];
     end do;
     RETURN(L2);
   fi;
   end if;
 end;


numericFunc:=proc(T,lvl,func)
	local v, vers;
	vers:=vertices(lvl);	
	for v in vers do
		T[v,func]:=evalf(T[v,func]);
	end do;
end;

LocNumericFunc:=proc(T,lvl,vers,func)
        local v;
        for v in vers do
                T[v,func]:=evalf(T[v,func]);
        end do; 
end;

funcQuotient:=proc(T,lvl,infunc1,infunc2,outfunc)
	local vers, vertex;
	vers:=vertices(lvl);
	for vertex in vers do
		if(abs(T[vertex,infunc2])<.000001) then#i didn't have a error 0 i just had 0 before
			T[vertex,outfunc]:=0;
			if(abs(T[vertex,infunc1])>.000001) then
				#print("the bottom was 0 but not the top"); 
			end if;
		else T[vertex,outfunc]:=T[vertex,infunc1]/T[vertex,infunc2];
		end if;
	end do;
end;


LocFuncQuotient:=proc(T,lvl,vers,infunc1,infunc2,outfunc)
        local vertex;
        for vertex in vers do
                if(abs(T[vertex,infunc2])<.000001) then#i didn't have a error 0 i just had 0 before
                        T[vertex,outfunc]:=0;
                        if(abs(T[vertex,infunc1])>.000001) then
                                #print("the bottom was 0 but not the top");
                        end if;
                else T[vertex,outfunc]:=T[vertex,infunc1]/T[vertex,infunc2];
                end if;
        end do;
end;


multiplyFuncByConstant:=proc(T,lvl,fnum,c)
	local vers, vertex;
	vers:=vertices(lvl);
	for vertex in vers do
		T[vertex,fnum]:=T[vertex,fnum]*c;
	end do;
end;

LocMultiplyFuncByConstant:=proc(T,lvl,vers,fnum,c)
        local vertex;
        for vertex in vers do
                T[vertex,fnum]:=T[vertex,fnum]*c;
        end do;
end;

funcDifference:=proc(T,lvl,infunc1,infunc2,outfunc)
	local vers, vertex;
	vers:=vertices(lvl);
	for vertex in vers do
		T[vertex,outfunc]:=T[vertex,infunc1]-T[vertex,infunc2];
	end do;
end;

LocFuncDifference:=proc(T,lvl,vers,infunc1,infunc2,outfunc)
        local vertex;
        for vertex in vers do
                T[vertex,outfunc]:=T[vertex,infunc1]-T[vertex,infunc2];
        end do;
end;

funcSumWithConstant:=proc(T,lvl,c,infunc1,infunc2,outfunc)
	local vers, vertex;
	vers:=vertices(lvl);
	for vertex in vers do
		T[vertex,outfunc]:=T[vertex,infunc1]+c*T[vertex,infunc2];
	end do;
end;

LocFuncSumWithConstant:=proc(T,lvl,vers,c,infunc1,infunc2,outfunc)
        local vertex;
        for vertex in vers do
                T[vertex,outfunc]:=T[vertex,infunc1]+c*T[vertex,infunc2];
        end do;
end;

funcProduct:=proc(T,lvl,infunc1,infunc2,outfunc)
	local vers,vertex;
	vers:=vertices(lvl);
	for vertex in vers do
		T[vertex,outfunc]:=T[vertex,infunc1]*T[vertex,infunc2];
	end do;
end;

LocFuncProduct:=proc(T,lvl,vers,infunc1,infunc2,outfunc)
        local vertex;   
        for vertex in vers do
                T[vertex,outfunc]:=T[vertex,infunc1]*T[vertex,infunc2];
        end do;
end;

funcAbs:=proc(T,lvl,infunc,outfunc)
        local vers,vertex;
        vers:=vertices(lvl);
        for vertex in vers do
                T[vertex,outfunc]:=abs(T[vertex,infunc]);
        end do;
end;

LocFuncAbs:=proc(T,lvl,vers,infunc,outfunc)
        local vertex;
        for vertex in vers do
                T[vertex,outfunc]:=abs(T[vertex,infunc]);
        end do;
end;

LocFuncRotate:=proc(T,lvl,ver,list,infunc,outfunc)
	local vertex,nvertex,i,temp;
	for vertex in ver do
		nvertex:=vertex;
		for i from 1 to nops(vertex) do
			nvertex[i]:=list[vertex[i]+1];
		end do;
		#print(cat("change ",vertex," to ",nvertex));
		if nops(nvertex)>1 then 
			if nvertex[nops(vertex)-1]>nvertex[nops(vertex)] then 
				temp:=nvertex[nops(nvertex)-1];
				nvertex[nops(nvertex)-1]:=nvertex[nops(nvertex)];
				nvertex[nops(nvertex)]:=temp;
			end if;
		end if;
		T[nvertex,outfunc]:=T[vertex,infunc];
	end do;
end;

funcRotate:=proc(T,lvl,list,infunc,outfunc)
        local ver,vertex,nvertex,i,temp;
        ver:=vertices(lvl);
        for vertex in ver do
                nvertex:=vertex;
                for i from 1 to nops(vertex) do
                        nvertex[i]:=list[vertex[i]+1];
                end do;
                #print(cat("change ",vertex," to ",nvertex));
                if nops(nvertex)>1 then
                        if nvertex[nops(vertex)-1]>nvertex[nops(vertex)] then
                                temp:=nvertex[nops(nvertex)-1];
                                nvertex[nops(nvertex)-1]:=nvertex[nops(nvertex)];
                                nvertex[nops(nvertex)]:=temp;
                        end if;
                end if;
                T[nvertex,outfunc]:=T[vertex,infunc];
        end do;
end;

discLaplac:=proc(T,lvl,infunc,outfunc)
	local ver,vertex,i;
	ver:=vertices(lvl);
	make_neighbors(T,lvl,ver);
	ver:=subsop(1=NULL,2=NULL,3=NULL,ver);
	for vertex in ver do
		T[vertex,outfunc]:=-4*T[vertex,infunc];
		for i from 1 to 4 do
			T[vertex,outfunc]:=T[vertex,outfunc]+T[T[vertex,i],infunc];
		end do;
	end do;
	T[[1],outfunc]:=0;
	T[[2],outfunc]:=0;
	T[[0],outfunc]:=0;
end;



discInnerProduct:=proc(T,lvl,infunc1,infunc2)
	local ver,vertex,ret;
	ret:=0;
	ver:=vertices(lvl);
	ver:=subsop(1=NULL,2=NULL,3=NULL,ver);
	for vertex in ver do
		ret:=ret+T[vertex,infunc1]*T[vertex,infunc2];
	end do;
	RETURN (ret);
end;

findFuncMax:=proc(T,lvl,funcnum)
	local ver,vertex,m;
	m:=T[[0],funcnum];
	ver:=vertices(lvl);
	for vertex in ver do
		if T[vertex,funcnum]>m then m:=T[vertex,funcnum];
		end if;
	end do;
	RETURN(m);
end;

LocFindFuncMax:=proc(T,lvl,ver,funcnum)
        local vertex,m;
        m:=T[ver[1],funcnum];
        for vertex in ver do
                if T[vertex,funcnum]>m then m:=T[vertex,funcnum];
                end if;
        end do;
        RETURN(m);
end;    

findFuncMin:=proc(T,lvl,funcnum)
	local ver,vertex,m;
	m:=T[[0],funcnum];
	ver:=vertices(lvl);
	for vertex in ver do
		if T[vertex,funcnum]<m then m:=T[vertex,funcnum];
		end if;
	end do;
	RETURN(m);
end;

LocFindFuncMin:=proc(T,lvl,ver,funcnum)
        local vertex,m;
        m:=T[ver[1],funcnum];
        for vertex in ver do
                if T[vertex,funcnum]<m then m:=T[vertex,funcnum];
                end if;
        end do;
        RETURN(m);
end;

findFuncMax2:=proc(T,lvl,funcnum)
        local ver,vertex,m;  
        m:=[T[[0],funcnum],[0]];
        ver:=vertices(lvl);
        for vertex in ver do
                if T[vertex,funcnum]>m[1] then 
			m[1]:=T[vertex,funcnum];
                	m[2]:=vertex;
		end if;
        end do;
        RETURN(m);
end;
                

findFuncMin2:=proc(T,lvl,funcnum)
        local ver,vertex,m;
        m:=[T[[0],funcnum],0];
        ver:=vertices(lvl);
        for vertex in ver do
                if T[vertex,funcnum]<m[1] then 
			m[1]:=T[vertex,funcnum];
                	m[2]:=vertex;
		end if;
        end do;
        RETURN(m);
end;

LocFindFuncMax2:=proc(T,lvl,ver,funcnum)
        local vertex,m;  
        m:=[T[ver[1],funcnum],ver[1]];
        for vertex in ver do
                if T[vertex,funcnum]>m[1] then 
			m[1]:=T[vertex,funcnum];
                	m[2]:=vertex;
		end if;
        end do;
        RETURN(m);
end;
                

LocFindFuncMin2:=proc(T,lvl,ver,funcnum)
        local vertex,m;
        m:=[T[ver[1],funcnum],ver[1]];
        for vertex in ver do
                if T[vertex,funcnum]<m[1] then 
			m[1]:=T[vertex,funcnum];
                	m[2]:=vertex;
		end if;
        end do;
        RETURN(m);
end;



normalizeFunc:=proc(T,lvl,funcnum)
	local m;
	m:=discInnerProduct(T,lvl,funcnum,funcnum);
	m:=sqrt(m);
	multiplyFuncByConstant(T,lvl,funcnum,1/m);
end;

compareFuncList:=proc(T,lvl,funcnum,valuelist)
	local i, v, vers, length, minvar, maxvar;
	
	length:=nops(valuelist);
	vers:=ourVertices(lvl);	
	maxvar:=0;
	
	#print("vers= "); print(vers);	

	for v in vers do
		minvar:=abs(valuelist[1] - T[v,funcnum]);

		for i from 2 to length do
			if (abs(valuelist[i] - T[v,funcnum]) < minvar) then
				minvar:=abs(valuelist[i] - T[v,funcnum]);
			end if;
		end do;

		if (minvar > maxvar) then maxvar:=minvar; end if;
	end do;

	return(maxvar);
end;

LocCompareFuncList:=proc(T,lvl,vers,funcnum,valuelist)
        local i, v, length, minvar, maxvar;
        
        length:=nops(valuelist);
        maxvar:=0;

        #print("vers= "); print(vers);
        
        for v in vers do   
                minvar:=abs(valuelist[1] - T[v,funcnum]);
                
                for i from 2 to length do
                        if (abs(valuelist[i] - T[v,funcnum]) < minvar) then
                                minvar:=abs(valuelist[i] - T[v,funcnum]);
                        end if;
                end do;

                if (minvar > maxvar) then maxvar:=minvar; end if;
        end do;

        return(maxvar);
end;

removeDuplicates:=proc(l)
	local i,j,bad;
	bad:=[];
	for i from 1 to nops(l) do
		for j from i+1 to nops(l) do
			if(l[i]=l[j]) then
				bad:=[op(bad),j=NULL];
			end if;
		end do;
	end do;
	RETURN(subsop(op(bad),l));
end;


copyRotateScaleChunk:=proc(T,lvl,phatlvl,infunc, outfunc, vertex)
	local l, badlist, sforbidden, tforbidden, i, j, p,  v, S, tversmn, tovers, tnvers, sversmn, sovers, snvers, tlvl, slvl;

	S:=gasket(2+phatlvl-lvl);

	tlvl:=lvl;
	slvl:=(tlvl - lvl + 2);
	
	tversmn:=ourVertices(tlvl);
        sversmn:=ourVertices(slvl);
        make_neighbors(T, tlvl,tversmn);
        make_neighbors(S, slvl,sversmn);

	tforbidden:=[T[vertex,1],T[vertex,2],T[vertex,3],T[vertex,4]];
	sforbidden:=[S[[1,2],1],S[[1,2],2],S[[1,2],3],S[[1,2],4]];
	
	tovers:=[op(tforbidden),vertex];
        sovers:=[op(sforbidden),[1,2]]; 
        tnvers:=[];
        snvers:=[];
	#print("tovers=");
	#print(tovers);
	#print("tforbidden=");
	#print(tforbidden);
	#print("sforbidden=");
	#print(sforbidden);
	#print("sovers=");
	#print(sovers);

	#print("mark 1");

	for tlvl from (lvl +1) to phatlvl do
	        #print("tovers.25=");
        	#print(tovers);
		slvl:=(tlvl - lvl + 2);

		tversmn:=vertices(tlvl);
		sversmn:=vertices(slvl);
		make_neighbors(T, tlvl,tversmn);
		make_neighbors(S, slvl,sversmn);

                #print("tovers.35=");
                #print(tovers);

		
		for j from 5 to nops(tovers) do
			tnvers:=[op(tnvers),tovers[j]];
			snvers:=[op(snvers),sovers[j]];
			#print(cat("tovers ",j," = ",tovers[j]));
			#print(cat("sovers ",j," = ",sovers[j]));
			

			for i from 1 to 4 do
				#print(cat("tovers ",j," neighbor ",i,"= ",T[tovers[j],i]));
	                        #print(cat("sovers ",j," neighbor ",i,"= ",S[sovers[j],i]));
				#if (sameCell(lvl,T[tovers[j],i],vertex)=true) then
					#print("tovers is in cell");
					tnvers:=[op(tnvers),T[tovers[j],i]];
					snvers:=[op(snvers),S[sovers[j],i]];
				#end if;
			end do;
		end do;

		#print("mark 2");
       		#print("tovers.5=");
        	#print(tovers);
		
		if (tlvl=lvl + 1) then 
			l:=nops(tnvers);
			for j from 1 to l do
                		for i from 1 to 4 do
                        		tnvers:=[op(tnvers),T[tnvers[j],i]];
                                	snvers:=[op(snvers),S[snvers[j],i]];
                        	end do;
                	end do;
		end if;

		#print("mark 3");

		badlist:=[];
		for j from 1 to nops(tnvers) do		
			#print("mark 3.1");
			if (tnvers[j]=tforbidden[1] or tnvers[j]=tforbidden[2] or tnvers[j]=tforbidden[3] or tnvers[j]=tforbidden[4]) then
				badlist:=[op(badlist),j=NULL];
			end if;
		end do;
		
		tnvers:=subsop(op(badlist),tnvers);
		snvers:=subsop(op(badlist),snvers);

		#for i from 1 to nops(tnvers) do
		#	if tnvers[i]=[0,0,0,0,0,1] then
		#		print("was a tnver at lvl");
		#		print(tlvl);
		#		print("while phatlvl=");
		#		print(phatlvl);
		#	end if;
		#end do;

		tovers:=[op(tforbidden),op(tnvers)];
		sovers:=[op(sforbidden),op(snvers)];
		tnvers:=[];
		snvers:=[];

		#if tlvl=lvl+1 then
		#	print("after lvl1, tovers=");
		#	print(tovers);
		#end if;
		#tovers:=removeDuplicates(tovers);
		#sovers:=removeDuplicates(sovers);
		#print("tlvl=");
		#print(tlvl);
		#print(cat("slvl= ",slvl));
		#print(sovers);

		#for i in vertices(slvl) do
		#	if(sameCell(2,[1,2],i)) then
		#		p:=0;
		#		for j from 1 to nops(sovers) do
		#			if (i=sovers[j]) then p:=p+1; end if;
		#		end do;
		#		if (p = 0) then 
		#			print(cat("Vertex ",i," not in sovers at level",slvl));
		#		end if;
		#	end if;
		#end do;

	end do;

        #print("tovers2=");
        #print(tovers);
	#for i from 1 to nops(tovers) do
	#	if tovers[i]=[0,0,0,0,0,0,1] then
	#		print("foudn it");
	#	end if;
	#end do;
	

	#print("mark 4");

	clearFunc(T,max(phatlvl,phatlvl-lvl+2),outfunc);

	for j from 1 to nops(tovers) do
		T[sovers[j],outfunc]:=T[tovers[j],infunc];
	end do;

	#clearFunc(T, phatlvl, 9);
	#clearFunc(T, phatlvl, 10);

	#for j from 1 to nops(tovers) do
        #        T[tovers[j],9]:=1;
	#	T[sovers[j],10]:=1;
        #end do;
	
	#print("mark 5");
	
	#for j from 1 to nops(tovers) do
	#	if (not sameCell(lvl,tovers[j],vertex)) then
	#		print(cat("vertex ",tovers[j], " out of cell in T\n"));
	#		fprintf(stderr,"out of Cell in T");
	#	end if;
	#	if (not sameCell(2,sovers[j],[1,2])) then
	#		print(cat("vertex ",sovers[j], "out of cell in S\n"));
	#		fprintf(stderr,"out of Cell in S");
	#	end if;
	#end do;

	return(nops(tovers));
end;


sameCell:=proc(lvl,vertex1,vertex2)
	local i, v1rep1, v1rep2, v2rep1, v2rep2, lv1, lv2, v1, v2;	

	if(lvl < 1) then return(true); end if;
	lv1:=nops(vertex1);
	lv2:=nops(vertex2);

	if (lv1=1) then
		v1:=[op(vertex1), op(vertex1)];
		lv1:=2;
	else 	v1:=vertex1;
	end if;

	if (lv2=1) then     
                v2:=[op(vertex2), op(vertex2)];
                lv2:=2;
        else    v2:=vertex2;
        end if;
	
	if (lv1 < lvl + 2) then
		v1rep1:=[op(subsop(lv1=NULL, (lv1 -1)=NULL, v1)), v1[lv1-1], seq(v1[lv1],   i=1..(lvl + 2 - lv1))]; 
		v1rep2:=[op(subsop(lv1=NULL, (lv1 -1)=NULL, v1)), v1[lv1],   seq(v1[lv1 -1],i=1..(lvl + 2 - lv1))];
	else	v1rep1:=v1; v1rep2:=v1; 
	end if;

	if (lv2 < lvl +2) then
		v2rep1:=[op(subsop(lv2=NULL, (lv2-1)=NULL, v2)), v2[lv2-1], seq(v2[lv2],    i=1..(lvl + 2 - lv2))];
		v2rep2:=[op(subsop(lv2=NULL, (lv2-1)=NULL, v2)), v2[lv2],   seq(v2[lv2 -1], i=1..(lvl + 2 - lv2))];
	else	v2rep1:=v2; v2rep2:=v2;
	end if;

	if 	([seq(v1rep1[i],i=1..lvl)]=[seq(v2rep1[i],i=1..lvl)]) then return(true);
	elif	([seq(v1rep1[i],i=1..lvl)]=[seq(v2rep2[i],i=1..lvl)]) then return(true);
	elif    ([seq(v1rep2[i],i=1..lvl)]=[seq(v2rep1[i],i=1..lvl)]) then return(true);
	elif    ([seq(v1rep2[i],i=1..lvl)]=[seq(v2rep2[i],i=1..lvl)]) then return(true);
	else					                           return(false);
	end if;
end;


clearPartFunc:=proc(T,lvl,funcnumber)
	local v, ver;
	
	ver:=ourVertices(lvl);

	for v in ver do
		if(not sameCell(2,v,[1,2])) then
			T[v,funcnumber]:=0;
		end if;
	end do;
end;

LocClearFunc:=proc(T,vers,funcnumber)
	local v;
	for v in vers do
		T[v,funcnumber]:=0;
	end do;
end;

copyFunc:=proc(T,lvl,infunc,outfunc)
	local vert,v;
	vert:=ourVertices(lvl);
	for v in vert do
		T[v,outfunc]:=T[v,infunc];
	end do;
end;