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.