Skip to content

Commit

Permalink
unitarity limit for s->infinity
Browse files Browse the repository at this point in the history
  • Loading branch information
@wkotlarski
wkotlarski committed Mar 5, 2023
1 parent 7152d15 commit 8d88abd
Show file tree
Hide file tree
Showing 5 changed files with 34 additions and 226 deletions.
13 changes: 5 additions & 8 deletions meta/FlexibleSUSY.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2502,15 +2502,12 @@ corresponding tadpole is real or imaginary (only in models with CP
];

WriteUnitarityClass[files_List] :=
Module[{},
matrix = Unitarity`GetScatteringMatrix[];
Module[{res},
res = Unitarity`GetScatteringMatrix[];
WriteOut`ReplaceInFiles[files,
{"@scatteringPairsLength@" -> ToString@matrix[[2]],
"@skipZeros@" -> "if (!(" <> StringRiffle["(i==" <> ToString[First@#-1] <> " && j==" <> ToString[Last@#-1] <> ")"& /@ First@matrix, "||"] <> ")) continue;",
"@scatteringElements@" -> TextFormatting`IndentText[Last@matrix],
"@generationSizes@" -> ToString[{{#}& /@ matrix[[3]]}],
"@infiniteSMatrix@" -> matrix[[4]],
Sequence @@ GeneralReplacementRules[]
{"@scatteringPairsLength@" -> ToString@res[[1]],
"@infiniteSMatrix@" -> res[[2]],
Sequence @@ GeneralReplacementRules[]
}];
];

Expand Down
123 changes: 24 additions & 99 deletions meta/Unitarity.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,23 +20,13 @@
*)

BeginPackage["Unitarity`", {"SARAH`", "Parameters`", "CConversion`", "TreeMasses`"}];
BeginPackage["Unitarity`", {"SARAH`", "Parameters`", "CConversion`", "TreeMasses`", "TextFormatting`", "Vertices`", "CXXDiagrams`"}];

GetScatteringMatrix::usage = "";

Begin["`Private`"];

replacePMASS[expr_] := Module[{temp},
temp = ToString@CForm[expr];
(* in expressions from SARAH masses are denoted as pmass(X)
this converts them to proper C++ form *)
temp = StringReplace[temp, "pmass(" ~~ f:Except[")"].. ~~ "(" ~~ i:DigitCharacter.. ~~ "))" :> "context.mass<" <> f <> ">({" <> ToString[ToExpression[i]-1] <> "})"];
temp = StringReplace[temp, "pmass(" ~~ f:Except[")"].. ~~ "(" ~~ i:Except[")"].. ~~ "))" :> "context.mass<" <> f <> ">({" <> ToString[i] <> "-1})"];
temp = StringReplace[temp, "pmass(" ~~ f:Except[")"].. ~~ ")" :> "context.mass<" <> f <> ">({})"];
temp
];

InfiniteS[a0Input_, generationSizes_] := Module[{params = Parameters`FindAllParametersClassified[a0Input], paramsCPP, mixingCPP, a0},
InfiniteS[a0Input_, generationSizes_, FSScatteringPairs_] := Module[{params = Parameters`FindAllParametersClassified[a0Input], paramsCPP, mixingCPP, a0},
(* CPP definitions of parameters present in the expression *)
paramsCPP =
StringJoin[
Expand All @@ -48,109 +38,44 @@

(* replace input parameters with their FS names *)
a0 = a0Input /. Thread[(Parameters`FSModelParameters /. params) -> CConversion`ToValidCSymbol /@ (Parameters`FSModelParameters /. params)];
a0 = a0 /. Susyno`LieGroups`conj -> Conj;
a0 = a0 //. SARAH`sum[idx_, start_, stop_, exp_] :> FlexibleSUSY`SUM[idx, start, stop, exp];
(* only for test
a0 = a0 /. Delta[c1_?SarahColorIndexQ, c2_?SarahColorIndexQ] :> 1;*)
resultInfinite = "";
For[i=1, i<=Length[a0], i++,
For[j=i, j<=Length[a0[[i]]], j++,
If[a0[[i,j]] =!= 0,
resultInfinite = resultInfinite <> "// " <> ToString[scatteringPairs[[i]]] <> "->" <> ToString[scatteringPairs[[j]]] <> "\n" <>
If[generationSizes[[i,j,1]] > 1, "for (int in1 = 1; in1 <= " <> ToString[generationSizes[[i,j,1]]] <> "; ++in1) {\n", ""] <>
If[generationSizes[[i,j,2]] > 1, "for (int in2 = 1; in2 <= " <> ToString[generationSizes[[i,j,2]]] <> "; ++in2) {\n", ""] <>
If[generationSizes[[i,j,3]] > 1, "for (int out1 = 1; out1 <= " <> ToString[generationSizes[[i,j,3]]] <> "; ++out1) {\n", ""] <>
If[generationSizes[[i,j,4]] > 1, "for (int out2 = 1; out2 <= " <> ToString[generationSizes[[i,j,4]]] <> "; ++out2) {\n", ""] <>
"matrix.coeffRef(" <> ToString[i-1] <> ", " <> ToString[j-1] <> ") = std::max(std::abs(" <> ToString@CForm[a0[[i,j]]] <> "), matrix.coeff(" <> ToString[i-1] <> ", " <> ToString[j-1] <> "));\n" <>
resultInfinite = resultInfinite <> "// " <> ToString[FSScatteringPairs[[i]]] <> "->" <> ToString[FSScatteringPairs[[j]]] <> "\n" <>
If[generationSizes[[i,j,1]] > 1, "for (int in1=1; in1<=" <> CXXNameOfField[First@FSScatteringPairs[[i]] /. Susyno`LieGroups`conj -> Identity] <> "::numberOfGenerations; ++in1) {\n", ""] <>
If[generationSizes[[i,j,2]] > 1, "for (int in2=1; in2<=" <> CXXNameOfField[Last@FSScatteringPairs[[i]] /. Susyno`LieGroups`conj -> Identity] <> "::numberOfGenerations; ++in2) {\n", ""] <>
If[generationSizes[[i,j,3]] > 1, "for (int out1=1; out1<=" <> CXXNameOfField[First@FSScatteringPairs[[j]] /. Susyno`LieGroups`conj -> Identity] <> "::numberOfGenerations; ++out1) {\n", ""] <>
If[generationSizes[[i,j,4]] > 1, "for (int out2=1; out2<=" <> CXXNameOfField[Last@FSScatteringPairs[[j]] /. Susyno`LieGroups`conj -> Identity] <> "::numberOfGenerations; ++out2) {\n", ""] <>
"matrix.coeffRef(" <> ToString[i-1] <> ", " <> ToString[j-1] <> ") = std::max(\n" <> TextFormatting`IndentText["std::abs(std::real(" <> ToString@CForm[FullSimplify[a0[[i,j]]]] <> ")),\nmatrix.coeff(" <> ToString[i-1] <> ", " <> ToString[j-1] <> ")\n"] <> ");\n" <>
If[generationSizes[[i,j,1]] > 1, "}\n", ""] <>
If[generationSizes[[i,j,2]] > 1, "}\n", ""] <>
If[generationSizes[[i,j,3]] > 1, "}\n", ""] <>
If[generationSizes[[i,j,4]] > 1, "}\n", ""]
If[generationSizes[[i,j,4]] > 1, "}\n", ""] <> "\n"
]
]
];
paramsCPP <> mixingCPP <> resultInfinite
TextFormatting`IndentText[paramsCPP <> "\n" <> mixingCPP <> "\n" <> resultInfinite]
];

(* return a C++ lambda computing a given scattering eigenvalue in function of sqrt(s) *)
ExpressionToCPPLambda[expr__, istates_List, fstates_List] := Module[{params = Parameters`FindAllParametersClassified[expr], paramsCPP, mixingCPP},
GetScatteringMatrix[] := Module[{generationSizes, a0, a0InfiniteS, FSScatteringPairs},
InitUnitarity[];

(* CPP definitions of parameters present in the expression *)
paramsCPP =
StringJoin[
("const auto " <> ToString@CConversion`ToValidCSymbol[#] <>
" = model_.get_" <> ToString@CConversion`ToValidCSymbol[#] <> "();\n")& /@ (Parameters`FSModelParameters /. params)
];
(* definition of mixing matrices *)
mixingCPP = ("auto " <> ToString[#] <> " = [&model_] (int i, int j) { return model_.get_" <> ToString@CConversion`ToValidCSymbol[#] <> "(i-1,j-1); };\n")& /@ (Parameters`FSOutputParameters /. params);
(* only color neutral final states *)
FSScatteringPairs = Select[scatteringPairs, (TreeMasses`GetColorRepresentation /@ #) == {S,S}&];

(* replace input parameters with their FS names *)
newExpr = expr /. Thread[(Parameters`FSModelParameters /. params) -> CConversion`ToValidCSymbol /@ (Parameters`FSModelParameters /. params)];
newExpr = newExpr /. Susyno`LieGroups`conj -> Conj;
newExpr = newExpr //. SARAH`sum[idx_, start_, stop_, exp_] :> FlexibleSUSY`SUM[idx, start, stop, exp];

"[&model" <> If[!FreeQ[expr, sChan], ",sChan", ""] <> If[!FreeQ[expr, tChan], ",tChan", ""] <> If[!FreeQ[expr, uChan], ",uChan", ""] <> If[!FreeQ[expr, qChan], ",qChan", ""] <> "](double sqrtS, int in1, int in2, int out1, int out2) {
auto model_ = model;
// couplings should be evaluated at the renormalization scale sqrt(s)
// see comment around eq. 20 of 1805.07306
model_.run_to(sqrtS);
model_.solve_ewsb();
const double s = Sqr(sqrtS);
const " <> FlexibleSUSY`FSModelName <> "_cxx_diagrams::context_base context {model_};\n" <>
"if (sqrtS < context.mass<" <> ToString[fstates[[1]] /. Susyno`LieGroups`conj->Identity] <> ">(" <> If[TreeMasses`GetDimension[fstates[[1]]]>1, "{out1-1}", "{}"] <> ") + context.mass<" <> ToString[fstates[[2]] /. Susyno`LieGroups`conj->Identity] <> ">(" <> If[TreeMasses`GetDimension[fstates[[2]]]>1, "{out2-1}", "{}"] <> ")) return 0.;\n" <>
paramsCPP <>
mixingCPP <>
(* TODO: can this really be complex? *)
"double res = 0.;\n" <>
With[{temp = Simplify@Coefficient[newExpr, qChan]},
If[temp =!= 0,
"if (qChan) {
res += " <> replacePMASS@temp <> ";
}\n",
""
]] <>
With[{temp = Coefficient[newExpr, sChan]},
If[temp =!= 0,
"if (sChan) {
res += " <> If[FreeQ[temp, Susyno`LieGroups`conj], "", "std::real("] <> replacePMASS@temp <> If[FreeQ[temp, Susyno`LieGroups`conj], "", ")"] <> ";
}\n",
""
]] <>
With[{temp = Coefficient[newExpr, tChan]},
If[temp =!= 0,
"if (tChan) {
res += " <> If[FreeQ[temp, Susyno`LieGroups`conj], "", "std::real("] <> replacePMASS@temp <> If[FreeQ[temp, Susyno`LieGroups`conj], "", ")"] <> ";
}\n",
""
]] <>
With[{temp = Coefficient[newExpr, uChan]},
If[temp =!= 0,
"if (uChan) {
res += " <> If[FreeQ[temp, Susyno`LieGroups`conj], "", "std::real("] <> replacePMASS@temp <> If[FreeQ[temp, Susyno`LieGroups`conj], "", ")"] <> ";
}\n",
""
]] <>
" return " <> If[FreeQ[expr, Susyno`LieGroups`conj], "res", "std::real(res)"] <> ";
};\n\n"
];
a0 = Outer[GetScatteringDiagrams[#1 -> #2]&, FSScatteringPairs, FSScatteringPairs, 1];
(* obviously 's' lives in Susyno`LieGroups` *)
a0InfiniteS = Map[Limit[# /. qChan -> 1 /. sChan|tChan|uChan -> 0, Susyno`LieGroups`s->Infinity]&, a0, {2}];

GetScatteringMatrix[] := Module[{result, generationSizes, a0, a0InfiniteS},
InitUnitarity[];
(*RemoveParticlesFromScattering={Se , Sv, Sd, Su};*)
generationSizes = Table[{i, j}, {i, Length[FSScatteringPairs]}, {j, Length[FSScatteringPairs]}];
generationSizes = Apply[Join[TreeMasses`GetDimension /@ FSScatteringPairs[[#1]], TreeMasses`GetDimension /@ FSScatteringPairs[[#2]]]&, generationSizes, {2}];

a0 = Outer[GetScatteringDiagrams[#1 -> #2]&, scatteringPairs, scatteringPairs, 1];
a0InfiniteS = Map[Limit[# /. qChan -> 1 /. sChan|tChan|uChan -> 0, Susyno`LieGroups`s->Infinity]&, a0, {2}];
generationSizes = Table[{i, j}, {i,1, Length[scatteringPairs]}, {j,1, Length[scatteringPairs]}];
generationSizes = Apply[Join[TreeMasses`GetDimension /@ scatteringPairs[[#1]], TreeMasses`GetDimension /@ scatteringPairs[[#2]]]&, generationSizes, {2}];
resultFinite = "";
resultInfinite = "";
For[i=1, i<=Length[a0], i++,
For[j=i, j<=Length[a0[[i]]], j++,
If[a0[[i,j]] =!= 0,
resultFinite = resultFinite <> "// " <> ToString[scatteringPairs[[i]]] <> "->" <> ToString[scatteringPairs[[j]]] <> "\n" <>
"matrix[" <> ToString[i-1] <> "][" <> ToString[j-1] <> "] = " <> ExpressionToCPPLambda[a0[[i,j]], scatteringPairs[[i]], scatteringPairs[[j]]];
resultInfinite = resultInfinite <> "// " <> ToString[scatteringPairs[[i]]] <> "->" <> ToString[scatteringPairs[[j]]] <> "\n" <>
"matrix[" <> ToString[i-1] <> "][" <> ToString[j-1] <> "] = " <> ToString@CForm[a0InfiniteS[[i,j]]] <> ";\n";
]
]
];
{SparseArray[a0]["NonzeroPositions"], Dimensions[a0][[1]], generationSizes, InfiniteS[a0InfiniteS, generationSizes], resultFinite}
{Length[FSScatteringPairs], InfiniteS[a0InfiniteS, generationSizes, FSScatteringPairs]}
];

End[];
Expand Down
3 changes: 1 addition & 2 deletions templates/run.cpp.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,7 @@ int run_solver(flexiblesusy::@ModelName@_slha_io& slha_io,

@calculateDecaysForModel@

@ModelName@_unitarity::max_scattering_eigenvalue(std::get<0>(models));
@ModelName@_unitarity::max_scattering_eigenvalue_infinite_s(std::get<0>(models));
std::cout << @ModelName@_unitarity::max_scattering_eigenvalue_infinite_s(std::get<0>(models)) << std::endl;
const bool show_result = !problems.have_problem() ||
spectrum_generator_settings.get(Spectrum_generator_settings::force_output);
// SLHA output
Expand Down
120 changes: 4 additions & 116 deletions templates/unitarity.cpp.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,12 @@

#include "@ModelName@_unitarity.hpp"
#include "@ModelName@_mass_eigenstates.hpp"
#include "@ModelName@_context_base.hpp"
#include "cxx_qft/@ModelName@_fields.hpp"

#include "minimizer.hpp"
#include "wrappers.hpp"
#include "sum.hpp"
#include "wrappers.hpp"

#include <gsl/gsl_min.h>
#include <algorithm>
#include <array>
#include <Eigen/SparseCore>
#include <functional>

namespace flexiblesusy {
namespace @ModelName@_unitarity {
Expand All @@ -51,119 +46,12 @@ inline double Sqrt2(int i, int j) {
double max_scattering_eigenvalue_infinite_s(const @ModelName@_mass_eigenstates& model) {
Eigen::SparseMatrix<double> matrix(size, size);

@infiniteSMatrix@
using namespace @ModelName@_cxx_diagrams::fields;

@infiniteSMatrix@
matrix.makeCompressed();
return matrix.coeffs().maxCoeff();
}

double max_scattering_eigenvalue(const @ModelName@_mass_eigenstates& model) {
static constexpr int size = @scatteringPairsLength@;
std::array<std::array<std::function<double(double, int, int, int, int)>, size>, size> matrix = {};
const int sChan = 1;
const int tChan = 1;
const int uChan = 1;
const int qChan = 1;

using namespace @ModelName@_cxx_diagrams::fields;

@scatteringElements@
std::array<std::array<double, size>, size> max_eigenvalues = {};

constexpr std::array<std::array<std::array<int, 4>, size>, size> generationSizes =
@generationSizes@;

for (int i=0; i<size; i++) {
for (int j=i; j<size; j++) {
@skipZeros@
std::array<int, 4> externalGenerationSizes = generationSizes[i][j];
// suply indices in mathematica counting (starting from 1)
for (int iIn1 = 1; iIn1 <= externalGenerationSizes[0]; ++iIn1) {
for (int iIn2 = 1; iIn2 <= externalGenerationSizes[1]; ++iIn2) {
for (int iOut1 = 1; iOut1 <= externalGenerationSizes[2]; ++iOut1) {
for (int iOut2 = 1; iOut2 <= externalGenerationSizes[3]; ++iOut2) {
std::function<double(double)> f = std::bind(matrix[i][j], std::placeholders::_1, iIn1, iIn2, iOut1, iOut2);

int status;
int iter = 0, max_iter = 100;
const gsl_min_fminimizer_type *T;
gsl_min_fminimizer *s;
double a = 100.0, b = 5000.0;

static constexpr int grid_size = 10;
std::array<double, grid_size> grid;
for (int idx = 0; idx < grid_size; ++idx) { grid[idx] = a + (idx+1)*(b-a)/(grid.size()+1); }
std::array<double, grid_size> f_grid;
for (int idx = 0; idx < grid_size; ++idx) { f_grid[idx] = -std::abs(f(grid[idx])); }
double m = grid[std::distance(f_grid.begin(), std::min_element(f_grid.begin(), f_grid.end()))];
// the minimum is between [a,b]
if (-std::abs(f(m)) < std::min(-std::abs(f(a)), -std::abs(f(b)))) {
gsl_function F = {
[](double d, void* vf) -> double {
auto& f = *static_cast<std::function<double(double)>*>(vf);
return -std::abs(f(d));
},
&f
};
T = gsl_min_fminimizer_brent;
s = gsl_min_fminimizer_alloc (T);
gsl_min_fminimizer_set (s, &F, m, a, b);

printf ("using %s method\n",
gsl_min_fminimizer_name (s));

printf ("%5s [%9s, %9s] %9s %10s %9s\n",
"iter", "lower", "upper", "min",
"err", "err(est)");

printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
iter, a, b,
m, b - a);

do
{
iter++;
status = gsl_min_fminimizer_iterate (s);

m = gsl_min_fminimizer_x_minimum (s);
a = gsl_min_fminimizer_x_lower (s);
b = gsl_min_fminimizer_x_upper (s);

// determine sqrt(s) with a relative error < 0.1%
status
= gsl_min_test_interval (a, b, 0.0, 1e-3);

if (status == GSL_SUCCESS)
printf ("Converged:\n");

printf ("%5d [%.7f, %.7f] "
"%.7f %.7f\n",
iter, a, b,
m, b - a);
}
while (status == GSL_CONTINUE && iter < max_iter);
max_eigenvalues[i][j] = std::abs(f(m));

gsl_min_fminimizer_free (s);
}
// the minimum is at a or b
else {
max_eigenvalues[i][j] = std::max(std::abs(f(m)), max_eigenvalues[i][j]);
}
}
}
}}}}
std::cout << "=======================================\n";
double max_elem = 0.;
for (int i=0; i<size; i++) {
for (int j=i; j<size; j++) {
std::cout << max_eigenvalues[i][j] << std::endl;
max_elem = std::max(max_elem, max_eigenvalues[i][j]);
}
}
std::cout << max_elem << std::endl;
return max_elem;
}

} // namespace @ModelName@_unitarity
} // namespace flexiblesusy
1 change: 0 additions & 1 deletion templates/unitarity.hpp.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@ namespace flexiblesusy {
class @ModelName@_mass_eigenstates;

namespace @ModelName@_unitarity {
double max_scattering_eigenvalue(const @ModelName@_mass_eigenstates&);
double max_scattering_eigenvalue_infinite_s(const @ModelName@_mass_eigenstates&);

} // namespace @ModelName@_unitarity
Expand Down

0 comments on commit 8d88abd

Please sign in to comment.