2
$\begingroup$

I am writing a simulation where I have an array with ten columns and on the order of millions of rows, for which I need to iterate on the order of ten thousand times, making computation time a major issue. At the moment I am using the Map[] function to apply calculations on every row (every row is independent of each other), but I realized that using the listability of the basic functions is significantly faster.

However, my problem is that I have one piecewise function that is not listable. I also cannot define a function outside of the calculations (as f[x_]:=...), since that incurs a huge bottleneck (factor ~20 slower computations). I am therefore looking for any tips or advice on how to solve this.

As an example, this is one of the calculations that are performed:

Q=0.18;
Qfactor=1.*10^-10;
Qspace=0.135;
sigma2289=7.83717;
Qspace17sigma2=0.0292835;
sigma17=3.7995;
sigma=1.64676;
TWOPI=2.*Pi;

Map[TWOPI*(Q + #[[7]]*Qfactor - 
  If[#[[5]]^2 + #[[6]]^2 < sigma2289, 
   Qspace17sigma2*Sqrt[#[[5]]^2 + #[[6]]^2], 
   Qspace/(Sqrt[#[[5]]^2 + #[[6]]^2] - sigma17)/sigma]) &, 
   particleArray];

My issue is this part:

 If[#[[5]]^2 + #[[6]]^2 < sigma2289, 
   Qspace17sigma2*Sqrt[#[[5]]^2 + #[[6]]^2], 
   Qspace/(Sqrt[#[[5]]^2 + #[[6]]^2] - sigma17)/sigma]

The conditional part of that function appears to not be listable. However, if I remove this and compare the computational times of the Map[] version, versus the Listability version, you can see that there is a significant improvement:

in: AbsoluteTiming[particleArray[[All, 10]] = 
       Map[TWOPI*(Q + #[[3]]*Qfactor - 
         Qspace17sigma2*Sqrt[#[[1]]^2 + #[[2]]^2]) &, particleArray];]
out: (0.14621, Null)

in: AbsoluteTiming[particleArray[[All, 8]] = 
     TWOPI*(Q + particleArray[[All, 3]]*Qfactor - 
     Qspace17sigma2*Sqrt[particleArray[[All, 1]]^2 + particleArray[[All, 2]]^2]);]
out: (0.01042, Null)

You can generate particleArray using the following code, for testing:

heavyGaussian = MixtureDistribution[{0.75, 0.25},                   
  {MultinormalDistribution[{0, 0}, {{1, 0}, {0, 1}}], 
   MultinormalDistribution[{0, 0}, {{1.8, 0}, {0, 1.8}}]}];

particleArray = ParallelTable[
   Flatten[{RandomVariate[heavyGaussian],
     RandomVariate[
       MultinormalDistribution[{0, 3.14159}, {{9.*10^-8,0},{0,0.0194882}}]],
 0, 0, 0, 0, 0, 0}], {x, 1, 100000}];
particleArray[[All, 3]] *= 2.6*10^10;

I also have CUDA enabled and a decent GPU, so any tips related to that would also be appreciated.

$\endgroup$
  • 1
    $\begingroup$ There’s a piecewise-to-unitstep converter in the Simplify context — something like PWToUnitStep[]. You can use PiecewiseExpand[] to convert If[] to Piecewise[] first if needed. $\endgroup$ – Michael E2 Feb 8 at 22:09
5
$\begingroup$

You may use Compile as follows:

Q = 0.18;
Qfactor = 1.*10^-10;
Qspace = 0.135;
sigma2289 = 7.83717;
Qspace17sigma2 = 0.0292835;
sigma17 = 3.7995;
sigma = 1.64676;
TWOPI = 2.*Pi;

cf = With[{
   TWOPI = TWOPI, Q = Q, Qfactor = Qfactor, Qspace = Qspace, 
   sigma = sigma, sigma17 = sigma17, sigma2289 = sigma2289, 
   Qspace17sigma2 = Qspace17sigma2},
  Compile[{{X, _Real, 1}},
   TWOPI*(Q + Compile`GetElement[X, 7]*Qfactor - If[
       Compile`GetElement[X, 5]^2 + Compile`GetElement[X, 6]^2 < 
        sigma2289, 
       Qspace17sigma2 Sqrt[
         Compile`GetElement[X, 5]^2 + Compile`GetElement[X, 6]^2], 
       Qspace/(Sqrt[
            Compile`GetElement[X, 5]^2 + Compile`GetElement[X, 6]^2] -
            sigma17)/sigma
       ]
     ),
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ]
  ]

And here is a speed comparison. First we have to generate particleArray. Notice how I generate it: I use RandomVariate with a second argument to produce many random numbers at once. Moreover, I use ConstantArray[0., {n, 6}] to generate all the zero (as machine precision 0., not as eaxct 0) and merge everything with Join[#,2]&. This way, particleArray is a packed array an can also be processed faster, in general. (The array woud be coerced to machine precision reals and packed when the CompiledFunction cf is applied to it anyways.)

n = 100000;
particleArray = Join[
     RandomVariate[heavyGaussian, n],
     RandomVariate[MultinormalDistribution[{0,3.14159}, {{9.*10^-8, 0}, {0, 0.0194882}}], n],
     ConstantArray[0., {n, 6}],
     2
     ]; // AbsoluteTiming // First

a = Map[
     TWOPI*(Q + #[[7]]*Qfactor - If[
          #[[5]]^2 + #[[6]]^2 < sigma2289, 
          Qspace17sigma2*Sqrt[#[[5]]^2 + #[[6]]^2], 
          Qspace/(Sqrt[#[[5]]^2 + #[[6]]^2] - sigma17)/sigma
          ]
        ) &, particleArray
     ]; // AbsoluteTiming // First

b = cf[particleArray]; // AbsoluteTiming // First
Max[Abs[a - b]]

0.020113

0.063061

0.003697

0.

$\endgroup$
  • $\begingroup$ This works perfectly, thank you! $\endgroup$ – bjorn Feb 8 at 20:48
  • 3
    $\begingroup$ You're welcome. Notice that I gave also some hints on how to speed up the random number generation by a factor of 600... $\endgroup$ – Henrik Schumacher Feb 8 at 20:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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