# Loop optimization

I wrote this to generate a 6D distribution of particles.

Bunch = {};
Do[
Particle = {0, 0, 0, 0, 0, 0};
tempX = { RandomVariate[NormalDistribution[0, Sqrt[EmitX]]],
RandomVariate[NormalDistribution[0, Sqrt[EmitX]]]};
tempY = { RandomVariate[NormalDistribution[0, Sqrt[EmitY]]],
RandomVariate[NormalDistribution[0, Sqrt[EmitY]]]};
coordS = {PhiS + \[Pi], BucketHeightP};
While[Abs[coordS[[1]]] >= (PhiS + \[Pi]) ||
Abs[coordS[[2]]] >= \[Delta][coordS[[1]]],
RandomVariate[NormalDistribution[0, MAXDPP ]]};];
coordX = {XCO + Sqrt[BetaX] tempX[[1]] + Dx coordS[[2]],
XPCO + (-(AlphaX/Sqrt[BetaX]) tempX[[1]] +
1/Sqrt[BetaX] tempX[[2]]) + Dpx coordS[[2]]};
coordY = {YCO + Sqrt[BetaY] tempY[[1]],
YPCO - AlphaY/Sqrt[BetaY] tempY[[1]] + 1/Sqrt[BetaY] tempY[[2]]};
Particle = Flatten[{coordX, coordY, coordS}];
Bunch = Append[Bunch, Particle],
{i, 1, 20000}];


I don't think the details are important. I call 6 times RandomVariate which I'm afraid is quite inefficient. Would you suggest some major improvement ?

-

Your best bet is to remove the procedural programming Do, While and Append Statements, building lists with Append is not quick. Then embrace a functional programming approach on which Mathematica thrives. Making use of constructs like Transpose, Part, Nest, Map, Table and Fold. These are generally much faster and lead to eventually to less buggy code.

If machine precision numbers are enough you can look at using Compile.

There are many approaches, but you could try something along these lines:-

There is no need to initialise Particle in Mathematica.

You can eliminate the Do loop and its 80,000 thousand calls to RandomVariate by just two calls.

dimensions=2; numParticles=20000;
XYs =
Transpose[{RandomVariate[NormalDistribution[0, Sqrt[EmitX]], {numParticles,dimensions}],
RandomVariate[NormalDistribution[0, Sqrt[EmitY]], {numParticles,dimensions}]}];


This gives you a list of pairs of Xs and Ys of length numParticles.

You could then write a set of functions to compute coordS, coordsX, coordsY, say C[], cX[], cY[].

Then combine Table and NestWhile.

The basic form for NestWhile is

NestWhile[function, initial state, terminating condition]


The terminating condition is checked, 'function' is is applied to 'initial state' producing a result.

The result is held in a variable known as #. # may have several components. These components are accessed by the part notation, #[[1]], #[[2]] and so on.

The terminating condition is then checked and if true, result is then given to function to evaluate.

This produces a new result ... and so on. It's basically recursion. The trick is to make the 'initial condition' have the same structure as the result of applying 'function'. Function may actually be a compound list of results gained by using several functions. And sometimes you need to 'pad' an element of the list, as it can be used to hold the parameters for the nest.

In your case perhaps something to compute the {coordX, coordY, coordS, XY }.

The shape of the code would look something like this

res =
Table[
NestWhile[{C[],cX[#[[1,1],#[[3,2]]],Cy[#[[2]]]}&,
{p[[1]],p[[2]],{PhiS + Pi, BucketHeightP}},
Abs[#[[3,1]]] >= (PhiS + Pi) ||  Abs[#[[3,2]]] >= d[#[[3,1]]] &],
{p,XYs}]


Assuming PhiS BucketHeightP do not change between iterations.

This is just a roughed out solution and probably skips over some details and could definitely be improved but should help you towards a faster implementation and give you a flavour of the functional approach.

I hope you can see that the form of the solution is essentially a Table command with a NestWhileand 3 auxiliary functions. You could even insert the initial creation of the variable XYs into Table command, but that would have just reduced comprehension in this case.

This is typical of what functional programming brings to the party. Very concise programs with a high level of abstraction.

One way of looking at it is that you start with your initial data and then apply some functional transforms to arrive at your final data. This avoids a lot of the bugs introduced by procedural programs that have lots of state information held in global variables. State information that often gets inadvertently corrupted.

Apologies, but this isn't debugged code, but a template solution - so there might be a few niggles to iron out before it flies.

Further performance tuning tips can be found at: tuning

-

Minor improvement: RandomVariate accepts a second argument with the number of elements you want to create. So your assignments to tempX, etc, are equivalent to

tempX = RandomVariate[NormalDistribution[0, Sqrt[EmitX]], 2];
tempY = RandomVariate[NormalDistribution[0, Sqrt[EmitY]], 2];

-
I did that, thanks ! –  Cedric H. Mar 28 '12 at 9:59

As you expect, it's much faster to call RandomVariate fewer times. In fact, it'd be much much faster to call it only once at the beginning, generating large arrays and then operating on them:

In[20]:= Bunch = {};
Timing[Do[
Particle = {RandomVariate[NormalDistribution[0, 1]],
RandomVariate[NormalDistribution[0, 1]],
RandomVariate[NormalDistribution[0, 1]],
RandomVariate[NormalDistribution[0, 1]],
RandomVariate[NormalDistribution[0, 1]],
RandomVariate[NormalDistribution[0, 1]]};
Append[Bunch, Particle],
{i, 1, 200000}]]
Out[21]= {2.92156, Null}

In[22]:= Bunch = {};
Timing[Do[
Particle = RandomVariate[NormalDistribution[0, 1], 6];
Append[Bunch, Particle],
{i, 1, 200000}]]
Out[23]= {0.873262, Null}

In[24]:= Timing[Bunch = RandomVariate[NormalDistribution[0, 1], {20000, 6}];]
Out[25]= {0.006057, Null}


In your example, it's slightly different because you're not calling it a predefined number of times due to your While loop. However, this optimization can still be done for then tempX and tempY variables, which can become fixed-sized arrays generated by one or two calls to RandomVariate before the loop.

After that, further optimization depends on how many iterations the While loop does on average. If it's rarely iterated (i.e. in most cases the condition is true after the first iteration), then you can perform your calculations on a big array, then go over it and only “fix” the values that need fixing. Anything you can put in array form, outside a Do loop, will go much faster.

-

In addition to the other answers, it's generally not very efficient to create a table using Append repeatedly. This is because every time you call p = Append[p, exp], Mathematica creates a new copy of the list. A more efficient way would be to use Table instead, i.e.

Bunch = Table[ code; Particle, {i, 20000}];

-
I did that too ! –  Cedric H. Mar 28 '12 at 9:59