Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

polyratfun and fill statements #454

Open
vsht opened this issue Jul 14, 2023 · 4 comments
Open

polyratfun and fill statements #454

vsht opened this issue Jul 14, 2023 · 4 comments

Comments

@vsht
Copy link

vsht commented Jul 14, 2023

Dear FORM experts,

I apologize if this has already been discussed/asked somewhere else, but I'm just wondering
whether there is a straightforward way to apply polyratfun to fill statements that are used to populate
a tablebase.

My use case is as follows: I do reduction with FIRE and employ a lightweight Mathematica script that
converts the content of the reduction tables into fill statements. Since the rational coefficients in front of the
masters can be quite complicated, I merely expand them (with Distribute) to make sure that each coefficient can be split into
a proper num and denom functions. So no factorization within Mathematica takes place and this task is
entirely offloaded to FORM. Here is an example of the output my script produces

*--#[ lsclFillStatements:
fill tabIBPtopology2(1, 1, 1, 1, 1, 1, -2) = lsclDen(gkin*(-2 + lsclD)*(-80 + 56*lsclD - 13*lsclD^2 + lsclD^3)*meta^3*u0b^3)*(lsclNum(-80) + lsclNum(102*lsclD) + lsclNum(-43*lsclD^2) + lsclNum(6*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0) + lsclDen((-2 + lsclD)*(20 - 9*lsclD + lsclD^2)*meta*u0b)*(lsclNum(10) + lsclNum(-9*lsclD) + lsclNum(2*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 1, 1, 1, 1, -1) = lsclDen(gkin^2*(-3 + lsclD)*(1280 - 1216*lsclD + 432*lsclD^2 - 68*lsclD^3 + 4*lsclD^4)*meta^4*u0b^4)*(lsclNum(-1680) + lsclNum(2174*lsclD) + lsclNum(-1051*lsclD^2) + lsclNum(225*lsclD^3) + lsclNum(-18*lsclD^4))*topology1(1, 0, 1, 0, 1, 0, 0) + lsclDen(gkin*(-3 + lsclD)*(-320 + 224*lsclD - 52*lsclD^2 + 4*lsclD^3)*meta^2*u0b^2)*(lsclNum(60) + lsclNum(-68*lsclD) + lsclNum(25*lsclD^2) + lsclNum(-3*lsclD^3))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 1, 1, 1, 1, 0) = lsclDen(gkin^3*(-4 + lsclD)*(160 - 72*lsclD + 8*lsclD^2)*meta^5*u0b^5)*(lsclNum(-720) + lsclNum(726*lsclD) + lsclNum(-243*lsclD^2) + lsclNum(27*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0) + lsclDen(gkin^2*(-4 + lsclD)*(-40 + 8*lsclD)*meta^3*u0b^3)*(lsclNum(90) + lsclNum(-57*lsclD) + lsclNum(9*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(-1, 1, 3, 1, 1, 1, 0) = lsclDen(2*gkin^3*(128 - 48*lsclD + 4*lsclD^2)*meta^3*u0b^3)*(lsclNum(240) + lsclNum(-162*lsclD) + lsclNum(27*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 2, 1, 1, 1, -1) = lsclDen(gkin^2*(-2 + lsclD)*(-384 + 256*lsclD - 56*lsclD^2 + 4*lsclD^3)*meta^3*u0b^3)*(lsclNum(-480) + lsclNum(564*lsclD) + lsclNum(-216*lsclD^2) + lsclNum(27*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(-1, 1, 2, 1, 1, 1, 0) = lsclDen(gkin^2*(-384 + 256*lsclD - 56*lsclD^2 + 4*lsclD^3)*meta^2*u0b^2)*(lsclNum(400) + lsclNum(-350*lsclD) + lsclNum(99*lsclD^2) + lsclNum(-9*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 1, 1, 1, 1, -1) = lsclDen(gkin*(-2 + lsclD)*(32 - 16*lsclD + 2*lsclD^2)*meta^2*u0b^2)*(lsclNum(-16) + lsclNum(14*lsclD) + lsclNum(-3*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 2, 1, 1, 1, 0) = lsclDen(gkin^3*(-3 + lsclD)*(192 - 80*lsclD + 8*lsclD^2)*meta^4*u0b^4)*(lsclNum(720) + lsclNum(-726*lsclD) + lsclNum(243*lsclD^2) + lsclNum(-27*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 1, 1, 1, 1, 0) = lsclDen(gkin^2*(-3 + lsclD)*(64 - 32*lsclD + 4*lsclD^2)*meta^3*u0b^3)*(lsclNum(-240) + lsclNum(242*lsclD) + lsclNum(-81*lsclD^2) + lsclNum(9*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, -1, 3, 1, 1, 1, -1) = 0;
fill tabIBPtopology2(1, 0, 2, 1, 1, 1, -2) = 0;
fill tabIBPtopology2(1, -1, 3, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 2, 1, 1, 1, -1) = 0;
fill tabIBPtopology2(1, -1, 2, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 1, 1, 1, 1, -1) = 0;
fill tabIBPtopology2(1, 0, 2, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 1, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 1, 0, 1, 1, 1, -1) = 0;
fill tabIBPtopology2(1, 1, -1, 1, 1, 1, 0) = lsclDen(gkin*(-4 + lsclD)*(-8 + 2*lsclD)*meta^3*u0b^3)*(lsclNum(16) + lsclNum(-14*lsclD) + lsclNum(3*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0) + lsclDen((-4 + lsclD)*meta*u0b)*(lsclNum(-6) + lsclNum(2*lsclD))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 0, 1, 1, 1, 0) = lsclDen(gkin^2*(-4 + lsclD)*(64 - 32*lsclD + 4*lsclD^2)*meta^4*u0b^4)*(lsclNum(240) + lsclNum(-242*lsclD) + lsclNum(81*lsclD^2) + lsclNum(-9*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0) + lsclDen(gkin*(-4 + lsclD)*(-8 + 2*lsclD)*meta^2*u0b^2)*(lsclNum(-30) + lsclNum(19*lsclD) + lsclNum(-3*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 2, 0, 1, 1, -1) = lsclDen(gkin*(-2 + lsclD)*(160 - 72*lsclD + 8*lsclD^2)*meta*u0b)*(lsclNum(-60) + lsclNum(48*lsclD) + lsclNum(-9*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 2, -1, 1, 1, 0) = lsclDen(gkin*(-2 + lsclD)*(-32 + 8*lsclD))*(lsclNum(-20) + lsclNum(16*lsclD) + lsclNum(-3*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 1, 0, 1, 1, -1) = lsclDen((-4 + lsclD)*(-2 + lsclD))*(lsclNum(-2) + lsclNum(lsclD))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 2, 0, 1, 1, 0) = lsclDen(gkin^2*(-3 + lsclD)*(-80 + 16*lsclD)*meta^2*u0b^2)*(lsclNum(90) + lsclNum(-57*lsclD) + lsclNum(9*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 1, 0, 1, 1, 0) = lsclDen(gkin*(-3 + lsclD)*(-16 + 4*lsclD)*meta*u0b)*(lsclNum(-30) + lsclNum(19*lsclD) + lsclNum(-3*lsclD^2))*topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 1, 2, 1, 0, 1, -1) = 0;
fill tabIBPtopology2(1, 1, 2, 1, 0, 1, 0) = 0;
fill tabIBPtopology2(1, 1, 1, 1, 0, 1, 0) = 0;
fill tabIBPtopology2(1, 1, 3, 1, 1, -1, -1) = lsclDen(2*gkin^2*(30720 - 38144*lsclD + 19520*lsclD^2 - 5264*lsclD^3 + 788*lsclD^4 - 62*lsclD^5 + 2*lsclD^6)*meta^2*u0b^2)*(lsclNum(16640) + lsclNum(-28992*lsclD) + lsclNum(19860*lsclD^2) + lsclNum(-6688*lsclD^3) + lsclNum(1107*lsclD^4) + lsclNum(-72*lsclD^5))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 2, 1, 1, 0, -2) = lsclDen(gkin*(-2 + lsclD)*(-1920 + 2144*lsclD - 952*lsclD^2 + 210*lsclD^3 - 23*lsclD^4 + lsclD^5)*meta^2*u0b^2)*(lsclNum(960) + lsclNum(-904*lsclD) + lsclNum(60*lsclD^2) + lsclNum(174*lsclD^3) + lsclNum(-61*lsclD^4) + lsclNum(6*lsclD^5))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 3, 1, 1, -1, 0) = lsclDen(gkin^3*(-4 + lsclD)*(-1920 + 944*lsclD - 152*lsclD^2 + 8*lsclD^3)*meta^3*u0b^3)*(lsclNum(1440) + lsclNum(-2172*lsclD) + lsclNum(1212*lsclD^2) + lsclNum(-297*lsclD^3) + lsclNum(27*lsclD^4))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 2, 1, 1, 0, -1) = lsclDen(gkin^2*(-3 + lsclD)*(-7680 + 8576*lsclD - 3808*lsclD^2 + 840*lsclD^3 - 92*lsclD^4 + 4*lsclD^5)*meta^3*u0b^3)*(lsclNum(15840) + lsclNum(-24852*lsclD) + lsclNum(15500*lsclD^2) + lsclNum(-4801*lsclD^3) + lsclNum(738*lsclD^4) + lsclNum(-45*lsclD^5))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 2, 1, 1, -1, 0) = lsclDen(gkin^2*(-4 + lsclD)*(-384 + 256*lsclD - 56*lsclD^2 + 4*lsclD^3)*meta^2*u0b^2)*(lsclNum(-480) + lsclNum(724*lsclD) + lsclNum(-404*lsclD^2) + lsclNum(99*lsclD^3) + lsclNum(-9*lsclD^4))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 1, 1, 1, 0, -1) = lsclDen(gkin*(-3 + lsclD)*(-128 + 96*lsclD - 24*lsclD^2 + 2*lsclD^3)*meta^2*u0b^2)*(lsclNum(-240) + lsclNum(242*lsclD) + lsclNum(-81*lsclD^2) + lsclNum(9*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 2, 1, 1, 0, 0) = lsclDen(gkin^3*(-4 + lsclD)*(240 - 88*lsclD + 8*lsclD^2)*meta^4*u0b^4)*(lsclNum(-720) + lsclNum(726*lsclD) + lsclNum(-243*lsclD^2) + lsclNum(27*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 1, 1, 1, 0, 0) = lsclDen(gkin^2*(-4 + lsclD)*(64 - 32*lsclD + 4*lsclD^2)*meta^3*u0b^3)*(lsclNum(240) + lsclNum(-242*lsclD) + lsclNum(81*lsclD^2) + lsclNum(-9*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 0, 3, 1, 1, 1, -1) = 0;
fill tabIBPtopology2(0, 0, 3, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(0, 0, 2, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(0, 1, 0, 1, 1, 1, 0) = lsclDen(gkin*(-3 + lsclD)*(-8 + 2*lsclD)*meta^2*u0b^2)*(lsclNum(-24) + lsclNum(17*lsclD) + lsclNum(-3*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 0, 0, 1, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 3, 0, 1, 1, -1) = 0;
fill tabIBPtopology2(1, 0, 3, 0, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 2, 0, 1, 1, 0) = 0;
fill tabIBPtopology2(1, 1, 0, 0, 1, 1, 0) = topology100(1, 1, 0, 0, 1, 1, 0);
fill tabIBPtopology2(1, 0, 3, 1, 0, 1, 0) = 0;
fill tabIBPtopology2(1, 0, 2, 1, 0, 1, 0) = 0;
fill tabIBPtopology2(0, 1, 3, 1, 1, 0, -1) = lsclDen(gkin^2*(16 - 2*lsclD)*(24 - 10*lsclD + lsclD^2)*meta^2*u0b^2)*(lsclNum(-240) + lsclNum(162*lsclD) + lsclNum(-27*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 3, 1, 1, 0, 0) = lsclDen(gkin^3*(32 - 8*lsclD)*(48 - 14*lsclD + lsclD^2)*meta^3*u0b^3)*(lsclNum(-960) + lsclNum(888*lsclD) + lsclNum(-270*lsclD^2) + lsclNum(27*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(0, 1, 2, 1, 1, 0, 0) = lsclDen(gkin^2*(28 - 8*lsclD)*(24 - 10*lsclD + lsclD^2)*meta^2*u0b^2)*(lsclNum(560) + lsclNum(-538*lsclD) + lsclNum(171*lsclD^2) + lsclNum(-18*lsclD^3))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 0, 3, 1, 1, 0, -1) = 0;
fill tabIBPtopology2(1, 0, 3, 1, 1, 0, 0) = 0;
fill tabIBPtopology2(1, 0, 2, 1, 1, 0, 0) = 0;
fill tabIBPtopology2(1, 1, 0, 1, 1, 0, 0) = lsclDen(gkin*(12 - 4*lsclD)*(-4 + lsclD)*meta^2*u0b^2)*(lsclNum(-48) + lsclNum(34*lsclD) + lsclNum(-6*lsclD^2))*topology1(1, 0, 1, 0, 1, 0, 0);
fill tabIBPtopology2(1, 1, 3, 0, 1, 0, -1) = 0;
fill tabIBPtopology2(1, 1, 3, 0, 1, 0, 0) = 0;
fill tabIBPtopology2(1, 1, 2, 0, 1, 0, 0) = 0;
fill tabIBPtopology2(1, 1, 3, 1, 0, 0, 0) = 0;
*--#] lsclFillStatements:

In this fashion I can apply polyratfun only after having substituted the reduction results into my amplitude. This is obviously not very efficient and I was thinking that it would have been much better to do this in advance.

On the other hand, converting each reduction rule into a local variable, applying polyratfun and then creating a fill statement out of it looks a bit too contrived to me. So perhaps there is a more natural way to integrate FIRE results into a FORM calculation that I'm not aware of.

Cheers,
Vlad

@jodavies
Copy link
Collaborator

jodavies commented Jul 14, 2023

Hi Vlad,

I usually do a little bit more in Mathematica before creating the Fill statements. I Collect against the integral function names on the RHS of the rules, and apply something like GCoeffSimplify2 to the coefficients.

ep does not end up in the prf, which helps stop the terms becoming too big for FORM (and I will later expand in ep in FORM with a routine which acts on the "denep" denominators).
If I am expanding in more symbols than just ep, I arrange that they end up in their own "den..." functions also. Then the prf functions contain only symbols which won't be expanded later.

Since performance on the Mathematica side can be an issue for large tables I split each reduction table into many parts and process each part in parallel.

(* Slightly modified, from Arnd Behring *)
crCollect[expr_,pat_,f_:Identity]:=
	Module[{patalt,vars,res},
		If[expr===0,
			Return[f[expr]]
		];
		patalt=If[Head[pat]===List,Alternatives@@pat,pat];
		vars=Union@Cases[expr,patalt,{0,Infinity}];
		If[Length[vars]===0,
			Return[f[expr]]
		];
		res=CoefficientRules[expr,vars];
		res=Map[Apply[Times,vars^#[[1]]]*f[#[[2]]]&,res];
		res=Total[res];
		Return[res]
	];


GCoeffSimplify2[x_] := Module[{numer,denom,coeff},
	coeff = Together[x];
	numer = Collect[Numerator[coeff], ep,
		num
	];
	denom = FactorList[Denominator[coeff]];
	denom = Map[
		(#/.{
			{a_/;ContainsAll[Variables[a],{ep}], p_}:>denep[a]^p,
                        {a_ /; ContainsNone[Variables[a], {ep}], p_} :> den[a]^p
		})&
		,
		denom
	];
	denom = denom/.List->Times;

	coeff = crCollect[numer*denom, {ep,_denep},
		(prf[Numerator[#],Denominator[#]]& @ Factor[Together[#/.{num[a_]->a,den[b_]->1/b}]])&
	]/.{denep[ep]->1/ep};
	Return[coeff];
];

As you suggest, you can alternatively create expressions instead of fill statements (wrapping the names in [] makes life much easier here...) and process them in FORM with InParallel; helping performance when processing many tiny expressions. Those results can be converted to fill statements with some simple regex.

@jodavies
Copy link
Collaborator

jodavies commented Jul 14, 2023

Assuming your lsclD = 4-2*ep, your second rule ends up something like:

((denep[1 + 2*ep]*prf[-45, 8*gkin^2*meta^4*u0b^4])/ep + 
   (denep[1 + 2*ep]*prf[-1, 4*gkin^2*meta^4*u0b^4])/ep^3 + 
   denep[1 + 2*ep]*prf[9, 2*gkin^2*meta^4*u0b^4] + 
   (denep[1 + 2*ep]*prf[17, 8*gkin^2*meta^4*u0b^4])/ep^2)*
  topology1[1, 0, 1, 0, 1, 0, 0] + 
 ((denep[1 + 2*ep]*prf[-1, gkin*meta^2*u0b^2])/ep + 
   (denep[1 + 2*ep]*prf[1, 4*gkin*meta^2*u0b^2])/ep^2 + 
   denep[1 + 2*ep]*prf[3, 4*gkin*meta^2*u0b^2])*
  topology100[1, 1, 0, 0, 1, 1, 0]

@vermaseren
Copy link
Owner

vermaseren commented Jul 15, 2023 via email

@vsht
Copy link
Author

vsht commented Jul 15, 2023

Hi Vlad,

I usually do a little bit more in Mathematica before creating the Fill statements. I Collect against the integral function names on the RHS of the rules, and apply something like GCoeffSimplify2 to the coefficients.

ep does not end up in the prf, which helps stop the terms becoming too big for FORM (and I will later expand in ep in FORM with a routine which acts on the "denep" denominators). If I am expanding in more symbols than just ep, I arrange that they end up in their own "den..." functions also. Then the prf functions contain only symbols which won't be expanded later.

Since performance on the Mathematica side can be an issue for large tables I split each reduction table into many parts and process each part in parallel.

Hi Josh,

many thanks for your input. At the moment I'm using FeynCalc`s Collect2 for organizing the prefactors and my own routine for loading the tables. The source codes for both can be found in the FeynCalc repo.

However, I noticed that the Mathematica parts of my setup are the one that require most time on the cluster, so I started to think of possibilities to minimize their amount in my code. At some point I would also like to make that code public and since not all universities actually have Mma site licenses, parallelizing Mathematica scripts might not be a viable option for everyone.

But InParallel indeed looks very interesting, I haven't really thought about that before.

Another way: : Do not construct the fill statements directly, but make an expression with the elements of the table mentioned as arguments of a function f; After that use the FillExpression statement to fill the table.

Dear Jos,

thanks a lot, this indeed sounds like an interesting solution, since somehow I overlooked fillexpression when browsing through the manual.

I think now I've been given enough food for thought to implement something that works reasonably well, thanks again to Josh and Jos :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants