creating tadpole vector more efficiently

by avoiding frequent calls to loop functions.
FlexibleSUSY · Mar 27, 2017 · 059719f · 059719f
1 parent 2f1d75c
commit 059719f
Show file tree

Hide file tree

Showing 3 changed files with 79 additions and 8 deletions.
diff --git a/meta/CConversion.m b/meta/CConversion.m
@@ -37,6 +37,8 @@
 AbsSqrt::usage="";
 FSKroneckerDelta::usage="";
 TensorProd::usage="";
+UV::usage="unit vector";
+Eval::usage="calls .eval()";
 
 HaveSameDimension::usage = "Checks if given types have same
 dimension";

diff --git a/meta/SelfEnergies.m b/meta/SelfEnergies.m
@@ -446,6 +446,56 @@ therefore not be accessed in the form Glu(gO2).
 ExpressionToStringSequentially[expr_, heads_, result_String] :=
     result <> " = " <> ExpressionToString[expr, heads] <> ";\n";
 
+PrepareExpr[expr_, vertexRules_] :=
+    expr /.
+    vertexRules /.
+    a_[List[i__]] :> a[i] /.
+    ReplaceGhosts[FlexibleSUSY`FSEigenstates] /.
+    C -> 1;
+
+RefactorSums[expr_] := SumOverToSum @ RecordSumCosts @ Expand @ SumToSumOver @
+		       expr;
+
+RecordSumCosts[expr_] := expr //.
+    SumOver[idx_, a_, b_] x_ :> SumOver[idx, a, b, IndexCost[idx, x]] x
+
+SumToSumOver[expr_] := expr //.
+    SARAH`sum[idx_, a_, b_, x_] :> SumOver[idx, a, b] x
+
+SumOverToSum[prod : SumOver[_,_,_,_] _] := Module[{
+	lst = List @@ prod,
+	sumOverToConvert,
+	idx, a, b, summand
+    },
+    sumOverToConvert = First @ SortBy[Cases[lst, SumOver[_,_,_,_]], Last];
+    {idx, a, b} = Take[List @@ sumOverToConvert, 3];
+    summand = Select[DeleteCases[lst, sumOverToConvert],
+		     !FreeQ[#, idx]&];
+    SumOverToSum[SARAH`sum[idx, a, b, Eval[SumOverToSum[Times @@ summand]]]
+		 Times @@ Complement[lst, {sumOverToConvert}, summand]]
+]
+
+SumOverToSum[x_Plus] := SumOverToSum /@ x;
+
+SumOverToSum[x_] := x;
+
+IndexCost[idx_, x:_[__]] /;
+    ValueQ[FunctionCost[x]] && !FreeQ[x, idx] := FunctionCost[x];
+
+IndexCost[idx_, x:_[__]] := Plus @@ (IndexCost[idx, #]& /@ List @@ x);
+
+IndexCost[idx_, _] := 0;
+
+(* cost function *)
+FunctionCost[SARAH`A0[_]]  := 1;
+FunctionCost[SARAH`B0[__]] := 2;
+FunctionCost[SARAH`B1[__]] := 2;
+FunctionCost[SARAH`B00[__]] := 2;
+FunctionCost[SARAH`B22[__]] := 2;
+FunctionCost[SARAH`F0[__]] := 2;
+FunctionCost[SARAH`G0[__]] := 2;
+FunctionCost[SARAH`H0[__]] := 2;
+
 CreateNPointFunction[nPointFunction_, vertexRules_List] :=
     Module[{decl, expr, prototype, body, functionName},
            expr = GetExpression[nPointFunction];
@@ -454,19 +504,36 @@ therefore not be accessed in the form Glu(gO2).
            prototype = type <> " " <> functionName <> ";\n";
            decl = "\n" <> type <> " CLASSNAME::" <> functionName <> "\n{\n";
            body = type <> " result;\n\n" <>
-                  ExpressionToStringSequentially[expr /.
-                                     vertexRules /.
-                                     a_[List[i__]] :> a[i] /.
-                                     ReplaceGhosts[FlexibleSUSY`FSEigenstates] /.
-                                     C -> 1,
+                  ExpressionToStringSequentially[
+                                     PrepareExpr[expr, vertexRules],
                                      TreeMasses`GetParticles[], "result"]  <>
                   "\nreturn result * oneOver16PiSqr;";
            body = IndentText[WrapLines[body]];
            decl = decl <> body <> "\n}\n";
            Return[{prototype, decl}];
           ];
 
-CreateNPointFunctionMatrix[_SelfEnergies`Tadpole] := { "", "" };
+CreateNPointFunctionMatrix[tadpole_SelfEnergies`Tadpole, vertexRules_List] :=
+    Module[{dim = GetDimension[GetField[tadpole]],
+            field = GetField[tadpole],
+            expr = GetExpression[tadpole],
+            type = CConversion`CreateCType[GetTadpoleType[GetField[tadpole]]],
+            functionName = CreateTadpoleFunctionName[GetField[tadpole], 1] <> "() const",
+            prototype, decl
+           },
+           If[dim == 1, Return[{ "", "" }]];
+           prototype = type <> " " <> functionName <> ";\n";
+           expr = RefactorSums[SARAH`sum[SARAH`gO1, 1, dim, UV[dim,SARAH`gO1] expr]];
+           decl = "\n" <> type <> " CLASSNAME::" <> functionName <> "\n{\n" <>
+              IndentText[type <> " result;\n\n" <>
+                         ExpressionToStringSequentially[
+                             PrepareExpr[expr, vertexRules],
+                             TreeMasses`GetParticles[], "result"]  <>
+                         "\nreturn result * oneOver16PiSqr;"
+                        ] <>
+           "\n}\n";
+           {prototype, decl}
+          ];
 
 FillHermitianSelfEnergyMatrix[nPointFunction_, sym_String] :=
     Module[{field = GetField[nPointFunction], dim, name},
@@ -501,7 +568,7 @@ therefore not be accessed in the form Glu(gO2).
                 ]
           ];
 
-CreateNPointFunctionMatrix[nPointFunction_] :=
+CreateNPointFunctionMatrix[nPointFunction_, vertexRules_List] :=
     Module[{dim, functionName, type, prototype, def},
            dim = GetDimension[GetField[nPointFunction]];
            If[dim == 1, Return[{ "", "" }]];
@@ -536,7 +603,7 @@ therefore not be accessed in the form Glu(gO2).
                {p,d} = CreateNPointFunction[nPointFunctions[[k]], vertexFunctionNames];
                prototypes = prototypes <> p;
                defs = defs <> d;
-               {p,d} = CreateNPointFunctionMatrix[nPointFunctions[[k]]];
+               {p,d} = CreateNPointFunctionMatrix[nPointFunctions[[k]], vertexFunctionNames];
                prototypes = prototypes <> p;
                defs = defs <> d;
               ];

diff --git a/templates/mass_eigenstates.cpp.in b/templates/mass_eigenstates.cpp.in
@@ -64,6 +64,8 @@ namespace flexiblesusy {
 #define LOCALINPUT(parameter) input.parameter
 #define MODELPARAMETER(parameter) model.get_##parameter()
 #define EXTRAPARAMETER(parameter) model.get_##parameter()
+#define Eval(expr) (expr).eval()
+#define UV(N,i) Eigen::Matrix<std::complex<double>,N,1>::Unit(i)
 
 #define HIGGS_2LOOP_CORRECTION_AT_AS     two_loop_corrections.higgs_at_as
 #define HIGGS_2LOOP_CORRECTION_AB_AS     two_loop_corrections.higgs_ab_as