#tries again to minimize the determinants of the Psi matrix for irregular sampling.
#Vertices vary within the cell around where they should be (excluding points of intersection of these cells. 
#Search is performed at increasing depths (i.e. phatlvls) to see how the minimum varies. 


read "itall.txt":

with(LinearAlgebra);

mindet:=proc(T,baselvl,lvl,funclist,listverts,tol)
	local l, m, n, d, oldlow, cycles, lowdet,i,j,k, newlist, tmpverts, nbs;
	

	#make_neighbors(T,lvl,T[lvl,-1]);
	l:=nops(listverts);

	nbs:=[seq([],i=1..l)];

	for i from 1 to l do
		for vertex in T[lvl,-1] do
			if(sameCell((baselvl+2),listverts[i],vertex)) then
			nbs[i]:=[op(nbs[i]),vertex];
			end if;
		end do;
	end do;
	
	newverts:=listverts;
	m:=Matrix([seq([seq(T[newverts[i],funclist[j]],j=1..l)],i=1..l)]);
	lowdet:=abs(Determinant(m));
	
	cycles:=0;
	oldlow:=lowdet+10*tol;
	
	while (lowdet < oldlow - tol) do	
		oldlow:=lowdet;
		for i from 1 to l do
			for vertex in nbs[i] do
				tmpverts:=subsop(i=vertex,newverts);
				m:=Matrix([seq([seq(T[tmpverts[n],funclist[k]],k=1..l)],n=1..l)]);	
				d:=abs(Determinant(m));
				if (d<lowdet) then
					newverts:=tmpverts;
					lowdet:=d;
				end if;
			end do;
		end do;
		cycles:=cycles+1;
	end do;

	return([lowdet,newverts]);

end;

PRINTVERTS:=true:
PRINTFIRSTLEVEL:=true:

lvl:=3;
phatlvl:=8;
T:=verGasket(phatlvl);
tol:=.000001;
startfunc:=12;

path:=cat("../data/level",lvl,"/");
lverts:=subsop(1=NULL,2=NULL,3=NULL,T[lvl,-1]);

for i from 1 to nops(lverts)  do
	readFunc(T,startfunc+i-1,cat(path,phatlvl,"phat",lvl,"psi.",i));
end do;

fd:=fopen(cat(path,lvl,"mindetssti.txt"),WRITE);
fprintf(fd,"level\t(base level %d)\tdeterminant\n",lvl);

if (PRINTFIRSTLEVEL) then
	st:=lvl+1:
else	st:=lvl+2:
end if;

for level from st to phatlvl do
	d:=mindet(T,lvl,level,[seq(12+i-1,i=1..nops(lverts))],lverts,tol);
	print(cat("mindet at level ",level," is ",d[1]));
	if (PRINTVERTS) then print(cat("at vertices ",d[2])); end if;
	fprintf(fd,"%d\t\t%20.20g\n",level,d[1]);
	if (PRINTVERTS) then fprintf(fd,"%a\n",d[2]); end if;
end do;
fclose(fd);