SimplifyEinsteinNotation.m

(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



BeginPackage["Susyno`SimplifyEinsteinNotation`",{"Susyno`ModelBuilding`"}]

(* Export the main functions: SimplifyEinsteinNotation and SimplifyModelExpression *)
{SimplifyEinsteinNotation,SimplifyModelExpression};

Begin["Private`"];


(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

FindMinHeadsPosition[list_]:=Module[{aux,cont},
cont=0;aux={};
While[aux==={},
cont++;aux=Position[list[[All,1]],cont,1];
];
Return[Flatten[aux]];
]

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

(* Result = {resIndex1,resIndex2, ...} with resIndexI = positions of indexI in the term: first coordinate is the head position (in the term) and the second is the position amongst the indices of that head *)
PositionOfIndices[term_,nIndices_]:=Module[{indices,result},
indices=Range[nIndices];
result=If[#==={},{},#[[All,{1,-1}]]]&/@Table[Position[term[[All,2]],index],{index,indices}];
Return[result];
]

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

(* Result = Index of the indices in the heads of a given term *)
IndexOfIndicesInTerm[term_]:=Module[{aux,sub,result},
Return[term[[All,2]]];
]

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

(* Chain loop - classify all heads connected to a leading head [=head which  is lexically ordered first] *)
BuildChain[init_,posI_,idxs_,scannedHeadsPositionsIn_]:=Module[{newHeadsPositions,scannedHeadsPositions,headP,headLabel,aux2,chain},
newHeadsPositions={init};
chain={};
scannedHeadsPositions=scannedHeadsPositionsIn;
While[Length[newHeadsPositions]>0,
headP=newHeadsPositions[[1]];
headLabel={headP};

(* Head indices/connections loop - get all the other heads that are connected with a given one *)
Do[
aux2=Flatten[DeleteCases[posI[[idxs[[headP,k]]]],{headP,k}]];
If[Length[aux2]==0,aux2={-1,idxs[[headP,k]]}]; (* Do this tricky if the index is free *)

(* IMPORTANT: this is not the correct procedure. We must not care about the position in the term of a given head (headP and posI) but instead we must convert this into the position of the heads along the chain which is being constructed. This conversion should be easy to do once the chain is completed but ideally (for speed reasons) we should not wait until the end of the chain's contruction to do this. *)
AppendTo[headLabel,aux2];
,{k,Length[idxs[[headP]]]}];

AppendTo[chain,headLabel];
AppendTo[scannedHeadsPositions,headP];

aux2=DeleteCases[headLabel[[2;;-1,1]],-1,{1}];
(* The Complement command cannot be used here because it sorts the results. With this solution (using Select) the problem is avoided *)
If[Length[aux2]>0,
newHeadsPositions=Select[DeleteDuplicates[Join[newHeadsPositions,aux2]],!MemberQ[scannedHeadsPositions,#]&],
newHeadsPositions=Select[newHeadsPositions,!MemberQ[scannedHeadsPositions,#]&]];
];
Return[chain];
]

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

FindChains[term_,nIndices_]:=Module[{scannedHeadsPositions,unscannedHeadsPositions,listOfheadLabels,init,posI,idxs,newHeadsPositions,headLabel,headP,firstHeadList,newChain,newListOfChains,aux,aux2, listOfTotalChains,newListOfTotalChains,headsInTerm},

(* Number of chains loop - a given term can have 1 or more disconnected set of heads *)
scannedHeadsPositions={};
listOfheadLabels={};

listOfTotalChains={{}};

posI=PositionOfIndices[term,nIndices];
idxs=IndexOfIndicesInTerm[term];

(* >>>>>>> Don't stop until all totalChain's in newListOfTotalChains are completed *)
While[Min[Length[Flatten[#,1]]&/@listOfTotalChains]<Length[term],
newListOfTotalChains={};

(* >>>>>>> Scan over the existing, incomplete totalChain's and try to complete them *)
Do[

If[Length[Flatten[totalChain,1]]==Length[term],
(* do noting to this totalChain because it is complete *)
AppendTo[newListOfTotalChains,totalChain];, 


scannedHeadsPositions=Flatten[totalChain[[All,All,1]]];
(* Use of complement is not a problem here *)
unscannedHeadsPositions=Complement[Range[Length[term]],scannedHeadsPositions];
firstHeadList=FindMinHeadsPosition[term[[unscannedHeadsPositions]]];

(* >>>>>>> Scan over possible first heads *)
Do[
init=unscannedHeadsPositions[[firstHead]]; 

(* Chain loop - classify all heads connected to a leading head [=head which  is lexically ordered first] *)
newChain=BuildChain[init,posI,idxs,scannedHeadsPositions];

AppendTo[newListOfTotalChains,Append[totalChain,newChain]];

,{firstHead,firstHeadList}];

];

,{totalChain,listOfTotalChains}];
listOfTotalChains=newListOfTotalChains;

];


(* XXXXXXXXXXXXXXXXXX  Adopt the final notation for the position of heads  XXXXXXXXXXXXXXXXXX *)

(* There are 3 different 'positions' for a head h: \!\(
\*SubsuperscriptBox[\(P\), \(H\), \(h\)] \[Equivalent] \(h\ is\ the\ position\ of\ the\ head\ in\ the\ heads' s\ list\ provided\ as\ input\)\);  \!\(
\*SubsuperscriptBox[\(P\), \(T\), \(h\)]\ is\ the\ position\ of\ the\ head\ in\ a\ given\ term\);  \!\(
\*SubsuperscriptBox[\(P\), \(C\), \(h\)]\ is\ the\ position\ of\ the\ head\ in\ a\ \(totalChain . \ In\)\ principle\), the chains start by having the notation {{\!\(
\*SubsuperscriptBox[\(P\), \(T\), \(h1\)], \({
\*SubsuperscriptBox[\(P\), \(T\), \(h1'\)], idx}\),  ... \)},{\!\(
\*SubsuperscriptBox[\(P\), \(T\), \(h2\)], \({
\*SubsuperscriptBox[\(P\), \(T\), \(h2'\)], idx}\),  ... \)}} but this must be converted to the form {{\!\(
\*SubsuperscriptBox[\(P\), \(H\), \(h1\)], \({
\*SubsuperscriptBox[\(P\), \(C\), \(h1'\)], idx}\),  ... \)},{\!\(
\*SubsuperscriptBox[\(P\), \(H\), \(h2\)], \({
\*SubsuperscriptBox[\(P\), \(C\), \(h2'\)], idx}\),  ... \)}} before any ordering or elimination of chains is done because only in this last notation are things independent of the dummy indices used to write the term originally *)

headsInTerm=term[[All,1]];
Do[
aux=Flatten[listOfTotalChains[[totalChainP]],1];

scannedHeadsPositions=Flatten[aux[[All,1]]];
aux2=Table[Position[scannedHeadsPositions,i][[1,1]],{i,Length[scannedHeadsPositions]}];
AppendTo[aux2,-1]; (* Trick to maintain the -1's *)

(* Do the first correction *)
aux[[All,1]]=headsInTerm[[#]]&/@aux[[All,1]];

(* Do the second correction *)
aux[[All,2;;-1,1]]=Map[aux2[[#]]&,aux[[All,2;;-1,1]],{2}];
(*
Print[totalChainP];
Print[aux[[All,2;;-1,1]]];
Print[];
*)
(* Replace the chain in the old format with the one in the new format *)
listOfTotalChains[[totalChainP]]=aux;

,{totalChainP,Length[listOfTotalChains]}];

(* Only the first chain matters (the minimal one) *)
listOfTotalChains=Sort[listOfTotalChains];

Return[Sort[listOfTotalChains][[1]]];
];

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

SimplifyEinsteinNotation[expr_,heads_,indices_]:=Module[{headsMod,saveThis,aux,numericalFactors,result,subIndices,removeBrakets},

(* Deal with lists of expressions *)
If[expr[[0]]===List,Return[SimplifyEinsteinNotation[#,heads,indices]&/@expr]];

aux=Expand[expr//.{Conjugate[a_+ b__]:>Conjugate[a]+Conjugate[Plus[b]],Conjugate[a_ b__]:>Conjugate[a]Conjugate[Times[b]]}];
If[aux[[0]]===Plus,aux[[0]]=List;removeBrakets=False;,aux={aux};removeBrakets=True;];
saveThis=aux;
headsMod=heads/.{BlankNullSequence:>x}/.x[]:>x___;

numericalFactors=Times@@Cases[#,Except[Alternatives@@heads],1]&/@aux;



aux=Flatten[Table[Cases[#,headsMod[[i]]:>{i,x},1],{i,Length[heads]}],1]&/@aux;

subIndices=MapThread[Rule,{indices,Range[Length[indices]]}];
aux=aux/.subIndices;

(* At this point, aux={term1,term2,...} where termI={head1,head2,...}. Each head is represented as {headPositionInListOfHeads,{<IndicesPositionInListOfIndices>}} so everything is converted to a position now *)
result=FindChains[#,Length[indices]]&/@aux;

(* TEMPORARY SOLUTION; IN FUTURE IT IS BEST TO BE ABLE TO GENERATE THE RESULTS FROM SCATCH AND IF NO 3-INDEX STRUCTURES ARE PRESENT THE RESULTS CAN BE SHOWN IN MATRIX FORM *)

aux=Table[{result[[i]],saveThis[[i]]/numericalFactors[[i]],i},{i,Length[result]}];
aux=Gather[aux,#1[[1]]==#2[[1]]&];
aux=Sum[Plus@@(numericalFactors[[aux[[i,All,3]]]])aux[[i,1,2]],{i,Length[aux]}];

(* If[removeBrakets,aux=aux[[1]]]; *)
Return[aux];
]

(* XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX *)

(* This is higher level function than SimplifyEinsteinNotation: it will automatically set the indices and heads parameters so that expressions involving parameters of a given model can be readily simplified *)
SimplifyModelExpression[expr_,model_]:=Module[{aux,heads,indices},

aux=Cases[Flatten[parameters[model],1],x_/;MemberQ[x,f[_],-1]];
aux[[All,-1]]=ConstantArray[BlankNullSequence[],Length[aux]];
heads=Append[Join[aux,Conjugate/@aux],KroneckerDelta[___]];
indices=Array[f,10];

Return[SimplifyEinsteinNotation[expr,heads,indices]];
]

End[];
EndPackage[];