#--------------------------------------------------------------------------------------------------------
# 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
##
#########################################################################################################


SimpleMultipleZeroRefiner:=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;


	#initialization
	nsys:=Vector(sys):
	nzero:=Vector(mzero):
	n:=nops(var):

	
	#coeffecient normalisation
	st:=time();
	for i from 1 to n do
		nsys[i]:=sys[i]/max(map(abs,{coeffs(sys[i])}));
	od:
	Jac:=VectorCalculus[Jacobian](nsys,var);
	Ev:=[seq(var[i]=nzero[i],i=1..n)]:
	userinfo(2,SimpleMultipleZeroRefiner,print(`session time for coeffecient normalisation:`,time()-st));
	userinfo(2,SimpleMultipleZeroRefiner,print(`coeffecient normalisation for original system:`,nsys));
	userinfo(2,SimpleMultipleZeroRefiner,print(`system evaluation at beginning:`,eval(nsys,Ev)));	


	#normalisation transformation #1
	st:=time();
 	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))]);
	userinfo(2,SimpleMultipleZeroRefiner,print(`session time for #1 normalisation transformation:`,time()-st));
	userinfo(2,SimpleMultipleZeroRefiner,print(`U and W for #1 normalisation transformation:`,UTrans,WTrans));


	#N1 iteration
	st:=time();
	Jac_nzero:=HermitianTranspose(UTrans).Jac_nzero.WTrans:	#use chain rule to avoid linear transformation 
	nzero:=nzero-SubMatrix(WTrans,1..-1,2..n).MatrixInverse(SubMatrix(Jac_nzero,1..n-1,2..n)).SubMatrix(HermitianTranspose(UTrans),1..n-1,1..-1).convert(eval(nsys,Ev),Vector):	#N1
	Ev:=[seq(var[i]=nzero[i],i=1..n)]:
	userinfo(2,SimpleMultipleZeroRefiner,print(`session time for N1:`,time()-st));
	userinfo(3,SimpleMultipleZeroRefiner,print(`Jacobian for #1 normalisation transformation:`,Jac_nzero));
	userinfo(1,SimpleMultipleZeroRefiner,print(`refined simple multiple zero for N1:`,nzero));
	userinfo(2,SimpleMultipleZeroRefiner,print(`system evaluation for N1:`,eval(nsys,Ev)));


	#normalisation transformation #2
	st:=time();
	Jac_nzero:=HermitianTranspose(UTrans).eval(Jac,[seq(var[i]=nzero[i],i=1..n)]).WTrans:	#updated Jacobian after N1
	VVTrans,UUTrans:=LinearAlgebra[SingularValues](Jac_nzero,output=['Vt','U']);
	WWTrans:=Matrix([HermitianTranspose(SubMatrix(VVTrans,-1,1..-1)),HermitianTranspose(SubMatrix(VVTrans,1..-2,1..-1))]);
	userinfo(2,SimpleMultipleZeroRefiner,print(`session time for #2 normalisation transformation:`,time()-st));
	userinfo(3,SimpleMultipleZeroRefiner,print(`Jacobian for N1:`,Jac_nzero));
	userinfo(2,SimpleMultipleZeroRefiner,print(`U and W for #2 normalisation transformation:`,UUTrans,WWTrans));


	#N2 itertation
	st:=time();
	Jac_nzero:=HermitianTranspose(UUTrans).Jac_nzero.WWTrans:	#use chain rule to avoid linear transformation
	userinfo(3,SimpleMultipleZeroRefiner,print(`Jacobian for #2 normalisation transformation:`,Jac_nzero));
	total_U:=HermitianTranspose(UTrans.UUTrans);	#linear transformation for equations
	total_W:=WTrans.WWTrans;	#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])
	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 >>:
 	od:
	dualsupp:=dualsupp/mul;
	nzero:=nzero-convert(total_U.subs(Ev,pol_last),Vector)[n]/convert(total_U.subs(Ev,dualsupp),Vector)[n]/mul*convert(SubMatrix(total_W,1..-1,1),Vector);	#N2
	userinfo(2,SimpleMultipleZeroRefiner,print(`session time for N2:`,time()-st));
	userinfo(1,SimpleMultipleZeroRefiner,print(`refined simple multiple zero for N2:`,nzero));
	userinfo(3,SimpleMultipleZeroRefiner,print(`system evaluation for N2:`,eval(nsys,[seq(var[i]=nzero[i],i=1..n)])));
	
	RETURN(nzero);

end proc: