2
$\begingroup$

I made a routine code CodeComp2 that performs algebraic evaluations using some data TablexgridWithEx, TablePhaseSpaceDecayProducts and a few analytic expressions: pproductLab1Spher, pproductLab3Spher, pproductLab2Spher, EproductLabSpher (the first group), and coordProj1, coordProj2 (the second group). Inside the code, some values that use the first group of expressions are evaluated first, and then they are inserted in the formulas from the second group. Because of some reason, the second step takes a lot of time, although the expressions from the first group itself look simpler than the first group.

What may be the reason for this slowdown, and how to fix the code?

These are the relevant analytic expressions:

phVal[px_, py_] = 
  If[py > 0, 
   ArcCos[px/(\[Sqrt](px^2 + py^2))], -ArcCos[
     px/(\[Sqrt](px^2 + py^2))]];
thVal[px_, py_, pz_] = ArcCos[pz/(\[Sqrt](px^2 + py^2 + pz^2))];
pMotherVecSpher[Ex_, mx_, thx_, phx_] = 
  Sqrt[Ex^2 - mx^2] {Sin[thx]*Cos[phx], Sin[thx]*Sin[phx], Cos[thx]};
vMotherVecSpher[Ex_, mx_, thx_, phx_] = +(
   pMotherVecSpher[Ex, mx, thx, phx]/Ex);
\[Gamma]Factor[Ex_, mx_] = Ex/mx;
\[CapitalGamma]factor[Ex_, mx_] = 
  Simplify[(\[Gamma]Factor[Ex, mx] - 
    1)/((v /. 
      Solve[\[Gamma]Factor[Ex, mx] == 1/Sqrt[1 - v^2], v])^2)[[1]]];
pproductRestVec[Eprodrest_, mprod_, thprodrest_, phprodrest_] = 
  Sqrt[Eprodrest^2 - mprod^2] {Sin[thprodrest]*Cos[phprodrest], 
    Sin[thprodrest]*Sin[phprodrest], Cos[thprodrest]};
pproductLabVecSpher[Ex_, mx_, thx_, phx_, Eprodrest_, mprod_, 
   thprodrest_, phprodrest_] = 
  Simplify[
   pproductRestVec[Eprodrest, mprod, thprodrest, 
     phprodrest] + \[Gamma]Factor[Ex, mx]*
     vMotherVecSpher[Ex, mx, thx, phx]*
     Eprodrest + \[CapitalGamma]factor[Ex, mx]*
     vMotherVecSpher[Ex, mx, thx, 
      phx] (vMotherVecSpher[Ex, mx, thx, phx] . 
       pproductRestVec[Eprodrest, mprod, thprodrest, phprodrest])];
EproductLabSpher[Ex_, mx_, thx_, phx_, Eprodrest_, mprod_, 
   thprodrest_, phprodrest_] = 
  Simplify[\[Gamma]Factor[Ex, 
     mx] (Eprodrest + 
      vMotherVecSpher[Ex, mx, thx, phx] . 
       pproductRestVec[Eprodrest, mprod, thprodrest, phprodrest])];
pproductLab1Spher[Ex_, mx_, thx_, phx_, Eprodrest_, mprod_, 
   thprodrest_, phprodrest_] = 
  pproductLabVecSpher[Ex, mx, thx, phx, Eprodrest, mprod, thprodrest, 
    phprodrest][[1]];
pproductLab2Spher[Ex_, mx_, thx_, phx_, Eprodrest_, mprod_, 
   thprodrest_, phprodrest_] = 
  pproductLabVecSpher[Ex, mx, thx, phx, Eprodrest, mprod, thprodrest, 
    phprodrest][[2]];
pproductLab3Spher[Ex_, mx_, thx_, phx_, Eprodrest_, mprod_, 
   thprodrest_, phprodrest_] = 
  pproductLabVecSpher[Ex, mx, thx, phx, Eprodrest, mprod, thprodrest, 
    phprodrest][[3]];
xyzPointxvertex[xLongx_, thx_, phx_] = {(xLongx/Sin[phx])*Cos[phx], 
   xLongx, (xLongx/(Sin[thx]*Sin[phx]))*Cos[thx]};
xyzPointxvertexX[xLongx_, thx_, phx_] = 
  xyzPointxvertex[xLongx, thx, phx][[1]];
xyzPointxvertexSecondCoordinate[xLongx_, thx_, phx_] = 
  xyzPointxvertex[xLongx, thx, phx][[3]];
xLongDecayProductFinalPlaneGivenExperiment = 85;
xyzPointDaughterProjection[xLongx_, thx_, phx_, px_, py_, 
   pz_] = {xyzPointxvertex[xLongx, thx, phx][[
     1]] + (xLongDecayProductFinalPlaneGivenExperiment - xLongx) (px/
       py), xLongDecayProductFinalPlaneGivenExperiment, 
   xyzPointxvertex[xLongx, thx, phx][[
     3]] + (xLongDecayProductFinalPlaneGivenExperiment - xLongx) (pz/
       py)};
coordProj1[xLongx_, thx_, phx_, px_, py_, pz_] = 
  If[py > 0, 
   Evaluate[
    xyzPointDaughterProjection[xLongx, thx, phx, px, py, pz][[1]]], 
   99999];
coordProj2[xLongx_, thx_, phx_, px_, py_, pz_] = 
  If[py > 0, 
   Evaluate[
    xyzPointDaughterProjection[xLongx, thx, phx, px, py, pz][[3]]], 
   99999];

Next, there is an example data:

mxval = 2;
TablexgridWithEx = 
  Tuples[{RandomReal[{mxval, 10^4}, 100] // Sort, 
    RandomReal[{0.33, 0.94}, 50] // Sort, 
    RandomReal[{60, 85}, 40] // Sort, 
    RandomReal[{-Pi, Pi}, 40] // Sort}];
indexEx = 1;
indexxLongx = 3;
indexthx = 3;
indexphx = 4;
indexEprod1rest = 1;
indexthprod1rest = 2;
indexphprod1rest = 3;
indexEprod2rest = 4;
indexthprod2rest = 5;
indexphprod2rest = 6;
mproduct1 = mproduct2 = 0.5;
nprod = 10;
RandomDir = RandomPoint[Sphere[], nprod];
TabProduct1Rest = 
  Join[Table[{(mxval^2 - mproduct2^2 + mproduct1^2)/(2 mxval)}, 
    nprod], {ArcCos[#[[3]]], phVal[#[[1]], #[[2]]]} & /@ RandomDir, 2];
TabProduct2Rest = 
  Join[Table[{(mxval^2 + mproduct2^2 - mproduct1^2)/(2 mxval)}, 
    nprod], {ArcCos[-#[[3]]], phVal[-#[[1]], -#[[2]]]} & /@ RandomDir,
    2];
TablePhaseSpaceDecayProducts = 
  Join[TabProduct1Rest, TabProduct2Rest, 2];

Finally, there is the block CodeComp2 compiling evaluations:

CodeComp2 = 
  Hold@Compile[{{TablexgridWithEx, _Real, 
                    1}, {TablePhaseSpaceDecayProducts, _Real, 
                    2}, {mx, _Real}, {mproduct1, _Real}, {mproduct2, \
_Real}},
                    
                    Module[{cond, count, FirstCoordinateProduct1, 
                    SecondCoordinateProduct1, FirstCoordinateProduct2,
                     SecondCoordinateProduct2, EnergyProduct1, 
                    EnergyProduct2, ex, xlongx, thx, phx, eprod1rest, 
                    thprod1rest, phprod1rest, eprod2rest, thprod2rest,
                     phprod2rest, pxprod1lab, pxprod2lab, pyprod1lab, 
                    pyprod2lab, pzprod1lab, pzprod2lab, \[Epsilon]Azx},
                    count = 0.;
                    ex = Compile`GetElement[TablexgridWithEx, indexEx];
                    
                    xlongx = 
                    Compile`GetElement[TablexgridWithEx, 
                    indexxLongx];
                    
                    thx = Compile`GetElement[TablexgridWithEx, 
                    indexthx];
                    
                    phx = Compile`GetElement[TablexgridWithEx, 
                    indexphx];
                    Do[
                    
                    eprod1rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexEprod1rest];
                    
                    thprod1rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexthprod1rest];
                    
                    phprod1rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexphprod1rest];
                    
                    eprod2rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexEprod2rest];
                    
                    thprod2rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexthprod2rest];
                    
                    phprod2rest = 
                    Compile`GetElement[TablePhaseSpaceDecayProducts, 
                    j, indexphprod2rest];
                    (*Where I use formulas that are evaluated fast*)
                    
                    pxprod1lab = 
                    pproductLab1Spher[ex, mx, thx, phx, eprod1rest, 
                    mproduct1, thprod1rest, phprod1rest];
                    
                    pyprod1lab = 
                    pproductLab2Spher[ex, mx, thx, phx, eprod1rest, 
                    mproduct1, thprod1rest, phprod1rest];
                    
                    pzprod1lab = 
                    pproductLab3Spher[ex, mx, thx, phx, eprod1rest, 
                    mproduct1, thprod1rest, phprod1rest];
                    
                    pxprod2lab = 
                    pproductLab1Spher[ex, mx, thx, phx, eprod2rest, 
                    mproduct2, thprod2rest, phprod2rest];
                    
                    pyprod2lab = 
                    pproductLab2Spher[ex, mx, thx, phx, eprod2rest, 
                    mproduct2, thprod2rest, phprod2rest];
                    
                    pzprod2lab = 
                    pproductLab3Spher[ex, mx, thx, phx, eprod2rest, 
                    mproduct2, thprod2rest, phprod2rest];
                    
                    EnergyProduct1 = 
                    EproductLabSpher[ex, mx, thx, phx, eprod1rest, 
                    mproduct1, thprod1rest, phprod1rest];
                    
                    EnergyProduct2 = 
                    EproductLabSpher[ex, mx, thx, phx, eprod2rest, 
                    mproduct2, thprod2rest, phprod2rest];
                    (*Where I use formulas that are evaluated slow*)
                    
                    FirstCoordinateProduct1 = 
                    coordProj1[xlongx, thx, phx, pxprod1lab, 
                    pyprod1lab, pzprod1lab];
                    
                    SecondCoordinateProduct1 = 
                    coordProj2[xlongx, thx, phx, pxprod1lab, 
                    pyprod1lab, pzprod1lab];
                    
                    FirstCoordinateProduct2 = 
                    coordProj1[xlongx, thx, phx, pxprod2lab, 
                    pyprod2lab, pzprod2lab];
                    
                    SecondCoordinateProduct2 = 
                    coordProj2[xlongx, thx, phx, pxprod2lab, 
                    pyprod2lab, pzprod2lab]
                    , {j, 1, Length[TablePhaseSpaceDecayProducts]}];
                    {count/Length[TablePhaseSpaceDecayProducts]}
                    ], CompilationTarget -> "C", 
                    RuntimeOptions -> "Speed", 
                    RuntimeAttributes -> {Listable}, 
                    Parallelization -> True] /. 
                   DownValues@coordProj1 /. DownValues@coordProj2 /. 
                 DownValues@pproductLab1Spher /. 
                DownValues@pproductLab2Spher /. 
               DownValues@pproductLab3Spher /. 
              DownValues@EproductLabSpher /. OwnValues@indexEx /. 
            OwnValues@indexxLongx /. OwnValues@indexthx /. 
          OwnValues@indexphx /. OwnValues@indexEprod1rest /. 
        OwnValues@indexthprod1rest /. OwnValues@indexphprod1rest /. 
      OwnValues@indexEprod2rest /. OwnValues@indexthprod2rest /. 
    OwnValues@indexphprod2rest // ReleaseHold;

The block has rows evaluating values using the functions pproductLab1Spher,pproductLab2Spher,pproductLab3Spher,EproductLabSpher, which are then inserted in coordProj1,coordProj2. The first functions look much more complicated than the second ones:

enter image description here

Let us first comment the following rows calling coordProj1, coordProj2:

FirstCoordinateProduct1 = 
 coordProj1[xlongx, thx, phx, pprod1lab, py1lab, pz1lab];
SecondCoordinateProduct1 = 
 coordProj2[xlongx, thx, phx, pprod1lab, py1lab, pz1lab];
FirstCoordinateProduct2 = 
 coordProj1[xlongx, thx, phx, pprod2lab, py2lab, pz2lab];
SecondCoordinateProduct2 = 
 coordProj2[xlongx, thx, phx, pprod2lab, py2lab, pz2lab]

Then

CodeComp2[TablexgridWithEx, TablePhaseSpaceDecayProducts, mxval, 
   mproduct1, mproduct2]; // AbsoluteTiming

{0.307338, Null}

Next, let us uncomment them:

{2.80388, Null}

The slowdown is ~10 times!

I have checked that the If condition in coordProj is not the main reason for the performance problem. Moreover, the rows using the first group of formulas (pproductLab1Spher, etc.; starting from pxprod1lab = ) take negligible time to evaluate. To see this, let us comment the rows

pxprod1lab = pproductLab1Spher[ex, mx, thx, phx, eprod1rest, mproduct1, thprod1rest, phprod1rest];
pyprod1lab = pproductLab2Spher[ex, mx, thx, phx, eprod1rest, mproduct1, thprod1rest, phprod1rest];
pzprod1lab = pproductLab3Spher[ex, mx, thx, phx, eprod1rest, mproduct1, thprod1rest, phprod1rest];
pxprod2lab = pproductLab1Spher[ex, mx, thx, phx, eprod2rest, mproduct2, thprod2rest, phprod2rest];
pyprod2lab = pproductLab2Spher[ex, mx, thx, phx, eprod2rest, mproduct2, thprod2rest, phprod2rest];
pzprod2lab = pproductLab3Spher[ex, mx, thx, phx, eprod2rest, mproduct2, thprod2rest, phprod2rest];
EnergyProduct1 = EproductLabSpher[ex, mx, thx, phx, eprod1rest, mproduct1, thprod1rest, phprod1rest];
EnergyProduct2 = EproductLabSpher[ex, mx, thx, phx, eprod2rest, mproduct2, thprod2rest, phprod2rest];
(*Where I use formulas that are evaluated slow*)
FirstCoordinateProduct1 = coordProj1[xlongx, thx, phx, pxprod1lab, pyprod1lab, pzprod1lab];
SecondCoordinateProduct1 = coordProj2[xlongx, thx, phx, pxprod1lab, 
pyprod1lab, pzprod1lab];
FirstCoordinateProduct2 = coordProj1[xlongx, thx, phx, pxprod2lab, 
pyprod2lab, pzprod2lab];
SecondCoordinateProduct2 = coordProj2[xlongx, thx, phx, pxprod2lab, pyprod2lab, pzprod2lab]

{0.284096, Null}

-- comparable with the evaluation time if commenting only FirstCoordinateProduct1 = ....

So the relative slowdown caused by FirstCoordinateProduct1 = ... is actually much larger.

Edit

Following the answer by Domen.

Even if replacing the trigonometric formulas for coordProj1, coordProj2 with simple

coordProj1[xLongx_, thx_, phx_, px_, py_, pz_] =   xLongx*thx*phx*px*py*pz; 
coordProj2[xLongx_, thx_, phx_, px_, py_, pz_] =   xLongx*thx*phx*px*py*pz;

the evaluation of CodeComp2 with uncommented rows still takes much longer time (0.7 seconds) than the evaluation with the commented rows (0.3 seconds), which again means a huge slowdown required to evaluate, e.g., FirstCoordinateProduct1 compared to the evaluation of, e.g., pxprod1lab (up to a 100 times longer).

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

Although it looks like your second formulae are "simpler" than the first ones, they contain two problematic functions: Csc[x] and Cot[x], which somehow happen to be significantly slower when compiled than all Sin, Cos, Sqrt that are present in the previous formulae.

Let me convince you:

opts = {CompilationTarget -> "C", RuntimeOptions -> "Speed"};

sin = Compile[{{xs, _Real, 1}}, Do[Sin[x], {x, xs}], Evaluate@{opts}];
sqrt = Compile[{{xs, _Real, 1}}, Do[Sqrt[x], {x, xs}], Evaluate@{opts}];
cos = Compile[{{xs, _Real, 1}}, Do[Cos[x], {x, xs}], Evaluate@{opts}];
cot = Compile[{{xs, _Real, 1}}, Do[Cot[x], {x, xs}], Evaluate@{opts}];
csc = Compile[{{xs, _Real, 1}}, Do[Csc[x], {x, xs}], Evaluate@{opts}];
recSin = Compile[{{xs, _Real, 1}}, Do[1/Sin[x], {x, xs}], Evaluate@{opts}];
recTan = Compile[{{xs, _Real, 1}}, Do[1/Tan[x], {x, xs}], Evaluate@{opts}];

xs = RandomReal[1, 10000000];
{#, Module[{f = ToExpression[#]}, 
     First@AbsoluteTiming[f[xs];]]} & /@ {"sin", "cos", "sqrt", "cot",
    "csc", "recSin", "recTan"} // Column

sin     0.0000221
cos     2.1*10^-6
sqrt    1.8*10^-6
cot     0.223804
csc     0.153003
recSin  4.*10^-6
recTan  2.1*10^-6

Therefore, use SetSystemOptions["SimplificationOptions" -> "AutosimplifyTrigs" -> False] and use replacements* /. Cot[x_] :> 1/Tan[x] and /. Csc[x_] :> 1/Sin[x] in your formulae. I tried it on your code, and it worked.

* I am not an expert in numerical calculations, so you should check yourself that using reciprocal values does not affect the stability and accuracy of your calculations.

Improve this answer
$\endgroup$
12
  • $\begingroup$ Thanks! I will check. However, what still disturbs me is that even if replacing these trigonometric formulas with simple coordProj1[xLongx_, thx_, phx_, px_, py_, pz_] = xLongx*thx*phx*px*py*pz; coordProj2[xLongx_, thx_, phx_, px_, py_, pz_] = xLongx*thx*phx*px*py*pz;, the evaluation of the corresponding rows in CodeComp2 still takes much-much longer than the evaluation of the rows pxprod1lab = ... (~0.01 second for the latter vs 0.4 seconds for the former+the latter on my machine). Do you have any ideas about the reason? $\endgroup$
    – John Taylor
    Commented Apr 24, 2023 at 14:48
  • $\begingroup$ By 0.01 second, I mean the relative timing between the computation of CodeComp2 with commented rows FirstCoordinateProduct1 = ... and similar, and the one with, in addition, commented rows pxprod1lab = ... and similar. $\endgroup$
    – John Taylor
    Commented Apr 24, 2023 at 14:50
  • $\begingroup$ Hm, are you sure you cleared the variables properly? If I use this simple definitions for coordProj1 and coordProj2, I get practically the same results for both bases (with or without "slow" formulas). (It's quite a long code, so it might also be me evaluating some wrongly ...) $\endgroup$
    – Domen
    Commented Apr 24, 2023 at 14:52
  • 1
    $\begingroup$ I have once again copied your whole code and used simplified coordProj1/2 without trigs. There is still only a slight increase in time when I compare funcs with or without the problematic four lines. Am I misunderstanding the problem? $\endgroup$
    – Domen
    Commented Apr 24, 2023 at 15:18
  • 1
    $\begingroup$ No, it seems that I had a problem on my side related to cache maybe. $\endgroup$
    – John Taylor
    Commented Apr 24, 2023 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.