#--------------------------------------------------------------------------------------------------------
# author : Nan Li
# date   : Feb. 28 2017
#--------------------------------------------------------------------------------------------------------
	
#####   Refine Simple Multiple Zeros   ##################################################################
##   
##	Calling sequence:  SimpleMultipleZeroRefiner(sys,var,mzero,mul)
##      
##	Input:  a polynomial system, a variable order, an approximate simple multiple zero and its multiplicity
##
##	Output:  the refined simple multiple zero
##
#########################################################################################################


DeltaOutput:=proc(sys,var,mzero,mul)

	local nsys,nzero,n,Jac,Ev,Jac_nzero,UTrans,VTrans,WTrans,UUTrans,VVTrans,WWTrans,total_U,total_W,inv_Jac,pol_last,pol_list,para_last,para_list,dualsupp,temp,i,j,k,st,SeqDelta;


	#initialization
	nsys:=sys:
	nzero:=mzero:
	n:=nops(var):
	SeqDelta:=[]:


	Jac:=VectorCalculus[Jacobian](nsys,var);
	Ev:=[seq(var[i]=nzero[i],i=1..n)]:
	
	
	Jac_nzero:=eval(Jac,Ev):
	SeqDelta:=[op(SeqDelta),Jac_nzero[n,1]]:
	inv_Jac:=Matrix(map(x->1/x,LinearAlgebra[Diagonal](SubMatrix(Jac_nzero,1..-2,2..-1))),shape=diagonal);
	trn_Jac:=SubMatrix(Jac,1..-1,2..-1);
	pol_last:=expand(SubMatrix(Jac,1..-1,1));	#polynomial after apply first order differential operator
	dualsupp:=diff(convert(pol_last,Vector),var[1]);	#equavilent to map(diff,pol_last,var[1])
	SeqDelta:=[op(SeqDelta),convert(subs(Ev,dualsupp),Vector)[n]/2]:
	pol_list:=<< pol_last >>:
	para_list:=Matrix(0,n-1);
	for k from 2 to mul-1 do	
  		para_last:=-inv_Jac.SubVector(eval(dualsupp,Ev),1..-2):	#calculate parameters
  		pol_last:=( dualsupp + trn_Jac . para_last )/k:	#polynomial after apply k-th order differential operator
  		para_list:=<< Transpose(para_last),para_list >>:
		dualsupp:=diff(convert(pol_last,Vector),var[1]);
  		temp:=pol_list . para_list;
  		for i from 1 to n-1 do                    
  			dualsupp:=dualsupp+diff(convert(SubMatrix(temp,1..-1,i),Vector),var[i+1]);	#equavilent to map(diff,temp[i],var[i+1])
  		od:
		pol_list:=<< pol_list | pol_last >>:
		SeqDelta:=[op(SeqDelta),convert(subs(Ev,dualsupp),Vector)[n]/(k+1)]:
 	od:
	
	RETURN(SeqDelta);


# 	Jac_nzero:=eval(Jac,Ev):
#	VTrans,UTrans:=LinearAlgebra[SingularValues](Jac_nzero,output=['Vt','U']);
#	WTrans:=Matrix([HermitianTranspose(SubMatrix(VTrans,-1,1..-1)),HermitianTranspose(SubMatrix(VTrans,1..-2,1..-1))]);
#
#
#	Jac_nzero:=HermitianTranspose(UTrans).Jac_nzero.WTrans:	#use chain rule to avoid linear transformation
	SeqDelta:=[op(SeqDelta),Jac_nzero[n,1]]:
#	total_U:=HermitianTranspose(UTrans);	#linear transformation for equations
#	total_W:=WTrans;	#linear transformation for variables
#	inv_Jac:=Matrix(map(x->1/x,LinearAlgebra[Diagonal](SubMatrix(Jac_nzero,1..-2,2..-1))),shape=diagonal);
#	pol_last:=expand(Jac.SubMatrix(total_W,1..-1,1));	#polynomial after apply first order differential operator
#	dualsupp:=VectorCalculus[Jacobian](convert(pol_last,Vector),var).SubMatrix(total_W,1..-1,1);	#equavilent to map(diff,pol_last,var[1])
#	SeqDelta:=[op(SeqDelta),convert(total_U.subs(Ev,dualsupp),Vector)[n]/2]:
#	pol_list:=<< pol_last >>:
#	para_list:=Matrix(0,n-1);
#	for k from 2 to mul-1 do	
#  		para_last:=-inv_Jac.SubMatrix(total_U.eval(dualsupp,Ev),1..-2,1..-1):	#calculate parameters
#  		pol_last:=( dualsupp + Jac . (SubMatrix(total_W,1..-1,2..-1).para_last) )/k:	#polynomial after apply k-th order differential operator
#  		para_list:=<< Transpose(para_last),para_list >>:
#		dualsupp:=VectorCalculus[Jacobian](convert(pol_last,Vector),var).SubMatrix(total_W,1..-1,1);
#  		temp:=pol_list . para_list;
#  		for i from 1 to n-1 do                    
#  			dualsupp:=dualsupp+VectorCalculus[Jacobian](convert(SubMatrix(temp,1..-1,i),Vector),var).SubMatrix(total_W,1..-1,i+1);	#equavilent to map(diff,temp[i],var[i+1])
#  		od:
#		pol_list:=<< pol_list | pol_last >>:
#		SeqDelta:=[op(SeqDelta),convert(total_U.subs(Ev,dualsupp),Vector)[n]/(k+1)]:
# 	od:
#	
#	RETURN(SeqDelta);

end proc: