(*  File is available at
    http://www.us.edu.pl/~gluza/csectors

   In case of any use, please, refer to
   J.Gluza, K. Kajda, T. Riemann, V. Yundin, 
   arXiv:1006.4728 [hep-ph] and arXiv:1010.1667

*)


Print["CSectors by K.Kajda and V. Yundin ver:1.0\nlast modified 22 sep 2009"]

BeginPackage["CSectors`"]

DoSectors::usage=""
invariants::usage = "A list of invariants e.g. {p1^2->s}"
PR::usage = "Propagator of the form: PR[momenta_,mass_,power of propagator_].\n 1/(k^2-m^2)^n=PR[k,m,n]"
eps::usage = "Dimensional regularization parameter for loop integrals."

(*--------------------------------------*)
(*      DEFAULT PROGRAM OPTIONS         *)
(*     can be modified if needed        *)
(*--------------------------------------*)
Options[DoSectors]={
	SetStrategy->C,
	SourceName->cpp_sectors	,
	TempFileDelete->True,
	TextDisplay->True,
        OutputForm->InputForm,
	LogFile->True,
	ShowErrors->True,
	IterationsLow->5,
	IterationsHigh->15,
	CallsLow->10000,
	CallsHigh->100000,
	compiler->"g++",
	cppflags->"-I/usr/local/include -I/usr/include",
	libs->"-L/usr/local/lib -lsector_decomposition -lginac -lcln -lgmp -L/usr/lib -lgsl -lgslcblas -lm"
};
(*-------------------------------------*)

Begin["`Private`"]

Xtranslate={x[n_]:>ToExpression["x"<>ToString[n]]};
MapTake[x_,n_]:=Map[Take[#,n]&, x] /.List->Plus;

OptChange[fun_,new___Rule]:=
  Module[{i},
     Map[If[Length[i=Position[{new},First[#]->_,1,1]]===0, 
	    #,{new}[[Sequence @@ First[i]]]]&, Options[fun]] ]



DoSectors[numerator_,integral_,momenta_,options___Rule][min_,order_]:=
 Module[{FUarray,U,F,Uexp,Fexp,dotinvariants,AboutNumerator,external,
	 tensors,gammas,IntComp,result},
	
	MQ[pos_,ind_]:=Part[mtilde,pos].q /.index->ind;
	MT[i_,j_]:=Part[mtilde,i,j];

	option=OptChange[DoSectors,options];
	Nup=Apply[List,First[integral]] /.PR[_,_,n_]->n;
	Nv=Plus @@ Nup;
	L=Length[momenta];
	
	FUarray=FUPolynomials[First[integral],momenta,invariants];
	U=FUarray[[1]];
	F=FUarray[[2]];

 If[numerator=!={1},
	mtilde=FUarray[[3]];
	q=FUarray[[4]];
	dotinvariants=invariants /.p_^2->p.p /.pa_*pb_->pa.pb;

	AboutNumerator=NumeratorStructure[numerator,momenta];
	external=AboutNumerator[[3]];
	tensors=TensorStructure[ AboutNumerator[[2]],AboutNumerator[[4]] ];
	
	tensors=Contract[tensors /.Xtranslate,external]/.dotinvariants;
	tensors=Simplify[ tensors*Table[1/(-2)^(r/2),{r,0,rank,2}] ];
	
(*--- Fix of polynomial going to zero...!!! ---*)
	tensors=Flatten[ { ZeroPolyFIX[#] } ]& /@ tensors,

	tensors={1};
	rank=0
  ];
	
	If[TextDisplay/.option,
	 WriteString[$Output,"Using strategy ",SetStrategy/.option,
               "\nU & F polynomials:\n   U = ",OutputForm[U]/.option,"\n   F = ",OutputForm[F]/.option,"\n" ]];

	gammas=Table[Nv-(2-eps)*L-r/2,{r,0,rank,2}];
	Uexp={Nv-2*(L+1)-rank,L+1};
	Fexp=Table[{-Nv+2*L+r/2,-L},{r,0,rank,2}];

	IntComp={tensors,gammas,Uexp,Fexp};


	numerics=CppRun[U,F,IntComp,min,order,option];


 Return[numerics /.eps^n_:>0 /;n>order]]



ZeroPolyFIX[polynomial_]:=
   Module[{chkpoly=Head[polynomial],condition,zeropoly,pickedpoly,samesigns,result},

	If[chkpoly===Integer || chkpoly===Rational,
		result=polynomial,

		condition[poly_]:=Map[{#>0,#<1}&, Variables[poly]] /. List->Sequence;

		zeropoly=Length[ FindInstance[{polynomial==0, condition[polynomial]}, Variables[polynomial]] ];

		(*Print["polynomial = ",polynomial];
		Print["zeropoly = ",zeropoly];  *)

		Switch[zeropoly,
			0,
		{result=polynomial},
			1,
		{
			pickedpoly=List @@ Expand[polynomial];
				(*Print["pickedpolys = ",pickedpoly];*)

			samesigns=Sign[pickedpoly] /. Sign[_]->1;
				(*Print["samesigns = ",samesigns];*)

			negative=Plus @@ Pick[pickedpoly,samesigns,-1];
				(*Print["negative = ",negative];*)

			positive=Plus @@ Pick[pickedpoly,samesigns, 1];
				(*Print["positive = ",positive];*)

			result={positive,negative}//Simplify
		}]

	];

	(*Print["result = ",result];*)

    Return[result]]



NumeratorStructure[numerator_,momenta_]:=
   Module[{IntMomMatrix,SeqRul,gmnRul,NoDotNum,
	   indices,ext,external,internal,IntPositions},
	
	IntMomMatrix=Flatten[Table[momenta[[i]].momenta[[j]],{i,L},{j,L}]];
        SeqRul=IntMomMatrix /.ka_ .kb_->(ka.kb->Sequence[ka,kb]);
	gmnRul=IntMomMatrix /.ka_ .kb_->(ka[mu_]**kb[nu_]->g[mu*nu]);
	NoDotNum=numerator /.Times->Dot /.x_^_->x.x /.SeqRul;

	rank=Length[NoDotNum];
	indices=Table[ToExpression["mu"<>ToString[i]],{i,1,rank}];

	ext=Flatten[NoDotNum /.Dot[x_,y_]:>DeleteCases[List[x,y],Alternatives @@ momenta]];
	external=Apply[NonCommutativeMultiply,Transpose[{ext, indices}]] 
			/.{k_,mu_}->k[mu] //. gmnRul /.NonCommutativeMultiply->Times;
	internal=Cases[NoDotNum /. Dot->Sequence, Alternatives @@ momenta];
	IntPositions=Map[Position[momenta,#]&,internal] //Flatten;


   Return[{internal,IntPositions,external,indices}]] 



TensorStructure[IntPositions_,indices_]:=
   Module[{UnsortedUnion,Drop2,poss,rul1,rul2,rul3,result,final},
	
	UnsortedUnion[x_]:=Reap[Sow[1,x],_,#1&][[2]];	
	Drop2[x_]:=Part[IntPositions,ToExpression[StringDrop[ToString[x], 2]] ];

	poss=Select[Permutations[indices], Signature[#]==1&];
	rul1=A[x_]*P[y_]:>Map[A[Take[#,x]]*P[Take[#,-y]]&,poss];
	rul2=A[x_]:>(Times @@ Map[g[Times @@ #]&, Partition[x, 2]]);
	rul3=P[x_]:>(Times @@ Map[P,x]);

	result=Map[UnsortedUnion[#]&,Table[A[r]P[rank-r],{r,0,rank,2}] /.rul1 /.rul2 /.rul3];
	final=result /.P[mu_]:>MQ[Drop2[mu],mu] /.g[mu_*nu_]:>g[mu*nu]MT[Drop2[mu],Drop2[nu]];

   Return[ Map[(Plus @@ #)&,final] ]]



Contract[expr_,external_]:=
      Module[{rulA,rulB,rulC,rulD,rulE,result},
        
    rulA=g[mu1_*mu3_]*g[mu2_*mu3_]:>g[mu1*mu2];
    rulB=g[__]^2->4-2*eps;
    rulC=g[mu1_*mu2_]*p_[mu2_]->p[mu1];
    rulD=p1_[x_]*p2_[x_]:>p1.p2/;p1=!=g;
    rulE=p1_[mu1_]^2:>p1.p1/;p1=!=x;
         
    result = Expand[expr*external] 
         //.rulA //.rulB //.rulC //.rulD //.rulE;

  Return[result]]



CppRun[U_,F_,IntComp_,min_,order_,option_]:=
 Module[{Qnozero={},xdim=Length[Nup],GiNaCexpr,tensors,gammas,Uexp,Fexp,R,ExtraLoop,
         xsymbols1,xsymbols2,chardim,ccresult,RawData,result={},errors={}},

	GiNaCexpr[x_] := StringReplace[ToString[CForm[x]], "Power" -> "pow"];
	chardim=ToString[Length[Characters[ToString[SourceName/.option]]]+6+2];
	strategy="STRATEGY_"<>ToString[SetStrategy /.option];

	tensors=IntComp[[1]];
	gammas=IntComp[[2]];
	Uexp=IntComp[[3]];
	Fexp=IntComp[[4]];
	R=Length[tensors];

	ExtraLoop=Map[# /.{List[x___]:>Length[List[x]],Times[___]->1}&,tensors];

	If[tensors=!={1},
	Do[
	  If[TextDisplay/.option,WriteString[$Output,"   Q"<>ToString[i]<>ToString[exlp]<>" = ", OutputForm[Part[tensors,i,exlp]]/.option,"\n" ]],
	{i,1,Length[tensors]}, {exlp,1,Part[ExtraLoop,i]} ]];
	
	xsymbols1=Sort[Variables[U]];
        xsymbols2=Map[(# /.x_:>ToString[x]<>"(\""<>ToString[x]<>"\")")&,xsymbols1];
	xsymbols2=Characters[ToString[xsymbols2]]//Flatten;

	If[TextDisplay/.option,WriteString[$Output,"Generating c++ source..."]];

(*------------------*)
(* LOOP-1 BEGIN *)
Do[

	If[tensors=!={1},
	  If[tensors[[int,exlp]]===0, Qnozero=Append[Qnozero,{int,exlp}];Goto["Qzero"]] ];

	strm = OpenWrite[ToString[SourceName/.option]<>ToString[int]<>ToString[exlp]<>".cc"];	
	If[TextDisplay/.option,
		WriteString[$Output,"Int",ToString[int],ToString[exlp],"..."]];

	WriteString[strm,"/*----------------------------\n",
                         "  This file was automatically generated by\n",
			 "  CSsectors on: ",Date[][[3]],".",Date[][[2]],".",Date[][[1]],
			 "\n -----------------------------*/\n\n",
	                 "#include <iostream>\n",
      If[LogFile/.option,"#include <fstream>\n",""],
                         "#include <vector>\n",
                         "#include <ginac/ginac.h>\n",
                         "#include \"sector_decomposition/sector_decomposition.h\" \n\n",
                         "int main()\n",
			 "{\n",
		         "      using namespace GiNaC;\n",
			 "      using namespace sector_decomposition;\n",
      If[LogFile/.option,"      char file["<>chardim<>"]=\""<>ToString[SourceName/.option]<>ToString[int]<>ToString[exlp]<>".log\";\n",""],
      If[LogFile/.option,"      std::ofstream ToFile;\n      ToFile.open (file,std::ios::trunc);\n",""],
			 "      CHOICE_STRATEGY = ",strategy,";\n\n",
			 "      int xdim = ",xdim,";\n",
                         "      symbol eps(\"eps\");\n",
                         "      symbol ",
    Apply[Sequence,Table[xsymbols2[[i]],{i,2,Length[xsymbols2]-1}]],";\n\n",
			 "/*-------------------------------------\n",
			 "  U, F and Q polynomials\n",
		         " -------------------------------------*/\n",
			 "      ex U = "<>GiNaCexpr[U]<>";\n",
			 "      ex F = "<>GiNaCexpr[F]<>";\n",
	If[tensors=!={1},"      ex Q = "<>GiNaCexpr[ tensors[[int, exlp]] ]<>";\n\n","\n"],
	                 "\n",
			 "/*-------------------------------------\n",
			 "  Gamma prefactor\n",
		         " -------------------------------------*/\n",
                         "      ex gams = pow(-1,",Nv,")",
    Apply[Sequence,Map["/tgamma("<>ToString[#]<>")"&,Nup]],
			 "*series(tgamma(",CForm[ gammas[[int]] ],"),eps==0,",order-min+1,");\n",
			 "      ex fact = evalf(gams.subs(Euler==0));\n\n",
			 "      std::vector<size_t> arg_index(xdim);\n",
                         "      for (int i=0; i<xdim ; i++) arg_index[i]=i;\n",
                         "\n\n",
			 "/*-------------------------------------\n",
			 "  Feynman parameters and exponents\n",
		         " --------------------------------------*/\n",
                         "      std::vector<ex> parameters;\n",
    Apply[Sequence,Table["      parameters.push_back("<>ToString[ xsymbols1[[i]] ]<>");\n",{i,xdim}]],
			 "\n",
                         "      std::vector<exponent> nu_minus_1;\n",
    Apply[Sequence,Table["      nu_minus_1.push_back(exponent("<>ToString[ Nup[[i]]-1 ]<>",0));\n",{i,1,xdim}]],
			 "\n\n",	
			 "/*-------------------------------------\n",
			 "  Polynomials and exponents\n",
		         " --------------------------------------*/\n",
                         "      std::vector<ex> poly_list;\n",
                         "      poly_list.push_back(U);\n",
                         "      poly_list.push_back(F);\n",
	If[tensors=!={1},"      poly_list.push_back(Q);\n\n","\n"],
                         "      std::vector<exponent> c;\n",
	                 "      c.push_back(exponent(",CForm[ Uexp[[1]] ],",",CForm[ Uexp[[2]] ],"));\n",
                         "      c.push_back(exponent(",CForm[ Fexp[[int,1]] ],",",CForm[ Fexp[[int,2]] ],"));\n",
	If[tensors=!={1},"      c.push_back(exponent(1,0));\n\n","\n"],
			 "\n",
		         "/*-------------------------------------\n",
			 "  Integrating\n",
		         " --------------------------------------*/\n",
		         "      integrand my_integrand(arg_index, nu_minus_1, poly_list, c, fact);\n",
			 "      monte_carlo_parameters mc_parameters( ",
			 IterationsLow /.option," , ",IterationsHigh /.option," , ",
                         CallsLow /.option," , ",CallsHigh /.option," );\n",
                         "\n",
      If[LogFile/.option,"      ToFile<<\" >> \"<<file<<\" << \"<<std::endl;\n",""],
			 "      for (int order=",min,"; order<=",order,"; order++)\n",
                         "      {\n",
                         "              integration_data global_data(parameters, eps, order);\n",
                         "              monte_carlo_result res = \n",
			 "              	 do_sector_decomposition(global_data, my_integrand, mc_parameters,0);\n\n",
                         "              std::cout<<res.get_mean()<<\"*\"<<pow(eps,order)<<\" \"<<res.get_error()<<\"*\"<<pow(eps,order)<<std::endl;\n",
      If[LogFile/.option,"              ToFile<<\"Order \"<<pow(eps,order)<<\": \"<<res.get_mean()<<\" +/- \"<<res.get_error()<<std::endl;\n",""],
			 "      }\n\n",
      If[LogFile/.option,"      ToFile.close();\n",""],
                         "      return 0;\n",
	                 "}\n"];

	Close[strm];

   Label["Qzero"],

{int,1,R}, {exlp,1,Part[ExtraLoop,int]}];
(* LOOP-1 END *)
(*------------------*)

	If[TextDisplay/.option,WriteString[$Output,"done\nCompiling source code..."]];

(*------------------*)
(* LOOP-2 BEGIN *)
Do[	
  If[FreeQ[Qnozero,{int,exlp}],
	If[TextDisplay/.option,WriteString[$Output,"Int",ToString[int],ToString[exlp],"..."]];
	Run[ToString[compiler/.option]<>" -o "<>ToString[SourceName/.option]<>ToString[int]<>ToString[exlp]<>" "<>ToString[cppflags/.option]<>" "<>ToString[SourceName/.option]<>ToString[int]<>ToString[exlp]<>".cc "<>ToString[libs/.option]]
  ],

{int,1,R}, {exlp,1,Part[ExtraLoop,int]}];
(* LOOP-2 END *)
(*------------------*)


  If[TextDisplay/.option,WriteString[$Output,"done\nRunning binary file....."]];


(*------------------*)
(* LOOP-3 BEGIN *)
Do[	
  If[FreeQ[Qnozero,{int,exlp}],
        If[ TextDisplay/.option,WriteString[$Output,"Int",ToString[int],ToString[exlp],"..."] ];
	RawData=Import["!./"<>ToString[SourceName/.option]<>ToString[int]<>ToString[exlp],"Table"];
	ccresult=ToExpression[StringReplace[#,{"e+"->"*10^","e-"->"*10^-"}]& /@ RawData];

	result=Join[result,{MapTake[ccresult, 1]}];
	errors=Join[errors,{MapTake[ccresult,-1]}/.Plus->List]
   ],	
{int,1,R}, {exlp,1,Part[ExtraLoop,int]}];
(* LOOP-3 END *)
(*------------------*)

	errors=Flatten[{Sqrt[Apply[Plus,Flatten[errors]^2]/.Plus->List]//PowerExpand}];
	If[TextDisplay/.option,WriteString[$Output,"done\n"]];
	If[TempFileDelete/.option,Run["rm -f ",ToString[SourceName/.option]<>"*"]];


 Return[ If[ShowErrors/.option,{Apply[Plus,result],errors},Apply[Plus,result]] ]]




FUPolynomials[integral_,momenta_,kinematics_] :=
 Module[{denom = List @@ integral, int, m, j, q, U, mtilde,inv},

  int = Sum[x[i]*denom[[i]] /.PR[k_,m_,_]->k^2-m^2,{i,Length[denom]}];
(*Print["int: ",int];*)

  m = Table[If[i != j, 1/2, 1]*
    Coefficient[int, momenta[[i]]*momenta[[j]]], {i,L}, {j,L}];
(*Print["m: ",m];*)

  j = int - momenta.m.momenta;
  q = Table[Expand[-1/2*Coefficient[j, momenta[[i]]]], {i,L}];
  j = Expand[j + 2*momenta.q] /. kinematics;
  U = Det[m];
  F = Simplify[U*(Expand[q.Inverse[m].q] - j) /. kinematics];
  invm=Inverse[m];	
  mtilde=Expand[U*invm] //Simplify;
(*Print["mtilde: ",mtilde];*)


 Return[{Simplify[U /.Xtranslate], F /.Xtranslate, mtilde, q /.a_.*p_*x[n_]:>a*p[index]*x[n]}]]


End[]

EndPackage[]

