Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Here is the physical problem I would like to simulate:

Consider the following variations to the 1d random walk:

  1. allow for certain, random sites to be excluded with probability 1/2.

  2. allow for the walker to be able to hops either 1 or 2 units to left or right

  3. whether hopping 1 or 2 units is weighted by hopping rate

  4. time in between hops is taken from exponential distribution with width being total hopping rate (to all 4 possible sites)

In the end, I would like to find out how much slower a walker takes to reach, say x = 100, when the exclusions are present as to when they are not.

The following will be my attempt at the problem. I think I have the basics down but I am hoping there are areas to decrease the computation time. These problems were alluded to earlier in the question and answers of 1D Random Walk variant.

First the 4 different hop rates are defined:

kBack1[k0_] := k0 Exp[- 2 ];
kFor1[k0_] := k0 Exp[- 2 ];
kBack2[k0_] := k0 Exp[- 2  * 2];
kFor2[k0_] := k0 Exp[- 2 *2];

I know these are redundant since there are actually only two unique rates. However, I would like to eventually add a bias to the random walk which will make the 4 rates all different.

Now we need to choose which site the walker will hop to and how long it takes to make the hop. So the excluded sites are accounted for by If statements. If an excluded site is nearby, there is a probability that that site is unreachable. The new site (newCoords) is chosen by weighting the different hopping rates. The time-step is sampled from the exponential distribution and the simulation time is updated.

FindRate[k0_] := Module[{},
  kB1 = If[({Xi} \[Intersection] blockedSites + 1) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kBack1[k0]}], kBack1[k0]];
  kB2 = If[({Xi} \[Intersection] blockedSites + 2) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kBack2[k0]}], kBack2[k0]];
  kF1 = If[({Xi} \[Intersection] blockedSites - 1) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kFor1[k0]}], kFor1[k0]];
  kF2 = If[({Xi} \[Intersection] blockedSites - 2) != {}, 
RandomChoice[{1/2, 1/2} -> {0, kFor2[k0]}], kFor2[k0]];

newCoords = 
RandomChoice[{kB1, kB2, kF1, kF2} -> {Xi - 1, Xi - 2, Xi + 1, 
  Xi + 2}];

dt = RandomReal[ExponentialDistribution[kB1 + kB2 + kF1 + kF2]];
t = simT;
{Xi, simT} = {newCoords, t + dt}
]

EDIT: More explanation on how the hop is chosen: Consider the walker to be at point 4. The next hop is chosen in the following way. Is the walker 1 or 2 spaces away from an excluded point? Let's say point 6 is excluded. The rate (kF2) to point 6 is then 0. The other points (2,3, and 5) can be hopped to with respective rate kB2,kB1, and kF1. Which one does it hop to? Well, compare the different rates. In the example we'll have kB2, kB1 = kF1, and kF2 = 0. Which rate is chosen is decided by RandomChoice where the choices are weighted by the rates. END EDIT

Now it's time to string the hops together

hopMod[k0_, Xf_] := Module[{},
Xi = 0; simT = 0;
blockedSites = RandomInteger[{-100, Xf-2}, 11];(*keep exclusion away from endpoint*)
NestWhileList[FindRate[k0] &, {0, 0}, Xi < Xf &]
]

where I define a random set of excluding sites. I would prefer to not have adjacent or next-adjacent points not both be exclusions (so separated by at least 2 normal sites). Is there a way to force RandomInteger to not sample consecutive integers?

Now calculate the time to reach a specified endpoint, Xf. And average this for n configurations.

With[{k0 = 1, Xf = 10, n = 20},(* n is # of runs to average over*)
 Print[AbsoluteTiming[
  Table[p[i] = hopMod[k0, Xf], {i, 1, n}];]];
 iTable = Flatten[Table[p[j][[-1]][[1]]/p[j][[-1]][[2]], {j, 1, n}]]]

This took me 1.33 seconds. The time is highly varying. The time increases dramatically as Xf or increase. For some reason, when I used ParallelTable, I got identical results every time I calculated another iTable. Must be something with the seed going wrong. To get convergence, n needs to be increased.

So how can this be sped up? Some different schemes for the hopMod module were discussed at 1D Random Walk variant. However I was not sure how to generalized all those approaches to my situation where I also need to keep track of the simulation time.

share|improve this question
1  
Could you clarify more in words, given that you are in a certain place, say 4, how do you get to the next step? I am sure it can be understood from the code, but it would be a lighter reading to understand it from words –  Rojo Jan 22 '13 at 22:34
    
@Rojo done. I hope that helps –  BeauGeste Jan 22 '13 at 22:53
    
Way better, thanks. So, as parameters, you get the weights and the exponential distribution parameter for the times. And where it hops is independent of how long it took. And you want to hop until you reach certain value, and receive as output the list of pairs {time, position}? –  Rojo Jan 22 '13 at 23:11
    
No, sorry, the exponential parameter is always the sum of the weights? –  Rojo Jan 22 '13 at 23:13
1  
I don't understand why t is localized in your Module and all the other variables assigned there are global. In fact, I don't think t is needed at all. Why not just return {newCoords, simT + dt} and make the assignment {Xi, simT} = FindRate[...] at the point of call? –  m_goldberg Jan 23 '13 at 1:33
show 4 more comments

1 Answer

Here are some tweaks, plus a way to generate non-adjacent blocked sites. (I wasn't sure whether you wanted two free sites between blocked sites or just one. There's just one below, and comment where it could be changed to two.) The run time (per hop) is shorter by about 2/3, only a modest improvement.

kBack1[k0_] := k0 Exp[-2];
kFor1[k0_] := k0 Exp[-2];
kBack2[k0_] := k0 Exp[-2*2];
kFor2[k0_] := k0 Exp[-2*2];

FindRate[k0_, {Xi_, simT_}] :=
   ({kB1, kB2, kF1, kF2} =
      {kBack1[k0], kBack2[k0], kFor1[k0], kFor2[k0]} ({Xi - 1, Xi - 2, Xi + 1, Xi + 2} /. blockedSites);
    {RandomChoice[{kB1, kB2, kF1, kF2} -> {Xi - 1, Xi - 2, Xi + 1, Xi + 2}], 
     simT + RandomReal[ExponentialDistribution[kB1 + kB2 + kF1 + kF2]]});

generateBlockedSites[min_, max_, nSites_] := Reap[Nest[
     DeleteCases[#, 
       With[{site = RandomChoice[#]}, Sow[site :> RandomInteger[]]; 
        Alternatives @@ Range[site - 1, site + 1]]] &, (* Or Range[site - 2, site + 2] ?? *)
     Range[min, max],
     nSites
     ]; Sow[_Integer -> 1]][[-1, -1]];

hopMod[k0_, Xf_] := Module[{}, Xi = 0; simT = 0;
   blockedSites = generateBlockedSites[-100, Xf - 2, 11];(*keep exclusion away from endpoint*)
   NestWhileList[FindRate[k0, #] &, {0, 0}, #[[1]] < Xf &]];

Clear[p];
With[{k0 = 1, Xf = 10, n = 20},(*n is # of runs to average over*)
 Print[AbsoluteTiming[Table[p[i] = hopMod[k0, Xf], {i, 1, n}];]];
 iTable = Table[p[j][[-1]][[1]]/p[j][[-1]][[2]], {j, 1, n}]] (* Flatten unnecessary *)

The variable blockedSites is a list of replacement rules of the form

In[216]:= blockedSites

Out[216]= {-68 :> RandomInteger[], -24 :> RandomInteger[], -62 :> RandomInteger[],
  -88 :> RandomInteger[], -82 :> RandomInteger[], -39 :> RandomInteger[],
  -53 :> RandomInteger[], -27 :> RandomInteger[], -65 :> RandomInteger[],
  -60 :> RandomInteger[], -97 :> RandomInteger[], _Integer -> 1}

Basically it takes a list of sites {Xi - 1, Xi - 2, Xi + 1, Xi + 2} and changes them to 1s or 0s, according to whether the current site can hop to the corresponding sites. This is then multiplied by the rates {kBack1[k0], kBack2[k0], kFor1[k0], kFor2[k0]} to get the rates {kB1, kB2, kF1, kF2} for the hop. I hope that preserves the model for the simulation.

share|improve this answer
    
thanks! one small point is that in FindRate, I think I'll want simT + RandomReal[ExponentialDistribution... in the last line so that the total time outputs from hopMod instead of just the time difference. As far as speed, I think the most gain will come from changing NestWhileLoop to something else (see link in original question). –  BeauGeste Jan 25 '13 at 16:04
    
@BeauGeste Fixed. Makes sense, too. Somehow I missed that I had dropped that part. –  Michael E2 Jan 25 '13 at 19:56
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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