(* ::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.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
BeginPackage["RGErun`"];


(* ::Input::Initialization:: *)
\[Beta]::usage="\[Beta][x,model,thresholds] is used to enter the beta function for the variable x in the model; the beta function may depend the thresholds \!\(\*
StyleBox[\"below\",\nFontSlant->\"Italic\"]\) the current energy scale if different effective field theories are used; if not, you can set it to _.";


(* ::Input::Initialization:: *)
run::usage="run[thresholds,variables,direction,model] runs the RGEs between thresholds in the direction (up or down), gives variables as functions of the energy scale and optionally matches variables at thresholds or prints out the values of the variables and their predefined functions at the thresholds together with a plot of running couplings.";


(* ::Input::Initialization:: *)
printOut::usage="printOut is an option for run which specifies whether the results and matched results at each scale should be printed out together with a plot of the running variables; the default value is printOut \[Rule] True";


(* ::Input::Initialization:: *)
matches::usage="matches is an option for run which specifies the matching conditions for the running variables at thresholds.";


(* ::Input::Initialization:: *)
functions::usage="functions is an option for run which specifies functions of the running variables (e.g. eigenvalues of a mass matrix) that should be printed out at thresholds.";


(* ::Input::Initialization:: *)
\[DoubleStruckCapitalI]::usage="\[DoubleStruckCapitalI] must be used in lieu of IdentityMatrix[n] in the beta functions.";


(* ::Input::Initialization:: *)
CirclePlus::usage="Use \[CirclePlus] instead of + to add a constant vector or matrix in a beta function.";


(* ::Input::Initialization:: *)
CircleTimes::usage="Use \[CircleTimes] instead of * to multiply element-wise by constant vector or matrix in a beta function.";


(* ::Input::Initialization:: *)
getRGEs::usage="getRGEs[name] imports the \[Beta]-functions defined in the output of SARAH (name is the directory of SARAH output) or PyR@TE (name is the file with PyR@TE Mathematica output).";


(* ::Input::Initialization:: *)
returnVariables::usage="returnVariables is an option for getRGEs that specifies whether to return a list of variables defined in a PyR@TE Mathematica output.";
replaceNames::usage="replaceNames is an option for getRGEs to replace variable names defined in the PyR@TE output.";


(* ::Input::Initialization:: *)
updateInstalled::usage="updateInstalled[directory] replaces the RGErun.m package in the user's Mathematica Applications folder with the one in the directory (the default is NotebookDirectory[]).";


(* ::Input::Initialization:: *)
Options[getRGEs]={returnVariables->False,replaceNames->{}};
getRGEs[name_:"",model_:"",OptionsPattern[]]:=Block[{variableSubs,AllRGEs,DimensionParameter,conj,Tp,Adj,trace,MatMul,pi},
Which[name==""&&And@@Map[Head[#]===List&,{BetaGauge,BetaLijkl,BetaYijk,BetaTijk,BetaMuij,BetaBij,BetaVEV}],(* SARAH definitions in the Global context *)
AllRGEs=
Join[BetaGauge,BetaLijkl,BetaYijk,BetaTijk,BetaMuij,BetaBij,BetaVEV],
DirectoryQ[name],(* Directory with SARAH output *)
AllRGEs=Flatten[Get[FileNameJoin[{name,#}]]&/@{"BetaGauge.m","BetaLijkl.m","BetaYijk.m","BetaTijk.m","BetaMuij.m","BetaBij.m","BetaVEV.m"},1],
True,(* PyR@TE output file *)
AllRGEs=Get[name]
];
If[Last[Dimensions[AllRGEs]]==3,(* Convert from SARAH format *){First[#],{1/(16 \[Pi]^2),(1/(16 \[Pi]^2))^2}.Rest[#]}&/@AllRGEs];
variableSubs=Map[#->Head[#]&,Select[AllRGEs[[All,1]],!AtomQ[#]&]];Scan[(\[Beta][First[#],model,_]:=Last[#])&,AllRGEs/.pi->\[Pi]/.variableSubs/.OptionValue[replaceNames]//.{conj->Conjugate,Tp->Transpose,Adj->ConjugateTranspose,trace[ms__]:>Tr[Dot[ms]],MatMul[ms__][__]:>Dot[ms]}];
If[OptionValue[returnVariables]==True,AllRGEs[[All,1]]/.variableSubs/.OptionValue[replaceNames]]
];


(* ::Input::Initialization:: *)
updateInstalled[directory_:NotebookDirectory[]]:=Block[{newVersion=FileNameJoin[{directory,"RGErun.m"}],
installedVersion=FileNameJoin[{$UserBaseDirectory,"Applications","RGErun.m"}]},
If[FileExistsQ[newVersion],
If[FileExistsQ[installedVersion],DeleteFile[installedVersion]];
CopyFile[newVersion,installedVersion];
Print[StringForm["Copied `1` to `2`.",newVersion,installedVersion]],
Print["No RGErun.m in ",directory]];
];


(* ::Input::Initialization:: *)
Begin["`Private`"];


(* ::Input::Initialization:: *)
SetAttributes[{\[DoubleStruckCapitalI]},Protected];
SetAttributes[{CirclePlus,CircleTimes},{HoldAll,Protected}];
chop\[Delta]=0 10^-20;
run\[Delta]=10^-30;


(* ::Input::Initialization:: *)
numericQ::usage="numericQ[expr] returns True if expr is a numeric quantity or a full array of numeric quantities and False otherwise.";
numericQ[x_]:=NumericQ[x]||ArrayQ[x,_,NumericQ];


(* ::Input::Initialization:: *)
flatten::usage="flatten[x] flattens x if x is a list, else returns {x}.";
flatten[xs_]:=If[Dimensions[N[xs]]=={},{xs},Flatten[xs]];


(* ::Input::Initialization:: *)
threadEqual::usage="threadEqual[xs,ys] returns equations formed from flatten'ed xs and ys";
threadEqual[xs_,ys_]:=Thread[flatten[xs]==flatten[ys]];


(* ::Input::Initialization:: *)
numericsNamesQ::usage="numericsNamesQ[variables,scale] returns the list of variables that are numeric at the scale.";
numericAtScaleQ[variables_,scale_]:=Select[variables,numericQ[#[scale]]&];


(* ::Input::Initialization:: *)
complexInterpolationFunctionQ::usage="complexInterpolationFunctionQ[if] returns True if the InterpolationFunction is complex anywhere and False otherwise.";
complexInterpolationFunctionQ[if_]:=Cases[Chop[Im[if[[4,3]]]],Except[0|0.]]!={};


(* ::Input::Initialization:: *)
split::usage="Split a list into parts such that the length of each part is given by the second argument.";
split[list_,parts_]:=Map[Take[list,#1+{1,0}]&,Partition[FoldList[Plus,0,parts],2,1]];
(* Taken from Programming With Mathematica, An Introduction by Paul Wellin *)


(* ::Input::Initialization:: *)
printSeparator::usage="printSeparator[n] prints a row of n \[LongDash]'s";
printSeparator[n_:50,c_:"\[LongDash]"]:=Print[StringJoin[ConstantArray[c,n]]];


(* ::Input::Initialization:: *)
style::usage="style[text,name] sets the text to a predefined style.";
style[text_,"header"]:=Style[text,FontSize->16];


(* ::Input::Initialization:: *)
display::usage="display[variable] returns the variable in MatrixForm if the variable is a vector or matrix, otherwise it returns the variable";
display[variable_]:=If[ArrayQ[variable],MatrixForm[Chop[N[variable],chop\[Delta]]],Chop[N[variable],chop\[Delta]]];


(* ::Input::Initialization:: *)
printIntermediateResults::usage="printIntermediateResults[thresholds, scale, variables, functions,done] prints out intermediate results in runup if printOut\[Rule]True; done specifies whether the variables are \"Run\" etc.";
printIntermediateResults[thresholds_,scale_,variables_,functions_,done_:"Run",unit_:"GeV"]:=Block[{},
Print["Results @ "<>StringJoin[Riffle[thresholds," = "]]<>" = ",N[scale]," ",unit];
Print[done<>" variables:"];
Print[Sequence@@Flatten[Riffle[{ToString[#]," = ",display[#[scale]]}&/@variables,
", "]]];
If[functions!={},
Print[done<>" functions:"];
Print[Sequence@@Flatten[Riffle[{ToString[#]," = ",display[#[scale]]}&/@functions,
", "]]];
];
];


(* ::Input::Initialization:: *)
checkInitialConditions::usage="checkInitialConditions[variables] checks whether all the initial conditions to run are defined and numerical at the scale.";
checkInitialConditions[variables_,scale_]:=With[{numerics=numericAtScaleQ[variables,scale]},If[numerics!=variables,Message[run::nonnum,ToString/@Complement[variables,numerics]]];
Complement[variables,numerics]];


(* ::Input::Initialization:: *)
groupThresholds::usage="groupThresholds[thresholds] takes as arguments the list of \"threshold name\" \[Rule] scale rules and returns the {sorted names, sorted mass scales} list.";
groupThresholds[thresholds_]:=Block[{thrs=GatherBy[SortBy[thresholds,Last],Last]},
{Map[First,thrs,{2}],Map[Last@First[#]&,thrs]}
];


(* ::Input::Initialization:: *)
ts::usage="ts[massScales] returns logs of ratios of mass levels from a vector.";
ts[massScales_]:={Log[#1/#2],run\[Delta]}&@@@Partition[massScales,2,1];(* The value 10^-30 instead of 0 as an upper bound avoids an NDSolve bug in v. 10.1 *)


(* ::Input::Initialization:: *)
runBetweenThresholds::usage="runBetweenThresholds[scale,tMin,tMax,thresholds,variables], where the initial values for the variables must be given at scale, solves the equations in the interval from tMin to tMax.";
runBetweenThresholds[scale_,tMin_,tMax_,model_,thresholds_,variables_,opts:OptionsPattern[]]:=Block[{eqs,vars,sub,t,sol,event,end},
sub=
ReleaseHold[With[{name=Symbol[ToString[#]<>"$"]},Hold[Sequence[#'[t]->Array[name[##]'[t]&,Dimensions[N[#[scale]]]],#->
Array[name[##][t]&,Dimensions[N[#[scale]]]]]]]&/@variables];
eqs=Flatten[ReleaseHold[Hold@Sequence[
threadEqual[#'[t]/.sub/.\[DoubleStruckCapitalI]:>IdentityMatrix[First[Dimensions[N[#[scale]]]]],\[Beta][#,model,thresholds]/.sub/.\[DoubleStruckCapitalI]:>IdentityMatrix[First[Dimensions[N[#[scale]]]]]/.{CirclePlus->Plus,CircleTimes->Times}],threadEqual[#/.sub/.t->tMin,#[scale]]]&/@variables]];
vars=Flatten[With[{name=Symbol[ToString[#]<>"$"]},Array[name[##]&,Dimensions[N[#[scale]]]]]&/@variables];
event[x_]:=Apply[And,Abs[#[x]]<4 \[Pi]&/@vars];
sol=First[NDSolve[eqs,vars,{t,tMin,tMax},Method->{"EventLocator","Event"->event[t]},FilterRules[{opts},Options[NDSolve]]]];
end =sol[[1,2,1,1,-1]];(* End of the domain of the interpolation function; running up! *)
{sol,end}
];


(* ::Input::Initialization:: *)
pieceTogether::usage="pieceTogether[variables,solutions,scales] pieces together each variable and its Derivatives as a function of the renormalisation scale from the solutions to the RGEs.";
pieceTogether[variables_,solutions_,scales_,model_]:=(
Scan[
(#[\[Mu]_?NumericQ,model]:=Block[{intervals=Interval/@Partition[scales,2,1]},
Piecewise[
With[{var=#,varS=Symbol[ToString[#]<>"$"]},
MapThread[{
 Hold[(* To avoid evaluation of functions in wrong intervals *)
With[{scale=#2[[1,2]]},
Array[
varS[##][Log[\[Mu]]-Log[scale]]&,
Dimensions[N[var[scales[[1]]]]]
]
]/.#1
],
IntervalMemberQ[#2,\[Mu]]
}&,
{solutions,intervals}]],
Which[\[Mu]<scales[[1]],#[scales[[1]]],\[Mu]>scales[[-1]],#[scales[[-1]]]]]//ReleaseHold
];
#[\[Mu]_?NumericQ]:=Evaluate[#][\[Mu],model];
)&,
variables];
Scan[
(Derivative[n_][#][\[Mu]_?NumericQ,model]:=Block[{intervals=Interval/@Partition[scales,2,1]},
Piecewise[
With[{var=#,varS=Symbol[ToString[#]<>"$"]},
MapThread[{
 Hold[(* To avoid evaluation of functions in wrong intervals *)
With[{scale=#2[[1,2]]},
Array[
1/\[Mu]^n StirlingS1[n,Range[n]].With[{is=##},Derivative[#][varS[is]][Log[\[Mu]]-Log[scale]]&/@Range[n]]&,
Dimensions[N[var[scales[[1]]]]]
]
]/.#1
],
IntervalMemberQ[#2,\[Mu]]
}&,
{solutions,intervals}]],
Array[0&,Dimensions[N[#[scales[[1]]]]]]]//ReleaseHold
];
Derivative[n_][#][\[Mu]_?NumericQ]:=Derivative[n][Evaluate[#]][\[Mu],model];
)&,
variables];
);


(* ::Input::Initialization:: *)
plotStyles={Black,ColorData["Legacy"]["LightCadmiumRed"],ColorData["Legacy"]["PermanentGreen"],ColorData["Crayola"]["Blue"],ColorData["Legacy"]["Apricot"],ColorData["Legacy"]["BlueViolet"],ColorData["Legacy"]["Gold"],Black};


(* ::Input::Initialization:: *)
(*optimalPosition::usage="optimalPosition[variables,minX,maxX] gives the rough position between minX and maxX where the variables are most distant from each other.";optimalPosition[variables_,areComplexes_,minX_,maxX_]:=Block[{dist},
dist[list_]:=Total[Abs[Subtract@@#]&/@Subsets[list,{2}]];
SortBy[Block[{x=#},{x,dist[Flatten@MapThread[If[And@@#2\[Equal]False,#1[x],ReIm[#1[x]]]&,{variables,areComplexes}]]}]&/@Rest[Most[Exp[Subdivide[Log[minX],Log[maxX],10]]]],Last]\[LeftDoubleBracket]-1,1\[RightDoubleBracket]
]*)


(* ::Input::Initialization:: *)
SetAttributes[{printOut},Protected];
Options[run]={printOut->True,matches->{},functions->{}}~Join~Options[NDSolve];
run::nonnum="The initial condition(s) for `1` is/are not numerical and/or not full arrays.";
run[thresholds_,variables_,direction_:"up",model_:"",opts:OptionsPattern[]]:=Block[{dir=Which[direction=="up",1,direction=="down",-1],thrs,scales,intervals,t,tEnd,thr,runthr,matchedHere,runBetween,piecewiseVariables,areComplexes,normalThresholds=Normal[thresholds](*,position*)},
{thrs,scales}=groupThresholds[normalThresholds];
intervals=Partition[scales,2,1];
intervals=Which[direction=="up",intervals,direction=="down",Reverse[intervals]];
Evaluate[thr/@scales]=thrs;
Evaluate[runthr/@scales]=FoldList[Join,First[thrs],Rest[thrs]];
(* Since there is no running at the first scale, it is undefined *)
Evaluate[t/@Which[direction=="up",Rest[scales],direction=="down",Most[scales]]]=ts[scales];
(* Run in each interval and print out results *)
runBetween=Block[{tMin,tMax,int=#2,runVariables},
If[checkInitialConditions[variables,int[[dir]]]=={},
{tMin,tMax}={t[int[[-dir]]][[dir]],t[int[[-dir]]][[-dir]]} ;
{runVariables,tEnd}=runBetweenThresholds[int[[dir]],tMin,tMax,model,runthr[int[[dir]]],variables,FilterRules[{opts},Options[NDSolve]]];

(*******************)
Print["tEnd = ",tEnd];
If[tEnd!=run\[Delta],Throw[tEnd]];


Scan[(#[int[[-dir]]]=With[{name=Symbol[ToString[#]<>"$"]},Array[name[##][tMax]&,Dimensions[N[#[scales[[dir]]]]]]]/.runVariables)&,variables];
If[OptionValue[printOut]==True,printSeparator[];
printIntermediateResults[thr[int[[-dir]]],int[[-dir]],variables,OptionValue[functions]]];
matchedHere={};

Scan[Module[{mtch=#[[2]][model,int[[-dir]],thr[int[[-dir]]],direction]},
If[numericQ[mtch],#[[1]][int[[-dir]]]=mtch;AppendTo[matchedHere,#[[1]]]]]&,OptionValue[matches]];


If[OptionValue[printOut]==True&&matchedHere!={},printIntermediateResults[thr[int[[-dir]]],int[[-dir]],matchedHere,OptionValue[functions],"Matched"];
printSeparator[]];
Append[#1,runVariables]
]
]&;

(* Do not evaluate anything if some initial conditions are ill defined: *)
If[checkInitialConditions[variables,scales[[dir]]]=={},
If[OptionValue[printOut]==True,
printSeparator[];
Print["RUNNING "<>ToUpperCase[direction]<>" IN MODEL "<>model];
 ];
If[OptionValue[printOut]==True,printIntermediateResults[thr[scales[[dir]]],scales[[dir]],variables,OptionValue[functions],"Initial"]];

Catch[piecewiseVariables=Fold[runBetween,{},intervals]]//Print;

If[OptionValue[printOut]==True,
printSeparator[]
 ];
If[direction=="down",piecewiseVariables=Reverse[piecewiseVariables]];

pieceTogether[variables,piecewiseVariables,scales,model];

(* Which variables have complex values? *)
areComplexes=Map[Or@@Map[complexInterpolationFunctionQ[Last[#]]&,#]&,Transpose[piecewiseVariables]];
areComplexes=split[areComplexes,Times@@Dimensions[#[scales[[1]]]]&/@variables];


If[OptionValue[printOut]==True&&checkInitialConditions[variables,scales[[-dir]]]=={},
(* Show a plot with all the running variables *)

(* Calculate the optimal position(s) for the labels *)
(*positions=Log[optimalPosition[variables,areComplexes,scales\[LeftDoubleBracket]1\[RightDoubleBracket],scales\[LeftDoubleBracket]-1\[RightDoubleBracket]]];*)
(*positions=optimalPositions[fs_,{{xMin_,xMax_},{yMin_,yMax_}}];*)

Print[
Show[Reverse@Sort@Flatten@MapThread[
Block[{variable=#1,areComplex=#2,colour=#3,position=RandomReal[Log[scales[[{1,-1}]]]],dimensions=Dimensions[N[#1[scales[[1]]]]],indices},
indices=If[dimensions=={},{""},Flatten[Array[Style[StringJoin[Riffle[ToString/@{##},","]],Red]&,dimensions]]];
MapIndexed[
{
LogLinearPlot[
With[{vari=flatten[variable[\[Mu]]][[Sequence@@#2]]},If[#1,{Re[vari],Im[vari]},vari]],{\[Mu],scales[[1]],scales[[-1]]},AxesLabel->{"\[Mu]/GeV",""},AxesOrigin->{scales[[1]],0},PlotStyle->colour,PlotRange->All,Method->{"AxesInFront"->False}],
Block[{index=#2},With[{vari=flatten[variable[E^position]][[Sequence@@index]],ind=indices[[Sequence@@index]]},
Graphics[Text[If[dimensions=={},First[#],Subscript[First[#],ind]],{position,Last[#]},Background->White]]&/@If[#1,{{StringForm["Re\[ThinSpace]``",variable],Re[vari]},{StringForm["Im\[ThinSpace]``",variable],Im[vari]}},{{ToString[variable],vari}}]]
]
}&,
areComplex]
]&,
{variables,areComplexes,Flatten[ConstantArray[plotStyles,Quotient[Length[variables],Length[plotStyles]]]]~Join~Take[plotStyles,Length[variables]~Mod~Length[plotStyles]]}
],
PlotRange->All]
];
printSeparator[]
 ];
];
];


(* ::Input::Initialization:: *)
End[];


(* ::Input::Initialization:: *)
EndPackage[];



