#--------------------------------------------------------------------------------------------------------
# author : Nan Li
# date   : Mar. 2nd 2010
#--------------------------------------------------------------------------------------------------------

#####   Approximative Nullspace  ########################################################################
##
##	Calling sequence:  AppNull(A,tol)
##      
##	Input:  One matrix A with finite elements, and a torlerance
##
##	Output:  The approximative nullspace of a matrix at a given torlerance
##
#########################################################################################################


AppNull:=proc(A,tol)
 
	local S1,V1,i,C;

	S1, V1 :=LinearAlgebra[SingularValues](A, output=['S', 'Vt']);
	for i from 1 to op(1,S1) do
		if S1[i]<tol then
			C:=SubMatrix(V1,i..RowDimension(V1),1..ColumnDimension(V1));
			break;
		fi;
		if i=op(1,S1) then
			C:=SubMatrix(V1,i+1..RowDimension(V1),1..ColumnDimension(V1));
			break;
		fi;
	od;

	RETURN(C);

end proc;


#####   Differential Function acts on a Polynomial System at Zero  ######################################
##   
##	Calling sequence:  DifOnPol(lsys,ord,dualsupport)
##      
##	Input:  An polynomial system, an order and a dualsupport
##
##	Output:  The polynomial after differential operator acts on the polynomial system
##
#########################################################################################################


DifOnPol:=proc(lsys,ord,dualsupport)

	local L,v,dualsupp,coefficient,i,dif,j,k,u;

	v:=Vector(nops(lsys));
	dualsupp:=dualsupport+1/ord[1];
	if nops(dualsupp)>2 then
		for i from 1 to nops(dualsupport) do
			L:=[];
			coefficient:=1;
			dif:=Conv(op(dualsupport)[i],ord);
			for j from 1 to nops(ord) do
				for k from 1 to dif[j] do
					L:=[op(L),ord[j]];
				od;
				coefficient:=coefficient*dif[j]!;
			od;
			for u from 1 to nops(lsys) do
				v[u]:=v[u]+coeffs(op(dualsupport)[i],ord)*diff(lsys[u],L)/coefficient;
			od;
		od;
	else
		L:=[];
		coefficient:=1;
		dif:=Conv(dualsupport,ord);
		for j from 1 to nops(ord) do
			for k from 1 to dif[j] do
				L:=[op(L),ord[j]];
			od;
			coefficient:=coefficient*dif[j]!;
		od;
		for u from 1 to nops(lsys) do
			v[u]:=v[u]+coeffs(dualsupport,ord)*diff(lsys[u],L)/coefficient;
		od;
	fi;

	RETURN(v);

end proc:
	
		
#####   Evalue Vector   #################################################################################
##   
##	Calling sequence:  EvalVector(ord)
##      
##	Input:  An order
##
##	Output:  EvalVector at the zero
##
#########################################################################################################


EvalVector:=proc(ord)

	local L,i;

	L:=[];
	for i from 1 to nops(ord) do
		L:=[op(L),ord[i]=0];
	od;

	RETURN(L);

end proc:
	
		
#####   Differential Operator   #########################################################################
##   
##	Calling sequence:  DifOperator(dualsupp,ord,i)
##      
##	Input:  A dualsupport, an order and ith
##
##	Output:  The differential functional
##
#########################################################################################################


DifOperator:=proc(dualsupp,ord,i)


	RETURN(expand(dualsupp*ord[i]));

end proc:
	
		
#####   Convert to Differential Functional   ############################################################
##   
##	Calling sequence:  Conv(mon,ord)
##      
##	Input:  A monomial and an order
##
##	Output:  The Differential Functional respect to the monomial
##
#########################################################################################################


Conv:=proc(mon,ord)

	local L,i;

	L:=[];
	for i from 1 to nops(ord) do
		L:=[op(L),degree(mon,ord[i])];
	od;

	RETURN(L);

end proc:
	
		
#####   Transform A Polynomial System Based an A Matrix   ###############################################
##   
##	Calling sequence:  TransPol(lsys,TransMatrix,ord)
##      
##	Input:  A polynomial system, a matrix and an order
##
##	Output:  The transformed system
##
#########################################################################################################


TransPol:=proc(lsys,TransMatrix,ord)

	local i,L,v;

	L:=[];
	v:=TransMatrix.convert(ord,Vector);
	for i from 1 to nops(ord) do
		L:=[op(L),ord[i]=v[i]];
	od;

	RETURN(expand(subs(L,lsys)));

end proc:
	
		
#####   Form A Full-Rank Matrix   #######################################################################
##   
##	Calling sequence:  FullRank(nullvector,tol)
##      
##	Input:  A vector and a tolerance
##
##	Output:  A full-rank matrix based on the vector 
##
#########################################################################################################


FullRank:=proc(nullvector,tol)

	local A;

	A:=AppNull(nullvector,tol);

	RETURN(Matrix([Transpose(nullvector),Transpose(A)]));
	
end proc:
	
		
#####   Separate the Differential Functional into Part   ################################################
##   
##	Calling sequence:  Part(rest,ord,j)
##      
##	Input:  A differential functional, an order and jth
##
##	Output:  The part of the differential functional
##
#########################################################################################################


Part:=proc(rest,ord,j)

	local L,k,r;

	L:=0;
	r:=rest+1/ord[1];
	if nops(r)>2 then
		for k from 1 to nops(rest) do
			if degree(op(rest)[k],ord[j])=0 then
				L:=L+op(rest)[k];
			fi;
		od;
	else
		if degree(rest,ord[j])=0 then
			L:=L+rest;
		fi;
	fi;

	RETURN(expand(L));

end proc:
	
		
#####   Compute the Coefficient Before A Differential Operator   ########################################
##   
##	Calling sequence:  MultiFactorial(ind)
##      
##	Input:  An index
##
##	Output:  The coefficient
##
#########################################################################################################


MultiFactorial:=proc(ind)

	local mf,i;

	mf:=1;
	for i from 1 to nops(ind) do
		mf:=mf*(ind[i]!);
	od;

	RETURN(mf);
	
end proc:
	
		
#####   Transform the Dual Basis Back to Original System   ##############################################
##   
##	Calling sequence:  TransBack(dualbasis,TransMatrix,ord)
##      
##	Input:  A dualbasis, a matrix and an order
##
##	Output:  The dualbasis respect to the original system
##
#########################################################################################################


TransBack:=proc(dualbasis,TransMatrix,ord)

	local i,j,k,DualBasis,poly,tmid,mon,mid,ttmid;

	DualBasis:=[];
	for i from 1 to nops(dualbasis) do
		poly:=0;
		tmid:=dualbasis[i]+1/ord[1];
		if nops(tmid)>2 then
			for j from 1 to nops(dualbasis[i]) do
				mon:=0;
				mid:=TransPol(op(dualbasis[i])[j],Transpose(TransMatrix),ord);
				ttmid:=mid+1/ord[1];
				if nops(ttmid)>2 then
					for k from 1 to nops(mid) do
						mon:=mon+MultiFactorial(Conv(op(mid)[k],ord))*op(mid)[k];
					od;
				else
					mon:=mon+MultiFactorial(Conv(mid,ord))*mid;
				fi;
				poly:=poly+mon/MultiFactorial(Conv(op(dualbasis[i])[j],ord));
			od;			
		else
			mon:=0;
			mid:=TransPol(dualbasis[i],Transpose(TransMatrix),ord);
			ttmid:=mid+1/ord[1];
			if nops(ttmid)>2 then
				for k from 1 to nops(mid) do
					mon:=mon+MultiFactorial(Conv(op(mid)[k],ord))*op(mid)[k];
				od;
			else
				mon:=mon+MultiFactorial(Conv(mid,ord))*mid;
			fi;
			poly:=poly+mon/MultiFactorial(Conv(dualbasis[i],ord));
		fi;
		DualBasis:=[op(DualBasis),poly];
	od;

	RETURN(DualBasis);
	
end proc:
	
		
#####   Refine the Mutiple Root   #######################################################################
##   
##	Calling sequence:  RootRefine(lsys,ord,appx,tol)
##      
##	Input:  One polynomial system, an order, an approxiate root and a tolerance
##
##	Output:  The multiplity, the index and the refined root 
##
#########################################################################################################


MultiStructure:=proc(pol,ord,mroot,tol)

	local lsys,F,s,t,L,Jac,nullvector,TransMatrix,dualbasis,DualBasis,
	adddual,depth,multiplicity,involutive,dualsupp,DualSupp,DualMatrix,
	DualSupport,G,i,j,k,rest,NewRow,dif,coe,NewMatrix;

	lsys:=pol;
	s:=nops(ord);
	t:=nops(lsys);
	L:=[];
	for i from 1 to s do
		L:=[op(L),ord[i]=mroot[i]];
	od;
	G:=EvalVector(ord);
	Jac:=eval(VectorCalculus[Jacobian](lsys,ord),L);
	nullvector:=AppNull(Jac,tol);
	TransMatrix:=FullRank(nullvector,tol);
	Jac:=Jac.TransMatrix;
	dualbasis:=[1,ord[1]];
	DualMatrix:=convert(ord,Matrix);
	depth:=1;
	multiplicity:=2;
	involutive:=0;
	while involutive=0 do
		NewRow:=[];
		dualsupp:=DifOperator(dualbasis[-1],ord,1);
		NewRow:=[op(NewRow),dualsupp];
		rest:=Part(dualbasis[-1],ord,1);
		for j from 2 to s do
			NewRow:=[op(NewRow),DifOperator(rest,ord,j)];
			rest:=Part(rest,ord,j);
		od;
		DualMatrix:=Matrix([[DualMatrix],[convert(NewRow,Matrix)]]);
		for i from multiplicity by -1 to 2 do
			F:=[];
			for k from 1 to i-1 do
				F:=[op(F),ord[1]];
			od;
			dif:=diff(dualsupp,F)/(i-1)!;
			coe:=coeff(dualbasis[i],ord[1]^(i-1));
			for j from 2 to s do
				dualsupp:=dualsupp+DualMatrix[i,j]*eval(coeff(dif,ord[j]),G)/coe;
			od;
		od;
		DualSupp:=[dualsupp,ord[2]];
		for i from 3 to s do
			DualSupp:=[op(DualSupp),ord[i]];
		od;
		NewMatrix:=Matrix([convert(eval(DifOnPol(lsys,ord,TransBack([dualsupp],TransMatrix,ord)[1]),L),Matrix),SubMatrix(Jac,1..-1,2..-1)]);
		nullvector:=AppNull(NewMatrix,tol);
		if RowDimension(nullvector)=0 then
			involutive:=1;
		fi;
		if RowDimension(nullvector)>0 then
			adddual:=nullvector.convert(DualSupp,Vector);
			dualbasis:=[op(dualbasis),adddual[1]];			
			depth:=depth+1;
			multiplicity:=multiplicity+1;
		fi;
	od;
	DualBasis:=TransBack(dualbasis,TransMatrix,ord);
	userinfo(2,MultiStructure,print(`The dual basis of new polynomial at isolated zero is:`,dualbasis));
	userinfo(1,MultiStructure,print(`The multiplicity is:`,multiplicity));
	userinfo(2,MultiStructure,print(`The depth is:`,depth));
	userinfo(1,MultiStructure,print(`The dual basis of original polynomial at isolated zero is:`,DualBasis));
	for i from 1 to multiplicity do
		userinfo(1,MultiStructure,print(`The dual basis `,i,` acts on original polynomial at isolated zero is:`,map(expand,eval(DifOnPol(pol,ord,DualBasis[i]),L))));
	od;

end proc: