(* AMBRE = Automatic Mellin Barnes REpresentation *)              
Print["by K.Kajda   ver: 1.1"]
Print["last modified 27 May 2007"]



BeginPackage["AMBRE`"]


Fauto::usage = "Fauto[0] allows manual modification of F-polynomial. Fauto[1] acitvates automatic mode, set by default"
Fullintegral::usage = "Fullintegral[{numerator},{propagators},{internal momenta}] is a function for input integrals"
invariants::usage = "A list of invariants e.g. {p1^2->s}"
IntPart::usage = "IntPart[iteration] prepares subintegral for given internal momentum."
numerator::usage = "Numerator for a given subloop. Currently it is {1}, no numerators. Please check for future upgrades."
integral::usage = "Propagators for a given subloop"
momentum::usage = "Momentum for a given subloop"
SubLoop::usage = "Subloop[integral]"
BarnesLemma::usage = "BarnesLemma[repr_,1] tries to apply 1st Barnes lemma on a given representation. \n BarnesLemma[repr_,2] is used to apply 2nd Barnes lemma.\n\n These two functions automatically apply lemmas on every M-B integral variable, only if a given variable appears in Gamma functions but not in functions related to kinematcal variables."
PR::usage = "PR[momenta,mass,power of propagator] is a propagator. Mass must be a symbol."
fupc::usage = "F-polynomial variable"
FX::usage = "FX[X[1]+X[2]+...]^2 means (X[1]+X[2]+...)^2"
X::usage = "Feynman parameters of a form X[1],X[2],..."
eps::usage = "Epsilon variable"
z::usage = "M-B integral variable, writen in a form z1,z2,..."
ARint::usage = "M-B integral for tensor diagram. ARint as a function is able to display given integral result: ARint[1,repr] - displays first integral in a result."


Fmanual::mode = "U and F polynomials will be calculated in MANUAL mode. Now you can modify F polynomial (fupc).";
Fauto::mode = "U and F polynomials will be calculated in AUTO mode. In order to use MANUAL mode execute Fauto[0].";


(***********************************************)
(*               Options                       *)
(***********************************************)
Option[SubLoop]={Result->False};
Option[BarnesLemma]={Shifts->False};

Begin["`Private`"]


(***********************************************)
(*               Functions                     *)
(***********************************************)
SFblue[x_]:=StyleForm[x, FontColor -> RGBColor[0, 0, 1]];
SFred[x_]:=StyleForm[x, FontColor -> RGBColor[1, 0, 0]];
ToExprStr[x_,y_]:=ToExpression[ToString[x]<>ToString[y]];
UnsUnApp[x_,y_]:=Reap[Sow[1,Append[x,y]],_,#1&][[2]];
ToTake[x_]:=If[Length[x]>iteration,Take[x,iteration],x];
ToList[x_]:=If[Head[x]===Plus,List @@ x,{x}];
ARint[num_,repr_]:=repr[[num]] /.ARint[x_]->x;


(***********************************************)
(*                 Fauto                       *)
(***********************************************)
Fauto[n_] :=
  Block[{},
    Fmode = n;
    If[Fmode == 0, 
    {
      Message[Fmanual::mode];
      fupc=FUpolyfun[integral,momentum,invariants];
      Print["fupc=",fupc]
    }];
    If[Fmode == 1, Message[Fauto::mode]]
 ]


(***********************************************)
(*              Fullintegral                   *)
(***********************************************)
Fullintegral[ARnumerators_,ARintegral_,ARmomentas_]:=
                (FullInt={ARnumerators,ARintegral,ARmomentas});


(***********************************************)
(*             Part of integral                *)
(***********************************************)
IntPart[iter_] :=
  Module[{mom,prop,rul1,rul2,den,a},

	 mom[j_]:= Last[FullInt][[j]];
	 prop[k_]:= 
	 Times @@ Complement[a[[k]],Flatten[Table[a[[i]],{i,1,k-1}]],{}];

	 Fmode=1;
	 iteration=iter;
         If[iter==1,
	 {
	    znumb={{1,0}};
	    PropList=FullInt[[2]];
	    ToAddList={1}
	 }];
	 
	 rul1 = PR[x_,y_,a_]PR[x_, y_, b_]->PR[x,y,a+b];
	 rul2 = PR[s1__,a_,n1_]PR[s2__,a_,n2_]:>
	         PR[s1,a,n1+n2]/;Expand[(s1)^2-(s2)^2]==0;
	 Propagators=Times @@ DeleteCases[
			Take[PropList,iteration],_Integer] //.{rul1,rul2};

         a = 
	 Table[
	       If[Head[Propagators]===PR,
		   {Propagators},
		   Cases[Propagators,PR[a_.*mom[i]+x_.,__]]],
	     {i,Length[Last[FullInt]]}];
      
    
	 integral=prop[iteration];
	 momentum=mom[iteration];
	 numerator=Cases[FullInt[[1]],a_.*momentum];
	 If[Length[numerator]==0,numerator=1];

	 Print["numerator=",numerator];
	 Print["integral=", integral];
	 Print["momentum=", momentum];

	 Fauto[Fmode]

 ]


(***********************************************)
(*        Main for given subloop               *)
(***********************************************)
  SubLoop[integral_,option___Rule] := 
  Module[{options,subres,sub0,NewProp,NewPropAux,toadd,
	   intover,Subpart,result},	 
	 
    
     options=If[Length[{option}]==0,Option[SubLoop],{option}];
     Print["Iteration nr",iteration,": >>Integrating over ",momentum,"<<"];


 (* Used for tadpoles *)

     If[(integral /. PR[__]->PR /.PR^n_.->n)==1 ,
      {
	{mass,lambda} = integral/. PR[_,b_,c_]->{b,c};
	subres = 
	   (-1)^lambda*Gamma[lambda+eps-2]/
	   (Gamma[lambda]*(mass^2)^(lambda+eps-2));
        znumb=UnsUnApp[znumb,{iteration+1,znumb[[iteration,2]]}];
	Goto[tadpole]
      },{}];

(********************)

     If[Fmode==1,
     {
       Print["Computing U & F polynomial in AUTO mode >>Fauto[",Fmode,"]<<"];
       fupc = Expand[FUpolyfun[integral,momentum,invariants]//Simplify];
       Print[fupc]
     },
     {
       Print["U & F polynomial was computed by user >>Fauto[",Fmode,"]<<"]
     }];


     subres = ARintegration[integral, fupc];    

Label[tadpole];
 
     PropList = ToTake[PropList];
     ToAddList = ToTake[ToAddList];

     sub0 = subres /.PR[__]->1;

     rul = PR[x_,y_,a_]->PR[x,y,-a];
     NewProp = PowerExpand[subres/sub0]/. PR[x_,y_]^a_->PR[x,y,a]; 
  
       
     If[NewProp==1,{NewPropAux=NewProp*iteration},{},{NewPropAux=NewProp}];


     PropList = UnsUnApp[PropList,NewPropAux /.rul];     
     ToAddList = UnsUnApp[ToAddList,sub0];

     toadd = Times @@ Take[ToAddList,iteration];

     If[iteration==Length[FullInt[[3]]],
     {
       Print["Final representation:"];
       result = toadd*subres;
     
       If[Head[result]=!=List,
	  Print[SFblue[result]/. ARint->SFred[ARint]];
	  Goto["ToResult"]],

       If[Head[result]===List && (Result /. options),
	  Print[SFblue[result]/. ARint->SFred[ARint]],
	  Print[SFblue[ARint[SFred[#]]]& /@ 
		Flatten[Position[result,ARint[x_]]]]]
     },
     {      
       intover = Take[Last[FullInt],iteration];
       intover = #<>##2& @@ (ToString[#]<>"..."& /@ intover);
       Print["Representation after integrating over: ",intover];
       Subpart = ToExprStr[SubLoop,iteration];
       result = Subpart[SFblue[toadd*sub0],SFred[NewProp]]
     }];

Label["ToResult"];    
       
 Return[result]]


(*************************************************)
(*      Evaluates given integral                 *)
(*************************************************)
ARintegration[integral_,fupc_] :=
  Module[{ipow, Nup,
          Spart,xint,Xpart,Gpart,
          result},


    ipow = List @@ integral /. PR[__, nu_]->nu;
    Nup = Total[ipow]-2+eps;

    Spart = Numk[integral, momentum, numerator, Nup];
    xint = xintegral[Spart,integral,ipow];

    Xpart = Spart*xint /.X[_]->1;
    Gpart = Times @@ (Gamma[#]& /@ ipow);
    result = Expand[(-1)^(Nup+2-eps)/Gpart*Xpart];

    If[Length[xint]==1,
    {
      result = Plus @@ Flatten[result];
    },
    {
      result = ARint[# /.List->Plus]& /@ result
    }];


 Return[result]]


(***********************************************)
(*            Mellin-Barnes formula            *)
(***********************************************)
MBarnes[fupc_, integral_, Nup_]:= 
  Module[{imin,fuplist,poly=fupc,nfu,
	  MB1,pow=Nup,MB2,barnes,
	  FXcheck,aux,temp={}},
	   
	
    znumb = Take[znumb,iteration];
    imin = znumb[[iteration,2]];
    
    Label[FXstart];
    fuplist = ToList[Expand[poly]];
 
    nfu = Length[fuplist];
    MB1 = Sum[ToExprStr[z,i+imin],{i,1,nfu-1}]+pow;
    MB2 = Product[fuplist[[i]]^ToExprStr[z,i+imin]*
		  Gamma[-ToExprStr[z,i+imin]],{i,1,nfu-1}];
    barnes = MB2*fuplist[[nfu]]^(-MB1)*Gamma[MB1]/Gamma[pow];
 
    znumb = UnsUnApp[znumb,{iteration+1,imin+nfu-1}];


    FXcheck = Cases[PowerExpand[barnes*aux],FX[_]^_];
    
    If[Length[FXcheck] != 0,
       {
	 {poly,pow} = FXcheck /.{FX[x_]^a_}->{x,-a};
	 temp = Append[temp,barnes/. FX[__]->1];
	 imin = znumb[[iteration+1,2]];
	 znumb = Delete[znumb,-1];
         Goto[FXstart]
       }];

    If[Length[temp]!=0,barnes=First[temp]*barnes];


 Return[barnes]]


(***********************************************)
(*            Integration over  x              *)
(***********************************************)
xintegral[fnu_, integral_, ipow_] := 
  Module[{tomult,expr,xonly,numG,denG},
     
    tomult = Product[X[i]^ipow[[i]], {i,1,Length[ipow]}];
    expr = PowerExpand[tomult*fnu/.Gamma[_]->1];
    xonly = expr/(expr /.X[_]->1);
   
    numG = xonly /. X[i_]^n_.:>Gamma[Expand[n]] //Expand;
    denG = numG //. Gamma[a_]*Gamma[b_]->Gamma[a+b];
   
       
 Return[numG/denG]]


(***********************************************)
(*         Vector & Tensor parts               *)
(***********************************************)
Numk[integral_, momentum_, numerator_, Nup_] :=
  Module[{mvar , Fpre, Flst, FF, AP, Qin, ext, 
	  indmu, murul,
	   APobject, gmn, Pm, APsum1, Spart},
             

     FF[i_] := MBarnes[fupc,integral,Nup-i];
     mvar = Length[numerator]; 

	  
     AP[x_,y_,ext_]:=
      Module[{A,result},

       If[x!=0,
	  A=Times @@ (gmn[#]& /@ Partition[Take[ext,x],2]),
	  A=1];
       result =
        A*Product[Pm[Take[ext,{i}]],{i,x+1,Length[ext]}];
       result = 
	result /. gmn[{a_,b_}]->gmn[a*b] 
	       /.{gmn[a_],gmn[b_]}->gmn[a]*gmn[b]
	       /. Pm[{a_}]->Pm[a];

      Return[result]];



     If[mvar == 0, 
     {
       Spart = {Gamma[Nup]*FF[0]};
       Goto[NumkEND];
     }];


     Qin = ToPoly[integral,momentum][[3]];
     ext = numerator/momentum;
	  indmu = Table[ToExpression["mu"<>ToString[i]],{i,1,mvar}];

     APobject = 
       Flatten[
        Table[
	  Gamma[Nup-r/2]/(-2)^(r/2)*
	  Fpre[r/2]*    
	  AP[r,mvar-r,#],{r,0,mvar,2}]& /@ Permutations[indmu]
         ]// Union;
  

     Flst = Flatten[Cases[#,Fpre[_]]& /@ APobject] //Union;

     APobject = (Cases[APobject,#*_.]& /@ Flst) /. Fpre->FF;
     murul = #[[1]]->#[[2]]& /@ (indmu*ext);
   
     APobject = APobject /. murul;
     APsum1 = APobject /.Pm[x_]:>Pm[Expand[Expand[x*Qin]]/.invariants];
     APsum1 = APsum1 //. Pm[a_]^n_.*Pm[b_]^m_.:>Pm[Expand[a^n*b^m]]
	             //. Pm[a_]^n_.:>Pm[Expand[a^n]];

     Spart = APsum1 /.Pm[x_]:>ToList[x] 
	             /.gmn[x_]:>g[x /.invariants]
	             /.g[x_]->x; 

  
     Label[NumkEND];
      
  Return[Spart]]  


(***********************************************)
(*           Qin, jpol, tt                     *)
(***********************************************)
ToPoly[integral_, momentum_] := 
  Module[{jpol, ttpol, Qin, denom, integrand},

    denom = List @@ integral;
    integrand =
      Sum[X[i]*denom[[i]] /.PR[k_,m_,_]->k^2-m^2,{i,Length[denom]}];
    ttpol = Coefficient[integrand, momentum*momentum];
    jpol = integrand - momentum*ttpol*momentum;
    Qin = Expand[-1/2*Coefficient[jpol, momentum]];
    jpol = Expand[jpol + 2*momentum*Qin];
   

 Return[{ttpol,jpol,Qin}]]


(***********************************************)
(*            Finds F function                 *)
(***********************************************)
FUpolyfun[integral_, momentum_, invariants_] := 
  Module[{jpol, ttpol, Qin, denom, integrand, fpoly},


    {ttpol,jpol,Qin}=ToPoly[integral,momentum];
    Upoly = ttpol;
    fpoly = Expand[Simplify[ttpol*(Expand[Qin*1/ttpol*Qin] - jpol)]];


	 If[momentum == Last[FullInt[[3]]],
	 fpoly = Expand[fpoly /. invariants],
	 {},
       fpoly = Expand[PRfind[fpoly,invariants] /.invariants]
    ];

    masslist = Union[List @@ integral /. PR[_, m_, _] -> m];
    masslist = DeleteCases[masslist, 0];
    Do[
  	caslist = Cases[fpoly, masslist[[i]]^2  X[n_]^2];
  	If[Length[caslist] >= 2,
	  {
      	    casexpr = caslist /. a_^2 -> a;
      	    expr = Simplify[Plus @@ casexpr];
      	    FXexpr = expr /. x_.*(y__) -> x^2*FX[y]^2;
      	    FXrule = {First[caslist]->FXexpr-(expr)^2+First[caslist]};
      	    fpoly = Expand[fpoly /. FXrule]; 
      	  },{}],
    {i, 1, Length[masslist]}];

    Print["U polynomial..."];
    Print[Upoly];
    Print["F polynomial..."];
 

 Return[fpoly]]


(*************************************************)
(*          Creates PRs in F polynomial          *)
(*************************************************)
PRfind[fpoly_,invariants_]:=
  Module[{intmom,masrul,nomass,preprop,p1,p2,ints,intrul,F,
	  PRnom,nom,PRmass,rul1,rul2,nm,PRadd,aux},
   
     
    intmom = Alternatives @@ (a_. Last[FullInt]);
    masrul = DeleteCases[Union[List @@ integral /.PR[_,a_,_]->a],0];
    nomass =
      Simplify[#]& /@ Collect[fpoly/.(#->0& /@ masrul),X[a_]X[b_]];


    preprop = Plus @@ Cases[nomass+aux,(intmom)^2*_|((intmom)+__)^2*_];
 
	   
    p1 = preprop /.(x_)^2->PR[x];
    p2 = Expand[fpoly-preprop];

    ints= Propagators;
    intrul = ints /.PR[a_,b_,_]->(PR[a]->PR[a,b]+b^2);
    intrul = DeleteCases[List @@ (intrul*(aux->0)),aux->0];

    F = Expand[p1+p2 /.intrul /.invariants];       
    
    PRnom = Cases[F,x_. PR[a_]];
    If[Length[PRnom] == 0,{},
    {
      nom = PRnom /. w_. PR[a_]X[n_]X[m_]->X[n]X[m];
      PRMass = Cases[F,Alternatives @@ (a_. m_^2*nom)];
      rul1 = PR[x_]-m_^2->PR[x,m];
      rul2 =-PR[x_]+m_^2->-PR[x,m];
      nm = Plus @@ Join[PRnom, PRMass];
      PRadd = Collect[nm,nom] /.{rul1,rul2};
      F=F-nm+PRadd
    }];
    
   
 Return[Expand[F /. PR[x_]->PR[x,0]]]]


(*************************************************)
(*          Small Barnes 1 Lemma                 *)
(*************************************************)
BarnesLemma[repres_,num_,option___Rule]:=
    Module[{options, 
	    ZETS, RSC, Zcheck, 
	    BL1, BL2, BL3, BL4, 
	    B2a,
	    Bmes1, Bmes2, Bmes3, BLru,
	    pow, Zpowers, ShiftRul, powlist, gamlist,
	    varlist, dimrepr,reprList,
	    before, after,
            counter=0, BList={},zel, denR, ifden, numR, ifnum, result},


       options=If[Length[{option}]==0,Option[BarnesLemma],{option}];

    ZETS[zlst_]:=
       Module[{FML,result},
	   FML[x_]:=Flatten[Map[(List @@ #)&, x]];
           result=DeleteCases[FML[FML[zlst]],_Integer]//Union;
       Return[result]];



       RSC[x_,y_]:=ToExprStr[z,#]& /@ 
	  Reverse[Sort[
            ToExpression[StringDrop[ToString[#], 1]]& /@ 
	    Complement[x,y,{}]]];

       Zcheck[x_,z_]:=Cases[x,Gamma[a_.*z+__]|Gamma[a_.*z]];

	  

    (* Using 1st or 2nd Barnes-Lemma *)

       Switch[num,
	      1,
	      
	 BL1[z_]:=
	   Gamma[z+(l1_.)]*Gamma[z+(l2_.)]*Gamma[-z+(l3_.)]*Gamma[-z+(l4_.)]
	   ->Gamma[l1+l3]*Gamma[l1+l4]*Gamma[l2+l3]*
	   (Gamma[l2+l4]/Gamma[l1+l2+l3+l4]);

	 BL2[z_]:=
	   Gamma[z+(l1_.)]^2*Gamma[-z+(l3_.)]*Gamma[-z+(l4_.)] 
	   ->Gamma[l1+l3]^2*Gamma[l1+l4]^2/Gamma[2*l1+l3+l4];

	 BL3[z_]:=
	   Gamma[z+(l1_.)]*Gamma[z+(l2_.)]*Gamma[-z+(l3_.)]^2 
	   ->Gamma[l1+l3]^2*Gamma[l2+l3]^2/Gamma[l1+l2+2*l3];

	 BL4[z_]:=
	   Gamma[z+(l1_.)]^2*Gamma[-z+(l3_.)]^2
	   ->Gamma[l1+l3]^4/Gamma[2*l1+2*l3];

	 Bmes1 = ">> Barnes 1st Lemma will be checked for: ";
	 Bmes2 = ">> Representation after 1st Barnes Lemma:  <<";
	 Bmes3 = "   1st Barnes Lemma was applied for: ";
	 BLrul = Dispatch[{BL1[zel],BL2[zel],BL3[zel],BL4[zel]}],
     (*----------------------------------------*)	      
	      2,

	 B2a[z_]:=
	   Gamma[(l1_.)+z]*Gamma[(l2_.)+z]*Gamma[(l3_.)+z]*
	   Gamma[(l4_.)-z]*Gamma[(l5_.)-z]/Gamma[(l6_.)+z]
	   :>Gamma[l1+l4]*Gamma[l2+l4]*Gamma[l3+l4]*Gamma[l1+l5]*
	   Gamma[l2+l5]*Gamma[l3+l5]/(Gamma[l1+l2+l4+l5]*
	   Gamma[l1+l3+l4+l5]*Gamma[l2+l3+l4+l5]) 
	   /;Expand[l1 + l2 + l3 + l4 + l5 - l6] === 0;
 

	 Bmes1 = ">> Barnes 2nd Lemma will be checked for: ";
	 Bmes2 = ">> Representation after 2nd Barnes Lemma:  <<";
	 Bmes3 = "   2nd Barnes Lemma was applied for: ";
	 BLrul = Dispatch[{B2a[zel]}]

       ];
      

    (* Preparing representation *)

       reprList=ToList[Expand[repres/.Gamma[x_]:>Gamma[Expand[x]]]];	
    

    (* Creating list of zs NOT to be used *)

       Zpowers=LSTofZ[reprList,1];
      

(*---------------------------*)
(*      Shifting part...     *)
       

       If[(Shifts /. options),
	  ShiftRul=SHIFT[Zpowers];
	  If[Length[ShiftRul]!=0,
	     Print[">>",SFred[" Shifting: "],ShiftRul];
	     reprList=(#/.ShiftRul/.Gamma[x_]:>
		       Gamma[Expand[x]])& /@ reprList; 
	     Zpowers=LSTofZ[reprList,1],{}]
	  ,{}];

(*---------------------------*)

       powlist=ZETS[Zpowers];


   (* Creating list of zs which appear in Gammas *)

       gamlist=ZETS[LSTofZ[reprList,2]];
	
	       
   (* Creating list of zs to check with B-L*)

       varlist=RSC[gamlist,powlist];   
       dimrepr=Length[gamlist];


   (* Beggining of checking *)

       Print[Bmes1,varlist," <<"];
       Print["   Starting with dim=",dimrepr," representation..."];
       
       result=reprList;
       Do[
	  zel= varlist[[i]];
	  WriteString[$Output,"\n"];
	  WriteString[$Output,i,". ","Checking ",zel];
	  
	  denR = Denominator[First[result]];
	  ifden = Length[Zcheck[denR,zel]];

	  numR = Numerator[First[result]];
	  ifnum = Length[Zcheck[numR,zel]];

	  If[num==1 && ifden!=0,Goto["break"]];
	  If[num==1 && ifnum>4,Goto["break"]];
	  If[num==2 && ifden==0,Goto["break"]];
	  If[num==2 && ifnum>5,Goto["break"]];

	  before=result;
	  result=(#/. BLrul)& /@ result;
	  after=result;

	  If[before=!=after,
	  {
	    counter=++counter;
	    BList=Prepend[BList,zel];
	    WriteString[$Output,"...Barnes Lemma was applied."]
	  },{}];
	 
  Label["break"],

       {i,1,Length[varlist]}];

       Print[Bmes2];

       If[Length[BList]==0,
	  Print[SFred["   Could not apply Barnes-Lemma"]],
       {
	  Print[SFred[Bmes3],BList];
	  Print["   Obtained representation has: dim=",dimrepr-counter]
       }];

 Return[Plus @@ result]]

LSTofZ[xli_,contr_] :=
  Module[{powers,pow,result},
    

     powers=First[FullInt[[2]]]/.PR[_,_,n_]:>Abs[n];
     powers=powers/. Abs[x_]->x;
     powers=#->0& /@ DeleteCases[(List @@ powers),_Integer];


   Switch[contr,
	  1,
     pow=First[xli]/.Gamma[_]->1 /.eps->0/.powers,
          2,
     rul={1/Gamma[x__]^a_.->x,Gamma[x_]^a_.->x};
     pow=Cases[First[xli],Gamma[__]^n_.]/.rul /.eps->0 /.powers
   ];


     result=DeleteCases[
		 If[Head[pow]===Times || Head[pow]===List,
		    List @@ pow,{pow}] /._^x_ :> 
			  Expand[x],_Integer|_Real|_Rational,2];

     If[result==1,result={},{},{}];


 Return[result]]

SHIFT[zlist_] := Module[{SD, prep, crosTab, comel, nocomel, ord, result},
    
    SD[x_] := ToExpression[StringDrop[ToString[x], 1]];
    INT[x_, a_] := 
      Table[If[i == j, {0}, 
          Intersection[Part[x, i], Expand[(a)*Part[x, j]]]], {i, 1, 
          Length[x]}, {j, 1, i}];
    
    	prep = DeleteCases[zlist, _Symbol];

    	crosTab = DeleteCases[INT[prep, -1] // Flatten, 0 | -1, 3];   
    	If[Length[Tab] == 1, GOTO["direct1"], {}];
    
    	comel = DeleteCases[INT[crosTab, 1] // Flatten, 0];
    	If[Length[comel] == 0, GOTO["direct1"], {}];
    
    	nocomel = 
      DeleteCases[
        DeleteCases[#, Alternatives @@ comel] & /@ crosTab, _Symbol];
    	ord = Sort[#, SD[#1] < SD[#2] &] & /@ nocomel;
  
    
    Label["direct1"];
    	result = (First[#] -> First[#] - (Take[#, 1 - Length[#]])) & /@ ord;
    
    Return[result]]

End[]


EndPackage[]
