# Fast and Listable Piecewise function

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 + #[]*Qfactor -
If[#[]^2 + #[]^2 < sigma2289,
Qspace17sigma2*Sqrt[#[]^2 + #[]^2],
Qspace/(Sqrt[#[]^2 + #[]^2] - sigma17)/sigma]) &,
particleArray];


My issue is this part:

 If[#[]^2 + #[]^2 < sigma2289,
Qspace17sigma2*Sqrt[#[]^2 + #[]^2],
Qspace/(Sqrt[#[]^2 + #[]^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 + #[]*Qfactor -
Qspace17sigma2*Sqrt[#[]^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.

• 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. Feb 8 '19 at 22:09

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 + CompileGetElement[X, 7]*Qfactor - If[
CompileGetElement[X, 5]^2 + CompileGetElement[X, 6]^2 <
sigma2289,
Qspace17sigma2 Sqrt[
CompileGetElement[X, 5]^2 + CompileGetElement[X, 6]^2],
Qspace/(Sqrt[
CompileGetElement[X, 5]^2 + CompileGetElement[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 + #[]*Qfactor - If[
#[]^2 + #[]^2 < sigma2289,
Qspace17sigma2*Sqrt[#[]^2 + #[]^2],
Qspace/(Sqrt[#[]^2 + #[]^2] - sigma17)/sigma
]
) &, particleArray
]; // AbsoluteTiming // First

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

0.020113

0.063061

0.003697

0.

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