Improving the run time of random walk simulation

Question

One of our professors told us that during his thesis he used Mathematica in simulating it since it is about networks, I don't know how he did it but all I know is he said it took about two weeks to run it and get a result. We use it on our research as well and we're new to it of course, one of our goal is to determine the estimated time or mean of two points to meet. The meet-up is defined position1-position2=1. It is repeated 10,000 times.

This code:

For[i = 1, i <= 10000, i++,
pos1 = 2;
pos2 = -2;
t = 1; While[Abs[pos1 - pos2] > 1,
If[2 > pos1 > -2, pos1 = pos1 + RandomChoice[{-1, 0, 1}], 
pos1 = pos1 + RandomChoice[{-1, 0}]];
If[2 > pos2 > -2, pos2 = pos2 + RandomChoice[{-1, 0, 1}], 
pos2 = pos2 + RandomChoice[{0, 1}]]; t++];
time[i] = t]
avetime = N[Mean[Table[time[i], {i, 1, 10000}]], 5]
stdev = N[StandardDeviation[Table[time[i], {i, 1, 10000}]]]

To simply explain this, what it does is there are two points on a one-dimensional integer line, on -2 and 2. What it does is move until they meet(not joined) as defined in our meet-up position1-position2=1, as put as a test inside the while code, "Abs[pos1 - pos2] > 1" it stops once this is false which is when it becomes pos1-pos2=1. That is done 10,000 times(that is why we put it inside the 'For' code). Now this actually ran with ease and fast but when we tried -40 to 40 and -50 to 50 it started to take a long time to get a result.

Now we are trying -500 and 500 and it's been 12 hours and still got no result. Do I wait like my professor turned his computer on for about 2 weeks to get a result or is the code just inefficient? It seems simple to take a long time to run. I'm using a laptop that is not high-end, does that matter? Also if I make my laptop sleep and the screen locks, when I turn it back on does it stop the running Mathematica and resume it when I turn it on again? If the code is inefficient any comments, suggestions, or improvements are open just bear that we don't necessarily need a graphic, getting the mean is enough.

 For[i = 1, i <= 10000, i++,
 pos1 = 500;
 pos2 = -500;
 t = 1; While[Abs[pos1 - pos2] > 1,
 If[500 > pos1 > -500, pos1 = pos1 + RandomChoice[{-1, 0, 1}], 
 pos1 = pos1 + RandomChoice[{-1, 0}]];
 If[500 > pos2 > -500, pos2 = pos2 + RandomChoice[{-1, 0, 1}], 
 pos2 = pos2 + RandomChoice[{0, 1}]]; t++];
 time[i] = t]
 avetime = N[Mean[Table[time[i], {i, 1, 10000}]], 5]
 stdev = N[StandardDeviation[Table[time[i], {i, 1, 10000}]]]

Answer 1

Expanding upon my comment, I would recommend using Compile:

First we define a CompiledFunction that takes the bounds as its argument

cf = Compile[{{n, _Integer}},
Module[
{t, pos1, pos2},
pos1 = n;
pos2 = -n;
t = 1;
    While[
     Abs[pos1 - pos2] > 1,
     If[500 > pos1 > -500,
      pos1 += RandomChoice[{-1, 0, 1}],
      pos1 += RandomChoice[{-1, 0}]
     ];
    If[
     500 > pos2 > -500,
     pos2 += RandomChoice[{-1, 0, 1}],
     pos2 += RandomChoice[{0, 1}]
    ];
    t++
   ];
t
]];

Then single calls to cf with a given bound tell us how many steps this iteration takes.

In[139]:= cf[500]

Out[139]= 186021

Then we can use Table rather than For to make repeated calls.

results = Table[cf[500], 100];

Now you can do whatever computations you would like to do with just some minor modifications:

avetime = N[Mean[results], 5]
stdev = N[StandardDeviation[results]]

Yields

4.2925*10^5

and

292485.

This seems to take about 10 seconds on my machine for the first 100 tests, so it seems increasing to 10,000 iterations should probably be faster.

You can play around with options like CompilationTarget to speed up evaluation.

Good luck!

Edit --

In fact, changing CompilationTarget->"C" increased speed by another 40% or so. You can follow the docs to see how to get this to work on your machine if it does not work natively.

Answer 2

There compiling the code is already an essential step. But there is much more room for fine-tuning:

• The checks for collision with the boundaries can be simplified a bit (the right particle cannot crash into the left wall before crashing into the left particle).

• We can make the function Listable and run it parallelized to process many simulations at once.

• We can turn off some unneeded security checks with RuntimeOptions -> "Speed".

The biggest effect however, has the observation that the calls to RandomChoice (as almost all calls to random number generators) has a huge calling overhead. This calling overhead can be reduced significantly by generating many random numbers at once and by putting them into a reservoir from which the numbers can be drawn consecutively. You can find an implementation of that below. Notice also that I replaced RandomChoice by RandomInteger because the latter performs slightly better. (Quite likely, RandomChoice effectively calls RandomInteger and additionally performs some replace or read operations.)

Block[{irand, nrand, rand, iboundaryrand, nboundaryrand, boundaryrand},
   (* Prepare the code for random number generation; this will be inlined into the code below. *)
   With[{
     newrand = If[irand < nrand,
       irand++;
       Compile`GetElement[rand, irand]
       ,
       rand = RandomInteger[{-1, 1}, {nrand}];
       irand = 1;
       Compile`GetElement[rand, irand]
       ],
     newboundaryrand = If[iboundaryrand < nboundaryrand,
       iboundaryrand++;
       Compile`GetElement[boundaryrand, iboundaryrand]
       ,
       boundaryrand = RandomInteger[{0, 1}, {nboundaryrand}];
       iboundaryrand = 1;
       Compile`GetElement[boundaryrand, iboundaryrand]
       ]
     },


    cf2 = Compile[{{n, _Integer}},
     Block[{t, pos1, pos2, irand, nrand, rand, iboundaryrand, nboundaryrand, boundaryrand},

      (* Initialize random number reservoirs. *)
      nrand = 10 n;
      rand = RandomInteger[{-1, 1}, {nrand}];
      irand = 1;
      nboundaryrand = n;
      boundaryrand = RandomInteger[{0, 1}, {nboundaryrand}];
      iboundaryrand = 1;

      (* Initialize simulation. *)
      pos1 = n;
      pos2 = -n;
      t = 1;

      (* Run simulation. *)
      While[
       Abs[pos1 - pos2] > 1,
       If[pos1 < n,
        pos1 += newrand,
        pos1 -= newboundaryrand
        ];
       If[pos2 > -n,
        pos2 += newrand,
        pos2 += newboundaryrand
        ];
       t++
       ];
      t
      ],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]


    ]
   ];

This is how this function performs versus enano9314's function cf (compiled with CompilationTarget->"C") on my Haswell Quad Core CPU:

time = AbsoluteTiming[results = Table[cf[500], {100}]][[1]];
time2 = AbsoluteTiming[results2 = Table[cf2[500], {100}]][[1]];
time2p = AbsoluteTiming[results2p = cf2[ConstantArray[500, {100}]]][[1]];
{Mean[N@results], time, time/Mean[N@results]}
{Mean[N@results2], time2, time2/Mean[N@results2]}
{Mean[N@results2p], time2p, time2p/Mean[N@results2p]}

{467725., 10.7932, 0.000023076}

{439379., 0.47879, 1.0897*10^-6}

{428576., 0.145688, 3.39935*10^-7}

So the simulation time per iteration has decreased by a factor of 67.9. Running the simulation 10000 times in the interval [-500,500] takes about 14 seconds. So no need for letting it run for 14 days...

PS.: Notice that the parallelization of Listable CompiledFunctions will take only effect, if it is called correctly. In our example,

Table[cf2[500], {100}]]

will be executed serially while

cf2[ConstantArray[500, {100}]]]

will be executed parallelized.

Answer 3

If you take into account that the distance between pos1 and pos2 is the minimum number of steps before meeting, you can reduce both the number of calls to RandomChoice[ ] and the number of While[ ] iterations and If[ ] calls. I will just adapt the function posted by enano9314, changing RandomChoice[ ] by RandomInteger[ ], as suggested by Schumacher.

cfsi[n_] := Module[{pos1, pos2, d, stride},
pos1 = n;
pos2 = -n;
t = 1;
While[(d = pos1 - pos2) > 1,
    stride = IntegerPart[d/2];
    If[pos1 + stride <= n, pos1 += Total[RandomInteger[{-1, 1}, stride] ],
       pos1 += Total[ Table[ If[n > pos1, RandomInteger[{-1, 1}], RandomInteger[{-1, 0}]],
                               {stride}] ] ];
    If[pos2 - stride >= -n, 
       pos2 += Total[ RandomInteger[{-1, 1}, stride] ], 
       pos2 += Total[ Table[ If[pos2 > -n, RandomInteger[{-1, 1}], RandomInteger[{0, 1}]], 
                               {stride}] ] ];
    t += stride];
t]

test = AbsoluteTiming[ Table[ cfsi[500] , {100}] ]; test[[1]]
7.47501
{Mean[test[[2]]] , StandardDeviation[test[[2]] ]} // N
295171., 219751.}

RandomChoice[ ] increases the procssing time in about 50 percent. The next gain in performance relies on that, when pos1+stride>n, the point is likely far from the wall, hence the stride can be divided in bits that don't cross the wall (bit=n-pos1). This code is more complicated to implement, but, since the upgraded code can be tested against a working code, it will not be so difficult to get it right.

I would like to call your attention also to the meaning of varname[ ], which is a form mainly reserved to functions. Mathematica allows you to give function values for specific arguments with the syntax you used: time[i]=t. Try ?time in your code to see that you defined 10000 specific values for time, and not a list with 10000 values. Elements of a list are accessed through time[[ ]], and elements can be added with AppendTo[time, t], or it can be build using Reap[ expression; Sow[t] ]