3
$\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
  • 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, 2019 at 22:09

1 Answer 1

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$
2
  • $\begingroup$ This works perfectly, thank you! $\endgroup$
    – a20
    Feb 8, 2019 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$ Feb 8, 2019 at 20:54

Your Answer

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

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